[GeoMechanicsApplication] Zero mass for geo interface (#12640)

* Return zero mass matrix for UPwSmallStrainInterface. * Removed obsolete AssembleDensityMatrix.
KratosMultiphysics · Aug 22, 2024 · 246ca6b · 246ca6b
1 parent b9c905a
commit 246ca6b
Show file tree

Hide file tree

Showing 3 changed files with 4 additions and 75 deletions.
diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
@@ -151,69 +151,8 @@ template <unsigned int TDim, unsigned int TNumNodes>
 void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::CalculateMassMatrix(MatrixType& rMassMatrix,
                                                                           const ProcessInfo& rCurrentProcessInfo)
 {
-    KRATOS_TRY
-
-    const unsigned int N_DOF = this->GetNumberOfDOF();
-
-    // Resizing mass matrix
-    if (rMassMatrix.size1() != N_DOF) rMassMatrix.resize(N_DOF, N_DOF, false);
-    noalias(rMassMatrix) = ZeroMatrix(N_DOF, N_DOF);
-
-    const PropertiesType&                           r_prop = this->GetProperties();
-    const GeometryType&                             r_geom = this->GetGeometry();
-    const GeometryType::IntegrationPointsArrayType& IntegrationPoints =
-        r_geom.IntegrationPoints(mThisIntegrationMethod);
-    const unsigned int NumGPoints = IntegrationPoints.size();
-
-    // Defining shape functions and the determinant of the jacobian at all integration points
-    const Matrix& NContainer = r_geom.ShapeFunctionsValues(mThisIntegrationMethod);
-    Vector        detJContainer(NumGPoints);
-    r_geom.DeterminantOfJacobian(detJContainer, mThisIntegrationMethod);
-    const auto integration_coefficients =
-        this->CalculateIntegrationCoefficients(IntegrationPoints, detJContainer);
-
-    // Element variables
-    InterfaceElementVariables Variables;
-    this->InitializeElementVariables(Variables, r_geom, r_prop, rCurrentProcessInfo);
-
-    RetentionLaw::Parameters RetentionParameters(r_prop);
-
-    // Defining necessary variables
-    BoundedMatrix<double, TDim, TNumNodes * TDim> aux_density_matrix = ZeroMatrix(TDim, TNumNodes * TDim);
-    BoundedMatrix<double, TDim, TDim> density_matrix = ZeroMatrix(TDim, TDim);
-
-    array_1d<double, TDim> LocalRelDispVector;
-    array_1d<double, TDim> RelDispVector;
-    const double&          MinimumJointWidth = r_prop[MINIMUM_JOINT_WIDTH];
-
-    // Loop over integration points
-    for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
-        InterfaceElementUtilities::CalculateNuMatrix(Variables.Nu, NContainer, GPoint);
-
-        noalias(RelDispVector) = prod(Variables.Nu, Variables.DisplacementVector);
-
-        noalias(LocalRelDispVector) = prod(Variables.RotationMatrix, RelDispVector);
-
-        this->CalculateJointWidth(Variables.JointWidth, LocalRelDispVector[TDim - 1], MinimumJointWidth, GPoint);
-
-        Variables.IntegrationCoefficient = integration_coefficients[GPoint];
-
-        CalculateRetentionResponse(Variables, RetentionParameters, GPoint);
-
-        const auto density =
-            GeoTransportEquationUtilities::CalculateSoilDensity(Variables.DegreeOfSaturation, r_prop);
-
-        GeoElementUtilities::AssembleDensityMatrix(density_matrix, density);
-
-        noalias(aux_density_matrix) = prod(density_matrix, Variables.Nu);
-
-        // Adding contribution to Mass matrix
-        GeoElementUtilities::AssembleUUBlockMatrix(
-            rMassMatrix, prod(trans(Variables.Nu), aux_density_matrix) * Variables.JointWidth *
-                             Variables.IntegrationCoefficient);
-    }
-
-    KRATOS_CATCH("")
+    // Resizing mass matrix, no mass in these interface elements.
+    rMassMatrix = ZeroMatrix(this->GetNumberOfDOF(), this->GetNumberOfDOF());
 }
 
 template <unsigned int TDim, unsigned int TNumNodes>

diff --git a/applications/GeoMechanicsApplication/custom_utilities/README.md b/applications/GeoMechanicsApplication/custom_utilities/README.md
@@ -54,9 +54,9 @@ File transport_equation_utilities.hpp includes
 ### Mass Matrix (M)
 
 The mathematical definition is:
-$$M = \int_\Omega N_{u}^T \rho N_u d\Omega$$
+$$M = \int_\Omega \rho N_{u}^T N_u d\Omega$$
 
-Where $\Omega$ is the domain, $N_u$ is the displacement shape function and $\rho$ is the density matrix that holds density for all directions.
+Where $\Omega$ is the domain, $N_u$ is the displacement shape function and $\rho$ is the density.
 
 ### Stiffness Matrix (K)
 

diff --git a/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp
@@ -191,16 +191,6 @@ class GeoElementUtilities
         KRATOS_CATCH("")
     }
 
-    template <typename MatrixType>
-    static inline void AssembleDensityMatrix(MatrixType& DensityMatrix, double Density)
-    {
-        for (unsigned int idim = 0; idim < DensityMatrix.size1(); ++idim) {
-            for (unsigned int jdim = 0; jdim < DensityMatrix.size2(); ++jdim) {
-                DensityMatrix(idim, jdim) = Density;
-            }
-        }
-    }
-
     template <typename MatrixType1, typename MatrixType2>
     static inline void AssembleUUBlockMatrix(MatrixType1& rLeftHandSideMatrix, const MatrixType2& rUUBlockMatrix)
     {