diff --git a/dumux/geomechanics/el1p2c/el1p2clocalresidual.hh b/dumux/geomechanics/el1p2c/el1p2clocalresidual.hh
index 1cfddda6d776523b5f76b8df21f3320c16adf4b9..edab240b0cedfb8617108e27e69add402e441c1d 100644
--- a/dumux/geomechanics/el1p2c/el1p2clocalresidual.hh
+++ b/dumux/geomechanics/el1p2c/el1p2clocalresidual.hh
@@ -177,8 +177,9 @@ namespace Dumux
                 Scalar stabilizationTerm(0.0);
                 if(withStabilization_){
                 // calculate distance h between nodes i and j
-                DimVector hVec = this->element_().geometry().corner(fluxVars.face().j)
-                                      - this->element_().geometry().corner(fluxVars.face().i);
+                const auto geometry = this->element_().geometry();
+                DimVector hVec = geometry.corner(fluxVars.face().j)
+                               - geometry.corner(fluxVars.face().i);
                 Scalar h = hVec.two_norm();
                 stabilizationTerm = (h * h) /
                                     (4 * (fluxVars.lambda()
diff --git a/dumux/geomechanics/el2p/el2plocaloperator.hh b/dumux/geomechanics/el2p/el2plocaloperator.hh
index 1733b120eeb57d591f1901e77335d5988f179c30..8ca2e1ce8bb342ca5cf8170c1dd34d4c2f29d9b8 100644
--- a/dumux/geomechanics/el2p/el2plocaloperator.hh
+++ b/dumux/geomechanics/el2p/el2plocaloperator.hh
@@ -219,7 +219,8 @@ public:
                         Traits::LocalBasisType::Traits::RangeType RT_P;
 
         // select quadrature rule for the element geometry type and with the order=qorder
-        Dune::GeometryType geomType = eg.geometry().type();
+        const auto geometry = eg.geometry();
+        Dune::GeometryType geomType = geometry.type();
         const Dune::QuadratureRule<DF,dim>& rule = Dune::QuadratureRules<DF,dim>::rule(geomType,qorder);
 
         // loop over quadrature points
@@ -232,7 +233,7 @@ public:
 
 
              // get inverse transposed jacobian for quadrature point
-             const JacobianInverseTransposed jacobian = eg.geometry().jacobianInverseTransposed(it->position());
+             const JacobianInverseTransposed jacobian = geometry.jacobianInverseTransposed(it->position());
 
              // calculate shape function gradients at the quadrature point in global coordinates. This is done
              // by multiplying the reference element shape functions with the inverse transposed jacobian
@@ -297,7 +298,7 @@ public:
              RT_P pn = pw + MaterialLaw::pc(materialParams, sw);
              RT_P pEff;
 
-             const GlobalPosition& globalPos = eg.geometry().global(it->position());
+             const GlobalPosition& globalPos = geometry.global(it->position());
 
              // calculate change in effective pressure with respect to initial conditions pInit (pInit is negativ)
              pEff = pw*sw + pn*sn + model_.problem().pInit(globalPos, it->position(), eg.entity());
@@ -339,7 +340,7 @@ public:
              RF rhoDiff = volVars.density(nPhaseIdx) - volVars.density(wPhaseIdx);
 
              // geometric weight need for quadrature rule evaluation (numerical integration)
-             RF qWeight = it->weight() * eg.geometry().integrationElement(it->position());
+             RF qWeight = it->weight() * geometry.integrationElement(it->position());
 
              // evaluate basis functions
              std::vector<RT_V> vBasis(dispSize);
@@ -399,7 +400,7 @@ public:
                 // position of quadrature point in local coordinates of element
                 DimVector local = isIt->geometryInInside().global(it->position());
 
-                GlobalPosition globalPos = eg.geometry().global(local);
+                GlobalPosition globalPos = geometry.global(local);
 
                 // evaluate boundary condition type
                 BoundaryTypes boundaryTypes;
@@ -465,7 +466,7 @@ public:
 //                             this doesn't work: DimVector local = isIt->geometryInInside().global(face_refElement.position(j,codim-1));
                             DimVector local = refElement.template geometry<1>(fIdx).global(face_refElement.position(j, codim-1));
 
-                            GlobalPosition globalPos = eg.geometry().global(local);
+                            GlobalPosition globalPos = geometry.global(local);
 
                             // evaluate boundary condition type
                             BoundaryTypes boundaryTypes;
diff --git a/dumux/geomechanics/el2p/el2pmodel.hh b/dumux/geomechanics/el2p/el2pmodel.hh
index 5ab223d8ab0609274b1d98d8d9b0a8831416aba0..9bc4d6aa30087f5b1023d3b43382aafd3496bb37 100644
--- a/dumux/geomechanics/el2p/el2pmodel.hh
+++ b/dumux/geomechanics/el2p/el2pmodel.hh
@@ -429,8 +429,10 @@ public:
                 effPorosity[eIdx] +=elemVolVars[scvIdx].effPorosity / numScv;
             };
 
-            const GlobalPosition& cellCenter = eIt->geometry().center();
-            const GlobalPosition& cellCenterLocal = eIt->geometry().local(cellCenter);
+            const auto geometry = eIt->geometry();
+
+            const GlobalPosition& cellCenter = geometry.center();
+            const GlobalPosition& cellCenterLocal = geometry.local(cellCenter);
 
             deltaEffPressure[eIdx] = effectivePressure[eIdx] + this->problem().pInit(cellCenter, cellCenterLocal, *eIt);
             // determin changes in effective stress from current solution
@@ -439,7 +441,7 @@ public:
             displacementLFS.child(0).finiteElement().localBasis().evaluateJacobian(cellCenterLocal, vRefShapeGradient);
 
             // get jacobian to transform the gradient to physical element
-            const JacobianInverseTransposed jacInvT = eIt->geometry().jacobianInverseTransposed(cellCenterLocal);
+            const JacobianInverseTransposed jacInvT = geometry.jacobianInverseTransposed(cellCenterLocal);
             std::vector < Dune::FieldVector<RF, dim> > vShapeGradient(dispSize);
             for (size_t i = 0; i < dispSize; i++) {
                 vShapeGradient[i] = 0.0;
diff --git a/dumux/implicit/1p2c/1p2cfluxvariables.hh b/dumux/implicit/1p2c/1p2cfluxvariables.hh
index d2c53ce290fe9097b8512db89cdfdd6005a65719..e36a7cc1e6c7b40169eb96e7a7b391be084967c9 100644
--- a/dumux/implicit/1p2c/1p2cfluxvariables.hh
+++ b/dumux/implicit/1p2c/1p2cfluxvariables.hh
@@ -331,8 +331,9 @@ protected:
         }
         else {
             // use two-point gradients
-            const GlobalPosition &globalPosI = element.geometry().corner(face().i);
-            const GlobalPosition &globalPosJ = element.geometry().corner(face().j);
+            const auto geometry = element.geometry();
+            const GlobalPosition &globalPosI = geometry.corner(face().i);
+            const GlobalPosition &globalPosJ = geometry.corner(face().j);
             tmp = globalPosI;
             tmp -= globalPosJ;
             Scalar dist = tmp.two_norm();
diff --git a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
index fa38e3b148f0b0fbfa3e2570ab22a35ed54a7273..63b1e621d95cb1970f49fddc57df3e1911e76455 100644
--- a/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
+++ b/dumux/implicit/2pdfm/2pdfmfluxvariables.hh
@@ -201,11 +201,13 @@ private:
             int i = this->face().i;
             int j = this->face().j;
 
+            const auto geometry = element.geometry();
+
             // compute sum of pressure gradients for each phase
             for (int phaseIdx = 0; phaseIdx < numPhases; phaseIdx++)
             {
-                const GlobalPosition localIdx_i = element.geometry().corner(i);
-                const GlobalPosition localIdx_j = element.geometry().corner(j);
+                const GlobalPosition localIdx_i = geometry.corner(i);
+                const GlobalPosition localIdx_j = geometry.corner(j);
 
                 isFracture_ = problem.spatialParams().isEdgeFracture(element, fIdx);
 
diff --git a/dumux/implicit/2pdfm/2pdfmlocalresidual.hh b/dumux/implicit/2pdfm/2pdfmlocalresidual.hh
index 00306707b9a281a15a4c098b3017b8c7f6700bf0..3d670abdf9a5510381574ad3eedc422c977db5cd 100644
--- a/dumux/implicit/2pdfm/2pdfmlocalresidual.hh
+++ b/dumux/implicit/2pdfm/2pdfmlocalresidual.hh
@@ -130,7 +130,8 @@ public:
         /*
          * Calculate the fracture volume fraction wf = 0.5 * Fwidth * 0.5 * Length
          */
-        Dune::GeometryType geomType = elem.geometry().type();
+        const auto geometry = elem.geometry();
+        Dune::GeometryType geomType = geometry.type();
         const ReferenceElement &refElement = ReferenceElements::general(geomType);
 
         Scalar vol; //subcontrol volume
@@ -147,8 +148,8 @@ public:
             {
                 Scalar fracture_width = this->problem_().spatialParams().fractureWidth(elem, fIdx);
 
-                const GlobalPosition global_i = elem.geometry().corner(i);
-                const GlobalPosition global_j = elem.geometry().corner(j);
+                const GlobalPosition global_i = geometry.corner(i);
+                const GlobalPosition global_j = geometry.corner(j);
                 GlobalPosition diff_ij = global_j;
                 diff_ij -=global_i;
                 Scalar fracture_length = 0.5*diff_ij.two_norm();
@@ -165,7 +166,7 @@ public:
         storageFracture[nPhaseIdx]    = 0.0;
         storageMatrix[wPhaseIdx]    = 0.0;
         storageMatrix[nPhaseIdx]    = 0.0;
-        //        const GlobalPosition &globalPos = elem.geometry().corner(scvIdx);
+        //        const GlobalPosition &globalPos = geometry.corner(scvIdx);
 
         Scalar dsm_dsf = volVars.dsm_dsf();
         if (!this->problem_().useInterfaceCondition())
diff --git a/dumux/implicit/2pdfm/2pdfmvolumevariables.hh b/dumux/implicit/2pdfm/2pdfmvolumevariables.hh
index 908bc7daa5775299b8d8eae313a099305b831862..362887e34a5af2a9a935803c5cf1bf821fee5a69 100644
--- a/dumux/implicit/2pdfm/2pdfmvolumevariables.hh
+++ b/dumux/implicit/2pdfm/2pdfmvolumevariables.hh
@@ -258,7 +258,8 @@ public:
                                         int scvIdx)
     {
         Scalar volSCVFracture;
-        Dune::GeometryType geomType = element.geometry().type();
+        const auto geometry = element.geometry();
+        Dune::GeometryType geomType = geometry.type();
         const ReferenceElement &refElement = ReferenceElements::general(geomType);
 
         for (int fIdx=0; fIdx<refElement.size(1); fIdx++)
@@ -272,8 +273,8 @@ public:
             {
                 Scalar fracture_width = problem.spatialParams().fractureWidth();
 
-                const GlobalPosition global_i = element.geometry().corner(i);
-                const GlobalPosition global_j = element.geometry().corner(j);
+                const GlobalPosition global_i = geometry.corner(i);
+                const GlobalPosition global_j = geometry.corner(j);
                 GlobalPosition diff_ij = global_j;
                 diff_ij -=global_i;
                 //fracture length in the subcontrol volume is half d_ij
diff --git a/dumux/implicit/cellcentered/ccfvelementgeometry.hh b/dumux/implicit/cellcentered/ccfvelementgeometry.hh
index bc38552da58067fbf8ff27073afe589c00fe6821..41e521673d3d0929bf0fca95f9fdd82f62e4c2fc 100644
--- a/dumux/implicit/cellcentered/ccfvelementgeometry.hh
+++ b/dumux/implicit/cellcentered/ccfvelementgeometry.hh
@@ -100,7 +100,7 @@ public:
 
     void updateInner(const Element& element)
     {
-        const Geometry& geometry = element.geometry();
+        const Geometry geometry = element.geometry();
 
         elementVolume = geometry.volume();
         elementGlobal = geometry.center();
@@ -126,7 +126,7 @@ public:
     {
         updateInner(element);
 
-        const Geometry& geometry = element.geometry();
+        const Geometry geometry = element.geometry();
 
         bool onBoundary = false;
 
@@ -134,6 +134,8 @@ public:
         IntersectionIterator isEndIt = gridView.iend(element);
         for (IntersectionIterator isIt = gridView.ibegin(element); isIt != isEndIt; ++isIt)
         {
+            const auto isGeometry = isIt->geometry();
+
             // neighbor information and inner cvf data:
             if (isIt->neighbor())
             {
@@ -147,11 +149,12 @@ public:
                 scvFace.i = 0;
                 scvFace.j = scvfIdx + 1;
 
-                scvFace.ipGlobal = isIt->geometry().center();
+                scvFace.ipGlobal = isGeometry.center();
                 scvFace.ipLocal =  geometry.local(scvFace.ipGlobal);
                 scvFace.normal = isIt->centerUnitOuterNormal();
-                scvFace.normal *= isIt->geometry().volume();
-                scvFace.area = isIt->geometry().volume();
+                Scalar volume = isGeometry.volume();
+                scvFace.normal *= volume;
+                scvFace.area = volume;
 
                 GlobalPosition distVec = elementGlobal
                                        - neighbors[scvfIdx+1]->geometry().center();
@@ -179,11 +182,12 @@ public:
                 int bfIdx = isIt->indexInInside();
                 SubControlVolumeFace& bFace = boundaryFace[bfIdx];
 
-                bFace.ipGlobal = isIt->geometry().center();
+                bFace.ipGlobal = isGeometry.center();
                 bFace.ipLocal =  geometry.local(bFace.ipGlobal);
                 bFace.normal = isIt->centerUnitOuterNormal();
-                bFace.normal *= isIt->geometry().volume();
-                bFace.area = isIt->geometry().volume();
+                Scalar volume = isGeometry.volume();
+                bFace.normal *= volume;
+                bFace.area = volume;
                 bFace.i = 0;
                 bFace.j = 0;
 
diff --git a/dumux/implicit/common/implicitvelocityoutput.hh b/dumux/implicit/common/implicitvelocityoutput.hh
index df04d4277a16e4d1520b1b10597625aea3b3f74f..d76bdde2f1c45be9aebf4e135c97095cec68ecc1 100644
--- a/dumux/implicit/common/implicitvelocityoutput.hh
+++ b/dumux/implicit/common/implicitvelocityoutput.hh
@@ -148,7 +148,9 @@ public:
     {
         if (velocityOutput_)
         {
-            Dune::GeometryType geomType = element.geometry().type();
+            const auto geometry = element.geometry();
+
+            Dune::GeometryType geomType = geometry.type();
             const ReferenceElement &referenceElement
                 = ReferenceElements::general(geomType);
 
@@ -157,7 +159,7 @@ public:
 
             // get the transposed Jacobian of the element mapping
             const typename Element::Geometry::JacobianTransposed jacobianT2 =
-                element.geometry().jacobianTransposed(localPos);
+                geometry.jacobianTransposed(localPos);
 
             if (isBox)
             {
@@ -172,7 +174,7 @@ public:
 
                     // Transformation of the global normal vector to normal vector in the reference element
                     const typename Element::Geometry::JacobianTransposed jacobianT1 =
-                        element.geometry().jacobianTransposed(localPosIP);
+                        geometry.jacobianTransposed(localPosIP);
 
                     FluxVariables fluxVars(problem_,
                                            element,
@@ -213,7 +215,7 @@ public:
                     Dune::FieldVector<CoordScalar, dimWorld> scvVelocity(0);
 
                     jacobianT2.mtv(scvVelocities[scvIdx], scvVelocity);
-                    scvVelocity /= element.geometry().integrationElement(localPos)*cellNum_[vIdxGlobal];
+                    scvVelocity /= geometry.integrationElement(localPos)*cellNum_[vIdxGlobal];
                     // add up the wetting phase subcontrolvolume velocities for each vertex
                     velocity[vIdxGlobal] += scvVelocity;
                 }
@@ -319,7 +321,7 @@ public:
                 Dune::FieldVector<Scalar, dimWorld> scvVelocity(0);
                 jacobianT2.mtv(refVelocity, scvVelocity);
 
-                scvVelocity /= element.geometry().integrationElement(localPos);
+                scvVelocity /= geometry.integrationElement(localPos);
 
 #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 4)
                 int eIdxGlobal = problem_.elementMapper().index(element);