diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp
index 0e00a0e700..7875ddd887 100644
--- a/opm/grid/CpGrid.hpp
+++ b/opm/grid/CpGrid.hpp
@@ -1743,7 +1743,7 @@ namespace Dune
         /// \brief Get the Position of a vertex.
         /// \param cell The index identifying the cell.
         /// \return The coordinates of the vertex.
-        const Vector& vertexPosition(int vertex) const;
+        Vector vertexPosition(int vertex) const;
 
         /// \brief Get the area of a face.
         /// \param cell The index identifying the face.
@@ -1751,12 +1751,12 @@ namespace Dune
 
         /// \brief Get the coordinates of the center of a face.
         /// \param cell The index identifying the face.
-        const Vector& faceCentroid(int face) const;
+        Vector faceCentroid(int face) const;
 
         /// \brief Get the unit normal of a face.
         /// \param cell The index identifying the face.
         /// \see faceCell
-        const Vector& faceNormal(int face) const;
+        Vector faceNormal(int face) const;
 
         /// \brief Get the volume of the cell.
         /// \param cell The index identifying the cell.
@@ -1764,7 +1764,7 @@ namespace Dune
 
         /// \brief Get the coordinates of the center of a cell.
         /// \param cell The index identifying the face.
-        const Vector& cellCentroid(int cell) const;
+        Vector cellCentroid(int cell) const;
 
         /// \brief An iterator over the centroids of the geometry of the entities.
         /// \tparam codim The co-dimension of the entities.
diff --git a/opm/grid/cpgrid/CpGrid.cpp b/opm/grid/cpgrid/CpGrid.cpp
index 8eaac4e59e..f980966834 100644
--- a/opm/grid/cpgrid/CpGrid.cpp
+++ b/opm/grid/cpgrid/CpGrid.cpp
@@ -1755,7 +1755,7 @@ const Dune::FieldVector<double,3> CpGrid::faceAreaNormalEcl(int face) const
     }
 }
 
-const Dune::FieldVector<double,3>& CpGrid::vertexPosition(int vertex) const
+Dune::FieldVector<double,3> CpGrid::vertexPosition(int vertex) const
 {
     return current_view_data_->geomVector<3>()[cpgrid::EntityRep<3>(vertex, true)].center();
 }
@@ -1765,12 +1765,12 @@ double CpGrid::faceArea(int face) const
     return current_view_data_->geomVector<1>()[cpgrid::EntityRep<1>(face, true)].volume();
 }
 
-const Dune::FieldVector<double,3>& CpGrid::faceCentroid(int face) const
+Dune::FieldVector<double,3> CpGrid::faceCentroid(int face) const
 {
     return current_view_data_->geomVector<1>()[cpgrid::EntityRep<1>(face, true)].center();
 }
 
-const Dune::FieldVector<double,3>& CpGrid::faceNormal(int face) const
+Dune::FieldVector<double,3> CpGrid::faceNormal(int face) const
 {
     return current_view_data_->face_normals_.get(face);
 }
@@ -1780,7 +1780,7 @@ double CpGrid::cellVolume(int cell) const
     return current_view_data_->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].volume();
 }
 
-const Dune::FieldVector<double,3>& CpGrid::cellCentroid(int cell) const
+Dune::FieldVector<double,3> CpGrid::cellCentroid(int cell) const
 {
     return current_view_data_->geomVector<0>()[cpgrid::EntityRep<0>(cell, true)].center();
 }
@@ -2857,10 +2857,10 @@ void CpGrid::refineAndProvideMarkedRefinedRelations( /* Marked elements paramete
                                                      int& cell_count,
                                                      std::vector<std::vector<int>>& preAdapt_level_to_leaf_cells_vec,
                                                      /* Additional parameters */
-                                                     const std::vector<std::array<int,3>>& cells_per_dim_vec) const 
+                                                     const std::vector<std::array<int,3>>& cells_per_dim_vec) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     // Each marked element for refinement (mark equal to 1), will be refined individuality, creating its own Lgr. The element index will
     // be also used to identify its lgr. Even though, in the end, all the refined entities will belong to a unique level grid.
     // For this reason, we associate "-1" with those elements that are not involved in any refinement and will appear
@@ -2878,7 +2878,7 @@ void CpGrid::refineAndProvideMarkedRefinedRelations( /* Marked elements paramete
             cell_count +=1;
             preAdapt_level_to_leaf_cells_vec[element.level()][element.getLevelElem().index()] = cell_count;
         }
-        
+
         // When the element is marked for refinement, we also mark its corners and faces
         // since they will get replaced by refined ones.
         if (getMark(element) ==  1) {
@@ -2914,7 +2914,7 @@ void CpGrid::refineAndProvideMarkedRefinedRelations( /* Marked elements paramete
                 refined_cell_count_vec[shiftedLevel] +=1;
 
             }
-            
+
             preAdapt_parent_to_children_cells_vec[element.level()][element.getLevelElem().index()] = std::make_pair( markedElemLevel, refinedChildrenList);
             for (const auto& [markedCorner, lgrEquivCorner] : parentCorners_to_equivalentRefinedCorners) {
                 cornerInMarkedElemWithEquivRefinedCorner[markedCorner].push_back({elemIdx, lgrEquivCorner});
@@ -2934,7 +2934,7 @@ CpGrid::defineChildToParentAndIdxInParentCell(const std::map<std::array<int,2>,s
                                               const int& cell_count) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     // ------------------------ Refined grid parameters
     // Refined child cells and their parents. Entry is {-1,-1} when cell has no father. Otherwise, {level parent cell, parent cell index}
     // Each entry represents a refined level.
@@ -3304,7 +3304,7 @@ void CpGrid::identifyLeafGridCorners(std::map<std::array<int,2>,int>& elemLgrAnd
                                      const std::vector<std::array<int,3>>& cells_per_dim_vec) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     // Step 1. Select/store the corners from the starting grid not involved in any (new) LGR.
     //         Replace the corners from level zero involved in LGR by the equivalent ones, born in LGRs.
     //         In this case, we avoid repetition considering the last appearance of the level zero corner
@@ -3458,8 +3458,8 @@ void CpGrid::identifyLeafGridFaces(std::map<std::array<int,2>,int>& elemLgrAndEl
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
 
     // Max level before calling adapt.
-    const int& preAdaptMaxLevel = this->maxLevel(); 
-    
+    const int& preAdaptMaxLevel = this->maxLevel();
+
     // Step 1. Add the LGR faces, for each LGR
     for (int elem = 0; elem < current_view_data_->size(0); ++elem) {
         if (markedElem_to_itsLgr[elem]!=nullptr)  {
@@ -3520,7 +3520,7 @@ void CpGrid::populateLeafGridCorners(Dune::cpgrid::EntityVariableBase<cpgrid::Ge
                                      const std::unordered_map<int,std::array<int,2>>& adaptedCorner_to_elemLgrAndElemLgrCorner) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     adapted_corners.resize(corner_count);
     for (int corner = 0; corner < corner_count; ++corner) {
         const auto& [elemLgr, elemLgrCorner] = adaptedCorner_to_elemLgrAndElemLgrCorner.at(corner);
@@ -3561,7 +3561,7 @@ void CpGrid::populateLeafGridFaces(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
                                    const int& preAdaptMaxLevel) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     adapted_faces.resize(face_count);
     mutable_face_tags.resize(face_count);
     mutable_face_normals.resize(face_count);
@@ -3741,7 +3741,7 @@ void CpGrid::populateLeafGridCells(Dune::cpgrid::EntityVariableBase<cpgrid::Geom
                                    const int& preAdaptMaxLevel) const
 {
     // If the (level zero) grid has been distributed, then the preAdaptGrid is data_[0]. Otherwise, preApaptGrid is current_view_data_.
-    
+
     // --- Adapted cells ---
     // Store the adapted cells. Main difficulty: to lookup correctly the indices of the corners and faces of each cell.
     adapted_cells.resize(cell_count);
@@ -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.
@@ -3871,7 +3871,7 @@ void CpGrid::populateRefinedCells(std::vector<Dune::cpgrid::EntityVariableBase<c
                                   const std::map<std::array<int,2>,int>& markedElemAndEquivRefinedCorn_to_corner,
                                   const std::vector<std::vector<std::array<int,2>>>& cornerInMarkedElemWithEquivRefinedCorner,
                                   const std::vector<std::array<int,3>>&  cells_per_dim_vec) const
-{   
+{
     // --- Refined cells ---
     for (std::size_t shiftedLevel = 0; shiftedLevel < refined_cell_count_vec.size(); ++shiftedLevel) {
 
@@ -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]);
@@ -4137,7 +4137,7 @@ void CpGrid::updateCornerHistoryLevels(const std::vector<std::vector<std::array<
             currentData()[refinedLevel]->corner_history_[refinedCorner] = preAdaptGrid_corner_history.empty() ? std::array<int,2>{{0, static_cast<int>(corner)}} :  preAdaptGrid_corner_history[corner];
         }
     }
-   
+
     // corner_history_ leaf grid view
     for ( int leafCorner = 0; leafCorner < corner_count; ++leafCorner){
         currentData().back()->corner_history_.resize(corner_count);
@@ -4664,5 +4664,3 @@ int CpGrid::replaceLgr1FaceIdxByLgr2FaceIdx(const std::array<int,3>& cells_per_d
 }
 
 } // namespace Dune
-
-
diff --git a/opm/grid/cpgrid/CpGridData.cpp b/opm/grid/cpgrid/CpGridData.cpp
index 4a2ec1a8e5..25df847a57 100644
--- a/opm/grid/cpgrid/CpGridData.cpp
+++ b/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}.
diff --git a/opm/grid/cpgrid/Entity.hpp b/opm/grid/cpgrid/Entity.hpp
index 169d5c5e63..deff9b7063 100644
--- a/opm/grid/cpgrid/Entity.hpp
+++ b/opm/grid/cpgrid/Entity.hpp
@@ -162,7 +162,7 @@ class Entity : public EntityRep<codim>
     }
 
     /// @brief Return the geometry of the entity (does not depend on its orientation).
-    const Geometry& geometry() const;
+    Geometry geometry() const;
 
     /// @brief Return the level of the entity in the grid hierarchy. Level = 0 represents the coarsest grid.
     int level() const;
@@ -387,7 +387,7 @@ unsigned int Entity<codim>::subEntities ( const unsigned int cc ) const
 }
 
 template <int codim>
-const typename Entity<codim>::Geometry& Entity<codim>::geometry() const
+typename Entity<codim>::Geometry Entity<codim>::geometry() const
 {
     return pgrid_->geomVector<codim>()[*this];
 }
@@ -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.");
diff --git a/opm/grid/cpgrid/Geometry.hpp b/opm/grid/cpgrid/Geometry.hpp
index ed21cbca1d..098cd941af 100644
--- a/opm/grid/cpgrid/Geometry.hpp
+++ b/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>
 
@@ -127,7 +129,7 @@ namespace Dune
             }
 
             /// Returns the position of the vertex.
-            const GlobalCoordinate& global(const LocalCoordinate&) const
+            GlobalCoordinate global(const LocalCoordinate&) const
             {
                 return pos_;
             }
@@ -172,7 +174,7 @@ namespace Dune
             }
 
             /// Returns the centroid of the geometry.
-            const GlobalCoordinate& center() const
+            GlobalCoordinate center() const
             {
                 return pos_;
             }
@@ -272,7 +274,7 @@ namespace Dune
             }
 
             /// This method is meaningless for singular geometries.
-            const GlobalCoordinate& global(const LocalCoordinate&) const
+            GlobalCoordinate global(const LocalCoordinate&) const
             {
                 OPM_THROW(std::runtime_error, "Geometry::global() meaningless on singular geometry.");
             }
@@ -319,20 +321,20 @@ namespace Dune
             }
 
             /// Returns the centroid of the geometry.
-            const GlobalCoordinate& center() const
+            GlobalCoordinate center() const
             {
                 return pos_;
             }
 
             /// This method is meaningless for singular geometries.
-            const FieldMatrix<ctype, mydimension, coorddimension>&
+            FieldMatrix<ctype, mydimension, coorddimension>
             jacobianTransposed(const LocalCoordinate& /* local */) const
             {
                 OPM_THROW(std::runtime_error, "Meaningless to call jacobianTransposed() on singular geometries.");
             }
 
             /// This method is meaningless for singular geometries.
-            const FieldMatrix<ctype, coorddimension, mydimension>&
+            FieldMatrix<ctype, coorddimension, mydimension>
             jacobianInverseTransposed(const LocalCoordinate& /*local*/) const
             {
                 OPM_THROW(std::runtime_error, "Meaningless to call jacobianInverseTransposed() on singular geometries.");
@@ -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.
@@ -526,7 +557,7 @@ namespace Dune
             }
 
             /// Returns the centroid of the geometry.
-            const GlobalCoordinate& center() const
+            GlobalCoordinate center() const
             {
                 return pos_;
             }
@@ -537,7 +568,7 @@ namespace Dune
             /// and {u_i} are the reference coordinates.
             /// g = g(u) = (g_1(u), g_2(u), g_3(u)), u=(u_1,u_2,u_3)
             /// g = map from (local) reference domain to global cell
-            const JacobianTransposed
+            JacobianTransposed
             jacobianTransposed(const LocalCoordinate& local_coord) const
             {
                 static_assert(mydimension == 3, "");
@@ -573,7 +604,7 @@ namespace Dune
             }
 
             /// @brief Inverse of Jacobian transposed. \see jacobianTransposed().
-            const JacobianInverseTransposed
+            JacobianInverseTransposed
             jacobianInverseTransposed(const LocalCoordinate& local_coord) const
             {
                 JacobianInverseTransposed Jti = jacobianTransposed(local_coord);
@@ -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,
diff --git a/opm/grid/cpgrid/GridHelpers.cpp b/opm/grid/cpgrid/GridHelpers.cpp
index 6f4f495599..a13fedbcf5 100644
--- a/opm/grid/cpgrid/GridHelpers.cpp
+++ b/opm/grid/cpgrid/GridHelpers.cpp
@@ -141,9 +141,9 @@ beginFaceCentroids(const Dune::CpGrid& grid)
     return FaceCentroidTraits<Dune::CpGrid>::IteratorType(grid, 0);
 }
 
-const double* cellCentroid(const Dune::CpGrid& grid, int cell_index)
+Vector cellCentroid(const Dune::CpGrid& grid, int cell_index)
 {
-    return &(grid.cellCentroid(cell_index)[0]);
+    return grid.cellCentroid(cell_index);
 }
 
 double cellVolume(const  Dune::CpGrid& grid, int cell_index)
@@ -161,8 +161,7 @@ CellVolumeIterator endCellVolumes(const Dune::CpGrid& grid)
     return CellVolumeIterator(grid, numCells(grid));
 }
 
-const FaceCentroidTraits<Dune::CpGrid>::ValueType&
-faceCentroid(const Dune::CpGrid& grid, int face_index)
+Vector faceCentroid(const Dune::CpGrid& grid, int face_index)
 {
     return grid.faceCentroid(face_index);
 }
@@ -184,14 +183,14 @@ face2Vertices(const Dune::CpGrid& grid)
     return Dune::cpgrid::FaceVerticesContainerProxy(&grid);
 }
 
-const double* vertexCoordinates(const Dune::CpGrid& grid, int index)
+Vector vertexCoordinates(const Dune::CpGrid& grid, int index)
 {
-    return &(grid.vertexPosition(index)[0]);
+    return grid.vertexPosition(index);
 }
 
-const double* faceNormal(const Dune::CpGrid& grid, int face_index)
+Vector faceNormal(const Dune::CpGrid& grid, int face_index)
 {
-    return &(grid.faceNormal(face_index)[0]);
+    return grid.faceNormal(face_index);
 }
 
 double faceArea(const Dune::CpGrid& grid, int face_index)
diff --git a/opm/grid/cpgrid/GridHelpers.hpp b/opm/grid/cpgrid/GridHelpers.hpp
index b910a06b8e..fc33f7662c 100644
--- a/opm/grid/cpgrid/GridHelpers.hpp
+++ b/opm/grid/cpgrid/GridHelpers.hpp
@@ -328,10 +328,10 @@ struct Cell2FacesTraits<Dune::CpGrid>
     typedef Dune::cpgrid::Cell2FacesContainer Type;
 };
 /// \brief An iterator over the cell volumes.
-template<const Dune::FieldVector<double, 3>& (Dune::CpGrid::*Method)(int)const>
+template<Dune::FieldVector<double, 3> (Dune::CpGrid::*Method)(int)const>
 class CpGridCentroidIterator
     : public Dune::RandomAccessIteratorFacade<CpGridCentroidIterator<Method>, Dune::FieldVector<double, 3>,
-                                              const Dune::FieldVector<double, 3>&, int>
+                                              Dune::FieldVector<double, 3>, int>
 {
 public:
     /// \brief Creates an iterator.
@@ -341,7 +341,7 @@ class CpGridCentroidIterator
         : grid_(&grid), cell_index_(cell_index)
     {}
 
-    const Dune::FieldVector<double, 3>& dereference() const
+    Dune::FieldVector<double, 3> dereference() const
     {
         return std::mem_fn(Method)(*grid_, cell_index_);
     }
@@ -349,7 +349,7 @@ class CpGridCentroidIterator
     {
         ++cell_index_;
     }
-    const Dune::FieldVector<double, 3>& elementAt(int n) const
+    Dune::FieldVector<double, 3> elementAt(int n) const
     {
         return  std::mem_fn(Method)(*grid_, n);
     }
@@ -430,7 +430,7 @@ double cellCentroidCoordinate(const Dune::CpGrid& grid, int cell_index,
 /// \brief Get the centroid of a cell.
 /// \param grid The grid whose cell centroid we query.
 /// \param cell_index The index of the corresponding cell.
-const double* cellCentroid(const Dune::CpGrid& grid, int cell_index);
+Vector cellCentroid(const Dune::CpGrid& grid, int cell_index);
 
 /// \brief Get vertical position of cell center ("zcorn" average).
 /// \brief grid The grid.
@@ -531,8 +531,7 @@ beginFaceCentroids(const Dune::CpGrid& grid);
 /// \param grid The grid.
 /// \param face_index The index of the specific face.
 /// \param coordinate The coordinate index.
-const FaceCentroidTraits<Dune::CpGrid>::ValueType&
-faceCentroid(const Dune::CpGrid& grid, int face_index);
+Vector faceCentroid(const Dune::CpGrid& grid, int face_index);
 
 template<>
 struct FaceCellTraits<Dune::CpGrid>
@@ -559,9 +558,9 @@ face2Vertices(const Dune::CpGrid& grid);
 /// \brief Get the coordinates of a vertex of the grid.
 /// \param grid The grid the vertex is part of.
 /// \param index The index identifying the vertex.
-const double* vertexCoordinates(const Dune::CpGrid& grid, int index);
+Vector vertexCoordinates(const Dune::CpGrid& grid, int index);
 
-const double* faceNormal(const Dune::CpGrid& grid, int face_index);
+Vector faceNormal(const Dune::CpGrid& grid, int face_index);
 
 double faceArea(const Dune::CpGrid& grid, int face_index);
 
diff --git a/opm/grid/cpgrid/Intersection.hpp b/opm/grid/cpgrid/Intersection.hpp
index e7d8376082..3c7290cbe6 100644
--- a/opm/grid/cpgrid/Intersection.hpp
+++ b/opm/grid/cpgrid/Intersection.hpp
@@ -161,7 +161,7 @@ namespace Dune
             /// @brief
             /// @todo Doc me!
             /// @return
-            const LocalGeometry& geometryInInside() const
+            LocalGeometry geometryInInside() const
             {
                 OPM_THROW(std::runtime_error, "This intersection class does not support geometryInInside().");
             }
@@ -171,7 +171,7 @@ namespace Dune
             /// @brief
             /// @todo Doc me!
             /// @return
-            const LocalGeometry& geometryInOutside() const
+            LocalGeometry geometryInOutside() const
             {
                 if (boundary()) {
                     OPM_THROW(std::runtime_error, "Cannot access geometryInOutside(), intersection is at a boundary.");
diff --git a/opm/grid/cpgrid/processEclipseFormat.cpp b/opm/grid/cpgrid/processEclipseFormat.cpp
index 4f7c51db04..a8e6373088 100644
--- a/opm/grid/cpgrid/processEclipseFormat.cpp
+++ b/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]);
             }
         };
 
diff --git a/opm/grid/transmissibility/TransTpfa_impl.hpp b/opm/grid/transmissibility/TransTpfa_impl.hpp
index 9b4f4b9c3d..c2cc1f9cb8 100644
--- a/opm/grid/transmissibility/TransTpfa_impl.hpp
+++ b/opm/grid/transmissibility/TransTpfa_impl.hpp
@@ -87,7 +87,7 @@ tpfa_htrans_compute(const Grid* G, const double *perm, double *htrans)
             s = 2.0*(face_cells(*f, 0) == c) - 1.0;
             n = faceNormal(*G, *f);
             const double* nn=multiplyFaceNormalWithArea(*G, *f, n);
-            const double* fc = &(faceCentroid(*G, *f)[0]);
+            const auto fc = faceCentroid(*G, *f);
             dgemv_("No Transpose", &nrows, &ncols,
                    &a1, K, &ldA, nn, &incx, &a2, &Kn[0], &incy);
             maybeFreeFaceNormal(*G, nn);
diff --git a/tests/cpgrid/geometry_test.cpp b/tests/cpgrid/geometry_test.cpp
index 74645353e0..d42d33cfc8 100644
--- a/tests/cpgrid/geometry_test.cpp
+++ b/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});