[GeoMechanicsApplication] Extract B matrix calculations #12337

rfaasse · 2024-05-02T09:49:27Z

📝 Description
In this PR we extract the calculation of the b-matrices for all integration points. This is done by:

Creating a CalculateBMatrices function, which calculates B for all integration points at once.
The calculation of the B matrix is removed from the CalculateKinematics function to avoid duplicate work.
Everywhere CalculateKinematics was called, a call to CalculateBMatrices is done before the loop, and in the loop Variables.B is set with the corresponding matrix in the list.
In a few places (see the following commit) the calculation of B was redundant, i.e. for the calculation of the outputs DEGREE_OF_SATURATION, EFFECTIVE_SATURATION, BISHOP_COEFFICIENT, DERIVATIVE_OF_SATURATION, RELATIVE_PERMEABILITY (since they are only related to pressure), GREEN_LAGRANGE_STRAIN_TENSOR/GREEN_LAGRANGE_STRAIN_VECTOR (because these use the deformation gradient instead of B) and the calculation of the mass matrix.

… made `CalculateBMatrix` const and used this construction in the `CalculateAll` of the UPwSmallStrain element

… we now calculate the b matrices up front and set the variables using that list

… functions

…rix-calculations

…'s only used in CalculateBMatrices

…since it's only used in CalculateBMatrices" This reverts commit 6d983ea.

…e strains, the mass matrix and water pressure based outputs

markelov208

Hi Richard, a very dedicated PR that is very easy to review. Just a small comment on possible improvement.

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

avdg81

This is a very clear and clean PR. In essence, you apply the same transformation in several locations. I have checked the usage of ElementVariables::B and as far as I could see you have modified the correct usages (i.e. where this member wasn't used in the first place, no additional changes are required).

avdg81 · 2024-05-02T11:46:31Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

 {
    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


Just to follow the Kratos Style Guide:

Suggested change

const Matrix& NContainer, const GeometryType::ShapeFunctionsGradientsType& DN_DXContainer) const

const Matrix& rNContainer, const GeometryType::ShapeFunctionsGradientsType& rDN_DXContainer) const

avdg81 · 2024-05-02T11:50:04Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

+        Matrix b_matrix;
+        this->CalculateBMatrix(b_matrix, DN_DXContainer[GPoint], row(NContainer, GPoint));
+        result.push_back(b_matrix);


I would suggest to let CalculateBMatrix return a Matrix object (rather than void). In that case, you can simplify the code as follows:

Suggested change

Matrix b_matrix;

this->CalculateBMatrix(b_matrix, DN_DXContainer[GPoint], row(NContainer, GPoint));

result.push_back(b_matrix);

result.push_back(this->CalculateBMatrix(DN_DXContainer[GPoint], row(NContainer, GPoint)));

Done (every reviewer suggested the same, such harmony 😄 )

avdg81 · 2024-05-02T11:58:05Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp

+    void                CalculateBMatrix(Matrix& rB, const Matrix& GradNpT, const Vector& Np) const;
+    std::vector<Matrix> CalculateBMatrices(const Matrix& NContainer,
+                                           const GeometryType::ShapeFunctionsGradientsType& DN_DXContainer) const;


The ordering of the function parameters for the shape functions (Np and NContainer) and the gradients of the shape functions (GradNpT and DN_DXContainer) is not consistent between CalculateBMatrix and CalculateBMatrices. Perhaps we should stick to the ordering adopted by CalculateBMatrix, since that order is consistent with the one adopted by StressStatePolicy::CalculateBMatrix.

Good point, done!

avdg81 · 2024-05-02T12:06:40Z

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

+        Matrix b_matrix;
+        this->CalculateBMatrix(b_matrix, DN_DXContainer[GPoint], row(NContainer, GPoint));
+        result.push_back(b_matrix);


Just to point out that here we may consider similar changes as proposed before:

Suggested change

Matrix b_matrix;

this->CalculateBMatrix(b_matrix, DN_DXContainer[GPoint], row(NContainer, GPoint));

result.push_back(b_matrix);

result.push_back(this->CalculateBMatrix(DN_DXContainer[GPoint], row(NContainer, GPoint)));

avdg81 · 2024-05-02T12:09:46Z

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp

@@ -287,7 +287,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub
    double CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const;
    double CalculateBiotCoefficient(const ElementVariables& rVariables, const bool& hasBiotCoefficient) const;

-    virtual void CalculateBMatrix(Matrix& rB, const Matrix& DNu_DX, const Vector& Np);
+    virtual void        CalculateBMatrix(Matrix& rB, const Matrix& DNu_DX, const Vector& Np) const;


Here, the virtual keyword is redundant, since this member is never overridden.

Ah, I only removed it in upw small strain, forgot in diff order, good catch!

WPK4FEM

Hi Richard,
Tried to go trough the CalculateKinematics places where B is not used afterwards. Could not find mistakes there. See my other remarks. Keeping similar function setup with return values and arguments in the same order would help me read faster.
For the b_matrices you correctly adopt snake_case, but the argument names for & lack the r's and ofter p is used where I see no reason for the p. Obviouly these are existing anomalies that only happen to coincide with the lines you had to adapt.
Regards, Wijtze Pieter

WPK4FEM · 2024-05-02T14:23:17Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp

+    void                CalculateBMatrix(Matrix& rB, const Matrix& GradNpT, const Vector& Np) const;
+    std::vector<Matrix> CalculateBMatrices(const Matrix& NContainer,
+                                           const GeometryType::ShapeFunctionsGradientsType& DN_DXContainer) const;


Probably only confusion for me:
the stress state policies have a CalculateBMatrix which return a matrix and take the 3 arguments rGradNpT, rNp and rGeometry
the elements have a CalculateBMatrix with return type void and 3 arguments rB, GradNpT and Np
and a CalculateBMatrices which returns a std::vector<Matrix> and 2 arguments NContainer and DN_DXContainer

The naming and order of the arguments is confusing for me. Please have the same order for GradNpT and Np
GradNpT and the elements of DN_DXContainer are the same thing.
The p I don't understand, we look at the N for displacement here, so Nu or just N

WPK4FEM · 2024-05-02T14:33:57Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

+        Matrix b_matrix;
+        this->CalculateBMatrix(b_matrix, DN_DXContainer[GPoint], row(NContainer, GPoint));
+        result.push_back(b_matrix);


With b_matrix as return value this becomes more readable for me:
result.push_back(this->CalculateBMatrix(DN_DXContainer[GPoint], row(NContainer, GPoint));

…nstead of using input args

…`CalculateBMatrix` function

…form kratos style guide

…-matrix-calculations

… PR)

…ck into this branch

WPK4FEM

The remarks about redundant or now incorrect comments can be made far more often than I did here. I think that changing that can be taken up in next PR's, because we are continuing to unravel these functionalities.

WPK4FEM · 2024-05-06T09:08:14Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp

@@ -91,6 +92,7 @@ void UPwUpdatedLagrangianElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef
    for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
        // Compute element kinematics B, F, GradNpT ...


Now this really has become a "rotting" comment.

WPK4FEM · 2024-05-06T09:09:34Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

@@ -149,10 +149,13 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeSolutionStep(const Proces
    // Create general parameters of retention law
    RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo);

+    const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
+
    // Loop over integration points
    for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
        // Compute Np, GradNpT, B and StrainVector


Maybe remove this comment, B is not computed here anymore.

WPK4FEM · 2024-05-06T09:10:27Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

@@ -253,9 +256,12 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeNonLinearIteration(const
    ElementVariables Variables;
    this->InitializeElementVariables(Variables, rCurrentProcessInfo);

+    const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
+
    // Loop over integration points


what else? Needless repeated comment on integration point loops.

WPK4FEM · 2024-05-06T09:13:41Z

applications/GeoMechanicsApplication/custom_elements/axisymmetric_stress_state.cpp

@@ -17,11 +17,11 @@
 namespace Kratos
 {

-Matrix AxisymmetricStressState::CalculateBMatrix(const Matrix&         rGradNpT,
-                                                 const Vector&         rNp,
+Matrix AxisymmetricStressState::CalculateBMatrix(const Matrix&         rDN_DX,


This is better, I did not understand the p in the old name GradNpT.
An alternative would have been GradNT ( I guess the T is for transposed ) DX also refers to a spatial gradient.

WPK4FEM · 2024-05-06T09:20:02Z

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

@@ -1666,9 +1683,6 @@ void SmallStrainUPwDiffOrderElement::CalculateKinematics(ElementVariables& rVari
    noalias(rVariables.DNu_DX) = rVariables.DNu_DXContainer[GPoint];


This now seems redundant. It was input for the B computation, which now already has taken place.

WPK4FEM · 2024-05-06T09:20:39Z

applications/GeoMechanicsApplication/custom_elements/steady_state_Pw_element.cpp

@@ -157,7 +157,7 @@ void SteadyStatePwElement<TDim, TNumNodes>::CalculateAll(MatrixType&        rLef

    const auto integration_coefficients =
        this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);
-
+    


Unintended addition of spaces?

rfaasse added 10 commits May 1, 2024 15:48

Added function for calculating B matrices for all integration points,…

1094267

… made `CalculateBMatrix` const and used this construction in the `CalculateAll` of the UPwSmallStrain element

Everywhere CalculateKinematics is used for the small strain elements,…

9288374

… we now calculate the b matrices up front and set the variables using that list

Added function for calculation of all bmatrices to diff order element

96d1ac5

Same exercise for the diff order elements

5f255ea

Removing the calculation of the bmatrix from the calculate kinematics…

f2e4b94

… functions

Merge remote-tracking branch 'origin/master' into 12319-extract-b-mat…

2adec13

…rix-calculations

Added this-> for classes who call a protected function of the base class

6d7210c

Removed CalculateBMatrix function, for UPwSmallStrainElement since it…

6d983ea

…'s only used in CalculateBMatrices

Revert "Removed CalculateBMatrix function, for UPwSmallStrainElement …

c8e7b42

…since it's only used in CalculateBMatrices" This reverts commit 6d983ea.

Removed redundant b matrix calculations for calculating Green Lagrang…

0550c35

…e strains, the mass matrix and water pressure based outputs

rfaasse self-assigned this May 2, 2024

rfaasse added GeoMechanics Issues related to the GeoMechanicsApplication Refactor When code is moved or rewrote keeping the same behavior labels May 2, 2024

rfaasse requested review from WPK4FEM, markelov208 and avdg81 May 2, 2024 09:57

rfaasse marked this pull request as ready for review May 2, 2024 10:01

markelov208 reviewed May 2, 2024

View reviewed changes

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp Outdated Show resolved Hide resolved

avdg81 reviewed May 2, 2024

View reviewed changes

WPK4FEM reviewed May 2, 2024

View reviewed changes

rfaasse changed the title ~~Geo/12319 extract b matrix calculations~~ [GeoMechanicsApplication] Extract B matrix calculations May 3, 2024

rfaasse added 6 commits May 3, 2024 15:25

Changed signature of CalculateBMatrix function to return the matrix i…

11e792d

…nstead of using input args

Switched arguments of CalculateBMatrices to be consistent with the …

3d02020

…`CalculateBMatrix` function

Made naming between elements/stress state policies consistent and con…

ff78ba9

…form kratos style guide

Merge remote-tracking branch 'origin/master' into geo/12319-extract-b…

ad007fd

…-matrix-calculations

Removed an unnecessary size check (non-blocking comment from previous…

496ced5

… PR)

Removed redundant statements which were added after merging master ba…

252bb0c

…ck into this branch

rfaasse requested review from avdg81, markelov208 and WPK4FEM May 6, 2024 08:59

rfaasse enabled auto-merge (squash) May 6, 2024 09:05

WPK4FEM approved these changes May 6, 2024

View reviewed changes

rfaasse merged commit e0b29dc into master May 6, 2024
11 checks passed

rfaasse deleted the geo/12319-extract-b-matrix-calculations branch May 6, 2024 09:25

rfaasse linked an issue May 6, 2024 that may be closed by this pull request

[GeoMechanicsApplication] Split the calculation of the inputs for the utility functions for element matrices and the calling of the functions #12319

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[GeoMechanicsApplication] Extract B matrix calculations #12337

[GeoMechanicsApplication] Extract B matrix calculations #12337

rfaasse commented May 2, 2024 •

edited

Loading

markelov208 left a comment

avdg81 left a comment

avdg81 May 2, 2024

rfaasse May 3, 2024

avdg81 May 2, 2024

rfaasse May 3, 2024

avdg81 May 2, 2024

rfaasse May 3, 2024

avdg81 May 2, 2024

rfaasse May 3, 2024

avdg81 May 2, 2024

rfaasse May 3, 2024

WPK4FEM left a comment

WPK4FEM May 2, 2024

rfaasse May 3, 2024

WPK4FEM May 2, 2024

rfaasse May 3, 2024

WPK4FEM left a comment

WPK4FEM May 6, 2024

WPK4FEM May 6, 2024

WPK4FEM May 6, 2024

WPK4FEM May 6, 2024

WPK4FEM May 6, 2024

WPK4FEM May 6, 2024

	const Matrix& NContainer, const GeometryType::ShapeFunctionsGradientsType& DN_DXContainer) const
	const Matrix& rNContainer, const GeometryType::ShapeFunctionsGradientsType& rDN_DXContainer) const

		@@ -91,6 +92,7 @@ void UPwUpdatedLagrangianElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef
		for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) {
		// Compute element kinematics B, F, GradNpT ...

		@@ -1666,9 +1683,6 @@ void SmallStrainUPwDiffOrderElement::CalculateKinematics(ElementVariables& rVari
		noalias(rVariables.DNu_DX) = rVariables.DNu_DXContainer[GPoint];

		@@ -157,7 +157,7 @@ void SteadyStatePwElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLef

		const auto integration_coefficients =
		this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);

[GeoMechanicsApplication] Extract B matrix calculations #12337

[GeoMechanicsApplication] Extract B matrix calculations #12337

Conversation

rfaasse commented May 2, 2024 • edited Loading

markelov208 left a comment

Choose a reason for hiding this comment

avdg81 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

WPK4FEM left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

WPK4FEM left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

rfaasse commented May 2, 2024 •

edited

Loading