New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

[GeoMechanicsApplication] Extract a static utility function for the calculation of the Mass Matrix (M) #12299

Merged

markelov208 merged 66 commits into master from geo/12279-extract-function-mass-matrix

May 10, 2024

Contributor

markelov208 commented Apr 19, 2024 •

edited

Loading

📝 Description
An utility function to calculate the mass matrix has been extracted.

🆕 Changelog
Please summarize the changes in one list to generate the changelog:
E.g.

The calculation of the mass matrix is available as a utility function CalculateMassMatrix. Its input are geometry, solid densities and integration coefficients for all integration points.
The solid densitis are calculated with utility function CalculateSoilDensities and the coefficients with utility function CalculateIntegrationCoefficientsInitialConfiguration.
CalculateMassMatrix function returns a matrix for UU part. This matrix shall be assembled to the mass matrix.
The utility function is documented and unit-tested for UPwSmallStrainElement and SmallStrainUPwDiffOrderElement.


          implemented CalculateMassMatrix for small_strain_U_Pw_diff_order_element

50f3582

markelov208 self-assigned this

Gennady Markelov added 15 commits

April 19, 2024 21:34


          removed CalculateMassMatrix

ea4cc17


          replaced rGeom with this->GetGeometry()

0dd0c5a


          added description of CalculateMassMatrix

fc82a72


          added Np, changed from Camel to snake_case

83f27a1


          removed CalculateMassMatrix

122e145


          removed commented lines

d1ea2e1


          removed Np vector

27b872a


          changed in accordance SOnarCloud:

561390d

- raw pointer
- index++
- 4 for loops into 2 for loops


          extracted CalculateMassMatrix for U_Pw_small_strain_element

2802e2b


          replaced local CalculateSoilDensity functions with the static function

db42ede


          fixed extracted CalculateMassMatrix for U_Pw_small_strain_element

f632a3b


          fixed return of MassMatrix

1fec4e6


          removed unique_ptr from function parameters

4e56457


          Revert "removed CalculateMassMatrix"

684103f

This reverts commit 122e145.

# Conflicts:
#	applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp
#	applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp

Revert "removed CalculateMassMatrix"

This reverts commit ea4cc17.


          extracted CalculateDegreeOfSaturation and CalculateIntegrationCoeffic…

bf235ca

…ientInitialConfiguration

markelov208 requested review from rfaasse, avdg81 and WPK4FEM

April 26, 2024 08:13

rfaasse reviewed

View reviewed changes

Contributor

rfaasse left a comment

I really like this PR, it cleaned up nicely and took away unnecessary responsibilities from the element (and the function is already on the element level as opposed to integration points)!

I have suggestions for improvements (apart from the discussed unit tests that indeed still need to be added) of the calculation code now that it's moved to the utilities and I think we might need a bit more discussion on the diff order/non-diff order calculations, but this is a very nice step!

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp Outdated

-                  MatrixType MassMatrix(ElementSize, ElementSize);
-                  this->CalculateMassMatrix(MassMatrix, rCurrentProcessInfo);
+                  MatrixType MassMatrix = GeoTransportEquationUtilities::CalculateMassMatrixDiffOrder(

Contributor

rfaasse Apr 26, 2024

Since this is a local variable, this should be snake case (problem of course is that it seems this class very inconsistent with that right now, maybe we can fix it only in the functions we touch?)

Suggested change

      
                MatrixType MassMatrix = GeoTransportEquationUtilities::CalculateMassMatrixDiffOrder(
          
                MatrixType mass_matrix = GeoTransportEquationUtilities::CalculateMassMatrixDiffOrder(

Contributor Author

markelov208 Apr 26, 2024

Thank you. Yesterday I made the same comments to Mohamed ;(

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

Comment on lines 167 to 169

+                      for (SizeType i = 0; i < NumberPNodes; ++i) {
+                          solution_vector[i] = rGeom[i].FastGetSolutionStepValue(SolutionVariable);
+                      }

Contributor

rfaasse Apr 26, 2024

This can be a transform 😃 :

Suggested change

      
                    for (SizeType i = 0; i < NumberPNodes; ++i) {
          
                        solution_vector[i] = rGeom[i].FastGetSolutionStepValue(SolutionVariable);
          
                    }
          
                    std::transform(rGeom.begin(), rGeom.end(), solution_vector.begin(),
          
                                   [&SolutionVariable](const auto& node) {
          
                        return node.FastGetSolutionStepValue(SolutionVariable);
          
                    });

Contributor Author

markelov208 Apr 26, 2024 •

edited

Loading

Great! Used it.

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

+                      return inner_prod(rNp, rPressureVector);
+                  }
+                  static Vector GetSolutionVector(SizeType NumberPNodes, const Geometry<Node>& rGeom, const Variable<double>& SolutionVariable)

Contributor

rfaasse Apr 26, 2024

Nitpicking:

Suggested change

      
                static Vector GetSolutionVector(SizeType NumberPNodes, const Geometry<Node>& rGeom, const Variable<double>& SolutionVariable)
          
                static Vector GetSolutionVector(SizeType NumberPNodes, const Geometry<Node>& rGeom, const Variable<double>& rSolutionVariable)

Contributor Author

markelov208 Apr 26, 2024

Thank you, changed.

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

+                      return mass_matrix;
+                  }
+                  static double CalculateFluidPressure(const Vector& rNp, const Vector& rPressureVector)

Contributor

rfaasse Apr 26, 2024

I think this function can be private

Contributor Author

markelov208 Apr 26, 2024

Indeed, made more functions private. No need to make unit tests!

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

+                      return inner_prod(rNp, rPressureVector);
+                  }
+                  static Vector GetSolutionVector(SizeType NumberPNodes, const Geometry<Node>& rGeom, const Variable<double>& SolutionVariable)

Contributor

rfaasse Apr 26, 2024

I think this function can be private

Contributor Author

markelov208 Apr 26, 2024

done

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

@@ @@ -79,6 +86,174 @@ class GeoTransportEquationUtilities @@
                   {
                       return -PORE_PRESSURE_SIGN_FACTOR * BiotModulusInverse * outer_prod(rNp, rNp) * IntegrationCoefficient;
                   }
+                  static Matrix CalculateMassMatrixDiffOrder(const Geometry<Node>&          rGeom,

Contributor

rfaasse Apr 26, 2024

As discussed, it is probably a good idea to let one of the functions call the other, such that we avoid duplications which are still there at this moment. Ideally the utilities shouldn't care about diff order/non-diff order, would it be possible to prepare the inputs for this functions in such a way in the diff order element that we can use the standard version? That way, all real 'diff-order' functionality stays in that element.

However, I haven't spend as much time as you in this piece of code and I think this is already a nice improvement, so we could leave it to a later PR in my opinion if that would be too much work for this one.

What is your opinion here @avdg81?

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

Comment on lines 151 to 155

+                      for (SizeType i = 0; i < number_U_nodes * dimension; ++i) {
+                          for (SizeType j = 0; j < number_U_nodes * dimension; ++j) {
+                              mass_matrix(i, j) += mass_matrix_contribution(i, j);
+                          }
+                      }

Contributor

rfaasse Apr 26, 2024

This looks very similar to this function in element utilities, maybe we can re-use it?

Kratos/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp

Lines 600 to 611 in f8087c7

    
           template <typename MatrixType1, typename MatrixType2> 
        
           static void AddMatrixAtPosition(const MatrixType1& rSourceMatrix, 
        
                                           MatrixType2&       rDestinationMatrix, 
        
                                           std::size_t        RowOffset, 
        
                                           std::size_t        ColumnOffset) 
        
           { 
        
               for (auto i = std::size_t{0}; i < rSourceMatrix.size1(); ++i) { 
        
                   for (auto j = std::size_t{0}; j < rSourceMatrix.size2(); ++j) { 
        
                       rDestinationMatrix(i + RowOffset, j + ColumnOffset) += rSourceMatrix(i, j); 
        
                   } 
        
               } 
        
           }

Although I see it's private, what is your opinion @avdg81 (since you added that function)?

Contributor Author

markelov208 Apr 26, 2024 •

edited

Loading

Thank for the info. I used AssembleUUBlockMatrix than both functions CalculateMassMatrix became more similar.

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Show resolved Hide resolved

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

+                          rGeom.IntegrationPoints(IntegrationMethod);
+                      const unsigned int number_G_points = integration_points.size();
+                      Matrix N_container = rGeom.ShapeFunctionsValues(IntegrationMethod);

Contributor

rfaasse Apr 26, 2024

Since the ShapeFunctionsValues function returns a const Matrix&, I'd suggest to make this const& to avoid the copy

Suggested change

      
                    Matrix N_container = rGeom.ShapeFunctionsValues(IntegrationMethod);
          
                    const Matrix& N_container = rGeom.ShapeFunctionsValues(IntegrationMethod);

Contributor Author

markelov208 Apr 26, 2024

changed.

applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp Outdated

+                      const unsigned int number_G_points = integration_points.size();
+                      Matrix N_container = rGeom.ShapeFunctionsValues(IntegrationMethod);
+                      Vector pressure_vector =

Contributor

rfaasse Apr 26, 2024

This can be const

Suggested change

      
                    Vector pressure_vector =
          
                    const Vector pressure_vector =

Contributor Author

markelov208 Apr 26, 2024

done

Gennady Markelov added 9 commits

April 26, 2024 15:19


          changed MassMatrix to mass_matrix

d8abb63


          added CalculateNuMatrix for Matrxi1&Matrix2 and clang-formatting

1e5a696


          changed in accordance with the review comments

7d7ac07


          added test for CalculateSoilDensity function

cfd8da6


          removed index=i-Dim

a14951c


          replaced pointer with reference

3d5bd1f


          replaced pointer with reference -cont

6417b44


          removed redundant NumberPNodes

f59640e


          changed RetentionParameters to snake_case

4d877be

Gennady Markelov added 9 commits

May 3, 2024 23:50


          deactivated 1d_consolidation test

e53ba8c


          used pragma once

419bbfe


          removed parentheses

35d8b43


          used arithmetic operations instead the final numbers

1e64d18


          removed commented includes

56b5828


          Merge branch 'master' into geo/12279-extract-function-mass-matrix

de3af9a


          renamed some local variables and removed unused includes

500d098


          removed unused includes in equations_of_motion_utilities.hpp

63cc3dd


          added s in CalculateIntegrationCoefficientsInitialConfiguration

f22c8a4

markelov208 requested review from rfaasse and WPK4FEM

May 5, 2024 07:36

Gennady Markelov added 12 commits

May 6, 2024 20:48


          fixed porosity and soil density description

3c04353


          restored CalculateMassMatrix for DiffOrder element

5879c6b


          used available elements' CalculateIntegrationCoefficients and removed…

2e1ed28

… rGeom from CalculateMassMatrix input


          split equation_of_motion_utilities.hpp into h and cpp files

fc15f3e


          removed unused includes in test_mass_matrix

8ed6b14


          removed unused dN_De in equation_of_motion_utilities

7f10d5f


          moved getting of pressure solution vector from transport_equation_uti…

99b0fba

…lities to elements


          removed extra qualification for GetPressureSolutionVector

19559c6


          removed private functions in transport_equation_utilities.hpp

383176a


          added this-> for CalculateIntegrationCoefficients call

c5c44ba


          extracted CalculateDegreesSaturation from CalculateSoilDensities


          changed setting of properties to remove a code smell

9012aeb

Contributor Author

markelov208 commented May 7, 2024 •

edited

Loading

Hi guys, following our conversation this Monday:
1.Now CalculateIntegrationCoefficients function that are developed by Anne for each element are used. Input is provided with CalculateDetJsInitialConfiguration.
2. CalculateMassMatrix function is restored for SmallStrainUPwDiffOrder
3. rGeom is removed from equation_of_motion CalculateMassMatrix
4. equation_of_motion_utilities.hpp is split into h and cpp files
5. transport_equation CalculateSolidDensities now calculates only soil densities.
6. A new function CalculateDegreesSaturation calculates degrees of saturation. This function is still in transport_equation. There are a lot of similar places with RetentionLaw in elements. One of CalculateDegreesSaturation input is pressure_vector which is provided with GetPressureSolution() now. This function is done for each element like CalculateIntegrationCoefficients.
7. Readme is corrected.
8. test_mass_matrix calls element CalculateMassMatrix; therefore all functions (soil densities, degrees saturation etc. ) are tested together.
9. It will be useful to create an issue to extract a single function CalculateFluidPressure. There are a number of its copies across elements.

rfaasse reviewed

View reviewed changes

Contributor

rfaasse left a comment

Thank you for incorporating the changes, I think the function is turning out to be very clear and readable and you extracted and cleaned up a lot of other duplications in this PR! I only have a few comments left about includes that could be removed/moved from the header.

applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.cpp Outdated

Comment on lines 21 to 41

+              Matrix GeoEquationOfMotionUtilities::CalculateMassMatrix(SizeType      dimension,
+                                                                       SizeType      number_U_nodes,
+                                                                       SizeType      NumberIntegrationPoints,
+                                                                       const Matrix& Nu_container,
+                                                                       const Vector& rSolidDensities,
+                                                                       const std::vector<double>& rIntegrationCoefficients)
+              {
+                  const SizeType block_element_size = number_U_nodes * dimension;
+                  Matrix         Nu                 = ZeroMatrix(dimension, block_element_size);
+                  Matrix         aux_density_matrix = ZeroMatrix(dimension, block_element_size);
+                  Matrix         density_matrix     = ZeroMatrix(dimension, dimension);
+                  Matrix         mass_matrix        = ZeroMatrix(block_element_size, block_element_size);
+                  for (unsigned int g_point = 0; g_point < NumberIntegrationPoints; ++g_point) {
+                      GeoElementUtilities::AssembleDensityMatrix(density_matrix, rSolidDensities(g_point));
+                      GeoElementUtilities::CalculateNuMatrix(dimension, number_U_nodes, Nu, Nu_container, g_point);
+                      noalias(aux_density_matrix) = prod(density_matrix, Nu);
+                      mass_matrix += prod(trans(Nu), aux_density_matrix) * rIntegrationCoefficients[g_point];
+                  }
+                  return mass_matrix;
+              }

Contributor

rfaasse May 8, 2024

This became a very easy-to-read and generic function!

Contributor Author

markelov208 May 8, 2024

Thank you, it is a real team work!

applications/GeoMechanicsApplication/custom_utilities/equation_of_motion_utilities.h Outdated

+              // Project includes
+              // Application includes
+              #include "custom_utilities/element_utilities.hpp"

Contributor

rfaasse May 8, 2024

Minor, but I think this include can move to the cpp file

Contributor Author

markelov208 May 8, 2024

I moved it to the cpp file but I had to add other h files. Nevertheless, a number of dependencies is decreased.

applications/GeoMechanicsApplication/tests/cpp_tests/test_mass_matrix.cpp

+              namespace Kratos::Testing
+              {
+              KRATOS_TEST_CASE_IN_SUITE(CalculateMassMatrix2D6NDiffOrderGivesCorrectResults, KratosGeoMechanicsFastSuite)

Contributor

rfaasse May 8, 2024

I really like that we have these unit tests to make sure these functions keep doing what they're doing now!

applications/GeoMechanicsApplication/tests/cpp_tests/test_soil_density.cpp Outdated

Comment on lines 13 to 16

+              #include "containers/model.h"
+              #include "custom_utilities/transport_equation_utilities.hpp"
+              #include "includes/checks.h"
+              #include "testing/testing.h"

Contributor

rfaasse May 8, 2024

I'm not sure if we need the containers/model.h, could we try to remove it?

Suggested change

      
            #include "containers/model.h"
          
            #include "custom_utilities/transport_equation_utilities.hpp"
          
            #include "includes/checks.h"
          
            #include "testing/testing.h"
          
            #include "custom_utilities/transport_equation_utilities.hpp"
          
            #include "includes/checks.h"
          
            #include "testing/testing.h"

Contributor Author

markelov208 May 8, 2024

You are right, yes we can and I made it now.

applications/GeoMechanicsApplication/tests/cpp_tests/test_utilities/model_setup_utilities.h Outdated

Comment on lines 13 to 16

               #pragma once
               #include "geo_aliases.h"
+              #include "includes/properties.h"

Contributor

rfaasse May 8, 2024

This include we might be able to move to the cpp file.

Contributor Author

markelov208 May 8, 2024

You are right, done.


          changed include statements following the review comments

cf9f0d9

WPK4FEM approved these changes

View reviewed changes

Contributor

WPK4FEM left a comment

Hi Gennady,
I think a large improvement is made with this pull request. Please judge what you think of the remarks is still useful and what is more suited for next actions.
Wijtze Pieter

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

Comment on lines +462 to +464

+                  if (rMassMatrix.size1() != element_size || rMassMatrix.size2() != element_size)
+                      rMassMatrix.resize(element_size, element_size, false);
+                  noalias(rMassMatrix) = ZeroMatrix(element_size, element_size);

Contributor

WPK4FEM May 10, 2024

The resize will test for the size itself, so no need for check before.
ZeroMatrix will set the size correctly, but I'm not certain about the noalias.
I think line 462 can be omitted

The same happens for the dampingmatrix below at 494

Contributor Author

markelov208 May 10, 2024 •

edited

Loading

noalias does not change the size. I think it will be a decision on all noalias'es later.

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

Comment on lines +492 to 501

    
                  if (rDampingMatrix.size1() != element_size)

                      rDampingMatrix.resize(element_size, element_size, false);

                  noalias(rDampingMatrix) = ZeroMatrix(element_size, element_size);

                  if (rProp.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_ALPHA] * MassMatrix;

                  else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * MassMatrix;

                  if (r_prop.Has(RAYLEIGH_ALPHA)) noalias(rDampingMatrix) += r_prop[RAYLEIGH_ALPHA] * mass_matrix;

                  else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_ALPHA] * mass_matrix;

                  if (rProp.Has(RAYLEIGH_BETA)) noalias(rDampingMatrix) += rProp[RAYLEIGH_BETA] * StiffnessMatrix;

                  if (r_prop.Has(RAYLEIGH_BETA))

                      noalias(rDampingMatrix) += r_prop[RAYLEIGH_BETA] * StiffnessMatrix;

                  else noalias(rDampingMatrix) += rCurrentProcessInfo[RAYLEIGH_BETA] * StiffnessMatrix;

Contributor

WPK4FEM May 10, 2024

When first the Rayleigh parameters are determined, the setting of the damping matrix becomes like the formula. e.g.

const double rayleigh_alpha = r_prop.Has(RAYLEIGH_ALPHA) ? r_prop[RAYLEIGH_ALPHA] : rCurrentProcessInfo[RAYLEIGH_ALPHA];
const double rayleigh_beta = r_prop.Has(RAYLEIGH_BETA) ? r_prop[RAYLEIGH_BETA] : rCurrentProcessInfo[RAYLEIGH_BETA];

rDampingMatrix = rayleigh_alpha * mass_matrix + rayleigh_beta + stiffness_matrix;

( I didn't check my syntax here, just trying to write down the intention )

As focus is on the mass matrix first, this is probably for a next pull request.

Contributor Author

markelov208 May 10, 2024

Agreed that it is better in any sense: visibility, understanding etc.

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

Comment on lines +2101 to +2109

+              Vector SmallStrainUPwDiffOrderElement::GetPressureSolutionVector()
+              {
+                  Vector result(mpPressureGeometry->PointsNumber());
+                  std::transform(this->GetGeometry().begin(),
+                                 this->GetGeometry().begin() + mpPressureGeometry->PointsNumber(), result.begin(),
+                                 [](const auto& node) { return node.FastGetSolutionStepValue(WATER_PRESSURE); });
+                  return result;
+              }

Contributor

WPK4FEM May 10, 2024

This is a very useful functionality, variants of it exist in custom_ulitities/element_utilities.hpp, the GetNodalVariableVector and GetNodalVariableMatrix. I think that is the place where these belong, but your implementation is more elegant.

Contributor Author

markelov208 May 10, 2024

I can't place GetPressureSolutionVector in element_utilities because it is a specific for DIffOrder and UPw. Thanks go to @avdg81 and @rfaasse

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

Comment on lines +1831 to +1839

+              template <unsigned int TDim, unsigned int TNumNodes>
+              Vector UPwSmallStrainElement<TDim, TNumNodes>::GetPressureSolutionVector()
+              {
+                  Vector result(TNumNodes);
+                  std::transform(this->GetGeometry().begin(), this->GetGeometry().end(), result.begin(),
+                                 [](const auto& node) { return node.FastGetSolutionStepValue(WATER_PRESSURE); });
+                  return result;
+              }

Contributor

WPK4FEM May 10, 2024

Very similar (but nicer written) to the GetNodalVariableVector or Matrix functions in element_utilities.hpp

Maybe this functionality belongs there.

rfaasse approved these changes

View reviewed changes

Contributor

rfaasse left a comment

Thank you for all your hard work on this PR! From my side, this is ready to go!

markelov208 merged commit 9c4a7b0 into master

11 checks passed

markelov208 deleted the geo/12279-extract-function-mass-matrix branch

May 10, 2024 09:25

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet