diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
index e587d898a82b..634838ad5b92 100644
--- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
+++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
@@ -769,16 +769,15 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const
         ElementVariables Variables;
         this->InitializeElementVariables(Variables, rCurrentProcessInfo);
 
+        const auto b_matrices = CalculateBMatrices(Variables.NContainer, Variables.DN_DXContainer);
+
         for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) {
             noalias(Variables.Np) = row(Variables.NContainer, GPoint);
-
+            Variables.B           = b_matrices[GPoint];
             Matrix J0, InvJ0;
             this->CalculateDerivativesOnInitialConfiguration(
                 Variables.detJInitialConfiguration, J0, InvJ0, Variables.GradNpTInitialConfiguration, GPoint);
 
-            // Calculating operator B
-            this->CalculateBMatrix(Variables.B, Variables.GradNpTInitialConfiguration, Variables.Np);
-
             // Compute infinitesimal strain
             this->CalculateCauchyStrain(Variables);
 
@@ -1236,21 +1235,14 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeElementVariables(ElementV
     KRATOS_CATCH("")
 }
 
-template <unsigned int TDim, unsigned int TNumNodes>
-void UPwSmallStrainElement<TDim, TNumNodes>::CalculateBMatrix(Matrix& rB, const Matrix& GradNpT, const Vector& Np) const
-{
-    rB = this->GetStressStatePolicy().CalculateBMatrix(GradNpT, Np, this->GetGeometry());
-}
-
 template <unsigned int TDim, unsigned int TNumNodes>
 std::vector<Matrix> UPwSmallStrainElement<TDim, TNumNodes>::CalculateBMatrices(
     const Matrix& NContainer, const GeometryType::ShapeFunctionsGradientsType& DN_DXContainer) const
 {
     std::vector<Matrix> result;
     for (unsigned int GPoint = 0; GPoint < DN_DXContainer.size(); ++GPoint) {
-        Matrix b_matrix;
-        this->CalculateBMatrix(b_matrix, DN_DXContainer[GPoint], row(NContainer, GPoint));
-        result.push_back(b_matrix);
+        result.push_back(this->GetStressStatePolicy().CalculateBMatrix(
+            DN_DXContainer[GPoint], row(NContainer, GPoint), this->GetGeometry()));
     }
 
     return result;
diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
index 1f33304574ec..4746593a537e 100644
--- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
+++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp
@@ -230,7 +230,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa
 
     void CalculatePermeabilityUpdateFactor(ElementVariables& rVariables);
 
-    void                CalculateBMatrix(Matrix& rB, const Matrix& GradNpT, const Vector& Np) const;
     std::vector<Matrix> CalculateBMatrices(const Matrix& NContainer,
                                            const GeometryType::ShapeFunctionsGradientsType& DN_DXContainer) const;
 
@@ -242,7 +241,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa
 
     virtual void CalculateAndAddCompressibilityMatrix(MatrixType& rLeftHandSideMatrix, ElementVariables& rVariables);
 
-    virtual void CalculateAndAddPermeabilityMatrix(MatrixType& rLeftHandSideMatrix, const ElementVariables& rVariables);
+    virtual void CalculateAndAddPermeabilityMatrix(MatrixType&             rLeftHandSideMatrix,
+                                                   const ElementVariables& rVariables);
 
     virtual void CalculateAndAddRHS(VectorType& rRightHandSideVector, ElementVariables& rVariables, unsigned int GPoint);