diff --git a/opm/grid/CpGrid.hpp b/opm/grid/CpGrid.hpp
index 74f18e6e0..eeaa89431 100644
--- a/opm/grid/CpGrid.hpp
+++ b/opm/grid/CpGrid.hpp
@@ -334,7 +334,7 @@ namespace Dune
         /// Set whether we want to have unique boundary ids.
         /// \param uids if true, each boundary intersection will have a unique boundary id.
         void setUniqueBoundaryIds(bool uids);
-       
+
 
         // --- Dune interface below ---
 
@@ -345,7 +345,7 @@ namespace Dune
         /// It's the same as the class name.
         /// What did you expect, something funny?
         std::string name() const;
-        
+
         /// Return maximum level defined in this grid. Levels are 0 and 1,  maxlevel = 1 (not counting leafview), 0 = the coarsest level.
         int maxLevel() const;
 
@@ -355,7 +355,7 @@ namespace Dune
         /// one past the end on this level
         template<int codim>
         typename Traits::template Codim<codim>::LevelIterator lend (int level) const;
-        
+
         /// Iterator to first entity of given codim on level and PartitionIteratorType
         template<int codim, PartitionIteratorType PiType>
         typename Traits::template Codim<codim>::template Partition<PiType>::LevelIterator lbegin (int level) const;
@@ -394,7 +394,7 @@ namespace Dune
 
         /// \brief Access to the LocalIdSet
         const Traits::LocalIdSet& localIdSet() const;
-        
+
         /// \brief Access to the LevelIndexSets
         const Traits::LevelIndexSet& levelIndexSet(int level) const;
 
@@ -626,10 +626,10 @@ namespace Dune
                                                                              const int& cell_count) const;
 
         /// @brief Define refined level grid cells indices and leaf grid view (or adapted grid) cells indices relations. Namely, level_to_leaf_cells_ for each new
-        ///        refined level grid, and leaf_to_level_cells_ for the updated leaf grid view. 
+        ///        refined level grid, and leaf_to_level_cells_ for the updated leaf grid view.
         ///
         /// @param [in] elemLgrAndElemLgrCell_to_refinedLevelAdRefinedCell:  Each marked element has been refined in its "own elemLgr". Refined entities should be stored in
-        ///                                                                  the corresponding assigned refined level grid. To keep track of the cell index relation, we 
+        ///                                                                  the corresponding assigned refined level grid. To keep track of the cell index relation, we
         ///                                                                  associate each
         ///                                                                  { marked element index ("elemLgr"), refined cell index in the auxiliary single-cell-refinement } with
         ///                                                                  { refined level grid assigned for the marked element, refined cell index in refined level grid }.
@@ -761,7 +761,7 @@ namespace Dune
         /// @param [out] adaptedFace_to_elemLgrAndElemLgrFace: Each marked element has been refined in its "own elemLgr". Refined corners should be stored in
         ///                                                    the corresponding assigned refined level grid. To keep track of the face index relation, we associate each
         ///                                                    face index in the leaf grid view (or adapted grid) with
-        ///                                                    { marked element index ("elemLgr"), refined corner index in the auxiliary single-cell-refinement }. 
+        ///                                                    { marked element index ("elemLgr"), refined corner index in the auxiliary single-cell-refinement }.
         /// @param [out] face_count:                           Total amount of faces on the leaf grid view (or adapted grid).
         /// @param [in] markedElem_to_itsLgr
         /// @param [in] assignRefinedLevel
@@ -814,7 +814,8 @@ namespace Dune
                                   const int& preAdaptMaxLevel,
                                   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;
+                                  const std::vector<std::array<int,3>>&  cells_per_dim_vec,
+                                  const std::vector<Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<0,3>>>& refined_corners_vec) const;
 
         /// @brief Set geometrical and topological attributes for each refined level grid.
         void setRefinedLevelGridsGeometries( /* Refined corner arguments */
@@ -1087,7 +1088,7 @@ namespace Dune
         int getParentFaceWhereNewRefinedFaceLiesOn(const std::array<int,3>& cells_per_dim, int faceIdxInLgr,
                                                    const std::shared_ptr<cpgrid::CpGridData>& elemLgr_ptr,
                                                    int elemLgr)  const;
-        
+
         /// --------------- Auxiliary methods to support Adaptivity (end) ---------------
 
         /// @brief Check if there are non neighboring connections on blocks of cells selected for refinement.
@@ -1712,7 +1713,7 @@ namespace Dune
         /// cell due the face being part of the grid boundary or the
         /// cell being stored on another process.
         int faceCell(int face, int local_index, int level = -1) const;
-      
+
         /// \brief Get the sum of all faces attached to all cells.
         ///
         /// Each face identified by a unique index is counted as often
@@ -1743,7 +1744,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 +1752,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 +1765,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 0945690bd..b7b828ffc 100644
--- a/opm/grid/cpgrid/CpGrid.cpp
+++ b/opm/grid/cpgrid/CpGrid.cpp
@@ -1761,7 +1761,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();
 }
@@ -1771,12 +1771,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);
 }
@@ -1786,7 +1786,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();
 }
@@ -2863,10 +2863,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
@@ -2884,7 +2884,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) {
@@ -2920,7 +2920,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});
@@ -2940,7 +2940,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.
@@ -3310,7 +3310,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
@@ -3464,8 +3464,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)  {
@@ -3526,7 +3526,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);
@@ -3567,7 +3567,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);
@@ -3747,7 +3747,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);
@@ -3851,7 +3851,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.
@@ -3876,8 +3876,9 @@ void CpGrid::populateRefinedCells(std::vector<Dune::cpgrid::EntityVariableBase<c
                                   const int& preAdaptMaxLevel,
                                   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
-{   
+                                  const std::vector<std::array<int,3>>&  cells_per_dim_vec,
+                                  const std::vector<Dune::cpgrid::EntityVariableBase<cpgrid::Geometry<0,3>>>& refined_corners_vec) const
+{
     // --- Refined cells ---
     for (std::size_t shiftedLevel = 0; shiftedLevel < refined_cell_count_vec.size(); ++shiftedLevel) {
 
@@ -3885,7 +3886,9 @@ void CpGrid::populateRefinedCells(std::vector<Dune::cpgrid::EntityVariableBase<c
         refined_cell_to_point_vec[shiftedLevel].resize(refined_cell_count_vec[shiftedLevel]);
         refined_global_cell_vec[shiftedLevel].resize(refined_cell_count_vec[shiftedLevel]);
 
-        const auto& allLevelCorners = refined_geometries_vec[shiftedLevel].geomVector(std::integral_constant<int,3>());
+        auto copyCorners = refined_corners_vec[shiftedLevel];
+        auto allLevelCorners = *refined_geometries_vec[shiftedLevel].geomVector(std::integral_constant<int,3>());
+        allLevelCorners.swap(copyCorners);
 
         for (int cell = 0; cell < refined_cell_count_vec[shiftedLevel]; ++cell) {
 
@@ -3972,7 +3975,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, indices_storage_ptr);
         } // refined_cells
         // Refined face to cell.
         refined_cell_to_face_vec[shiftedLevel].makeInverseRelation(refined_face_to_cell_vec[shiftedLevel]);
@@ -4048,7 +4051,8 @@ void CpGrid::setRefinedLevelGridsGeometries( /* Refined corner arguments */
                          preAdaptMaxLevel,
                          markedElemAndEquivRefinedCorn_to_corner,
                          cornerInMarkedElemWithEquivRefinedCorner,
-                         cells_per_dim_vec);
+                         cells_per_dim_vec,
+                         refined_corners_vec);
 
 }
 
@@ -4143,7 +4147,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);
@@ -4670,5 +4674,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 883d46a44..789a9b328 100644
--- a/opm/grid/cpgrid/CpGridData.cpp
+++ b/opm/grid/cpgrid/CpGridData.cpp
@@ -560,7 +560,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)
@@ -1602,8 +1602,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);
@@ -2090,9 +2090,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) {
@@ -2105,7 +2105,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]);
@@ -2331,7 +2331,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);
     }
 }
@@ -2724,7 +2724,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 84ed382ab..af2a016b0 100644
--- a/opm/grid/cpgrid/Entity.hpp
+++ b/opm/grid/cpgrid/Entity.hpp
@@ -163,7 +163,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;
@@ -382,7 +382,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];
 }
@@ -523,35 +523,23 @@ Dune::cpgrid::Geometry<3,3> Dune::cpgrid::Entity<codim>::geometryInFather() cons
         OPM_THROW(std::logic_error, "Entity has no father.");
     }
 
-    // Indices of corners in entity's geometry in father reference element.
-    static constexpr std::array<int,8> in_father_reference_elem_corner_indices = {0,1,2,3,4,5,6,7};
-    // 'static': The returned object Geometry<3,3> stores a pointer to in_father_reference_elem_corner_indices. Therefore,
-    // this variable is declared static to prolongate its lifetime beyond this function (static storage duration).
-
     auto idx_in_parent_cell = pgrid_ -> cell_to_idxInParentCell_[this->index()];
     if (idx_in_parent_cell !=-1) {
         const auto& cells_per_dim =  (*(pgrid_ -> level_data_ptr_))[this->level()] -> cells_per_dim_;
-        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;
-        // 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,
-                                                        corners_in_father_reference_elem_temp + 8);
+        std::array<FieldVector<double, 3>, 8> in_father_reference_elem_corners = pgrid_->getReferenceRefinedCorners(idx_in_parent_cell, cells_per_dim);
         // Compute the center of the 'local-entity'.
         Dune::FieldVector<double, 3> center_in_father_reference_elem = {0., 0.,0.};
         for (int corn = 0; corn < 8; ++corn) {
             for (int c = 0; c < 3; ++c)
             {
-                center_in_father_reference_elem[c] += corners_in_father_reference_elem_temp[corn][c]/8.;
+                center_in_father_reference_elem[c] += in_father_reference_elem_corners[corn][c]/8.;
             }
         }
         // Compute the volume of the 'local-entity'.
         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());
+                                           in_father_reference_elem_corners);
     }
     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 ed21cbca1..ea25711c1 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,7 @@ namespace Dune
 
             /// @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
@@ -413,17 +416,33 @@ 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)
+                : pos_(pos), vol_(vol)
             {
-                assert(allcorners_ && corner_indices);
+                assert(allcorners_view && corner_indices);
+                for(std::size_t i = 0; i < 8; ++i)
+                    allcorners_[i] = allcorners_view->data()[corner_indices[i]];
+            }
+
+            /// @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 corner positions of the geometry
+            Geometry(const GlobalCoordinate& pos,
+                     ctype vol,
+                     const std::array<GlobalCoordinate, 8>& allcorners)
+                : pos_(pos), vol_(vol)
+            {
+                for (std::size_t i = 0; i != 8; ++i)
+                    allcorners_[i] = cpgrid::Geometry<0, 3>{allcorners[i]};
             }
 
             /// 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_()
             {
             }
 
@@ -511,8 +530,7 @@ 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 allcorners_[cor].center();
             }
 
             /// Cell volume.
@@ -526,7 +544,7 @@ namespace Dune
             }
 
             /// Returns the centroid of the geometry.
-            const GlobalCoordinate& center() const
+            GlobalCoordinate center() const
             {
                 return pos_;
             }
@@ -537,7 +555,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 +591,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 +1040,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 +1061,7 @@ 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
+            std::array<cpgrid::Geometry<0, 3>, 8> allcorners_;
 
             /// @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 6f4f49559..a13fedbcf 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 b910a06b8..fc33f7662 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/processEclipseFormat.cpp b/opm/grid/cpgrid/processEclipseFormat.cpp
index 4f7c51db0..a8e637308 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 9b4f4b9c3..0f66bdf4d 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]);
+            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 74645353e..d42d33cfc 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});