Skip to content
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

[GeoMechanicsApplication] Extract list for constitutive matrices and stress vectors #12380

Merged
merged 15 commits into from
May 16, 2024
Merged

[GeoMechanicsApplication] Extract list for constitutive matrices and stress vectors #12380

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
6536da2
Extracted a function to calculate constitutive response
rfaasse May 14, 2024
2cb0d31
Formatting
rfaasse May 14, 2024
69d643e
Simplified constitutive matrix function
rfaasse May 14, 2024
687004c
Simplified constitutive matrix function more by calculating determina…
rfaasse May 14, 2024
7fa8f04
Changed signature such that it's clear what is changed by the function
rfaasse May 14, 2024
eb20b9f
Created new function for list of constitutive matrices in upw small s…
rfaasse May 15, 2024
2295b5d
Used new function in derived classes as well
rfaasse May 15, 2024
76acc15
Removed redundant determinant calculations
rfaasse May 15, 2024
f7315bd
Change signature of function to make it more symmetric
rfaasse May 15, 2024
1ffdc87
Merge remote-tracking branch 'origin/master' into geo/12319-extract-l…
rfaasse May 15, 2024
f7bf495
Remove unused variable and some redundant comments
rfaasse May 16, 2024
0c3ca03
Changed implementation of CalculateAnyOfMaterialResponse function to …
rfaasse May 16, 2024
4a68a07
Make sure the inputs are filled with zero vectors/matrices of the rig…
rfaasse May 16, 2024
155f48b
Merge remote-tracking branch 'origin/master' into geo/12319-extract-l…
rfaasse May 16, 2024
8007aaa
Processing review comments
rfaasse May 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
71 changes: 25 additions & 46 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -139,26 +139,20 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::InitializeNonLinearIteration(con

const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = this->CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients);
const auto strain_vectors = this->CalculateStrains(
auto strain_vectors = this->CalculateStrains(
deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain);
std::vector<Matrix> constitutive_matrices;
this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

At very first sight this function name is not very descriptive for me. I have to learn what it does reading the function itself.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree, I was struggling with the name. Would CalculateAnyOfStressesStrainsOrConstitutiveMatrices be better?

Variables.NContainer, Variables.DN_DXContainer,
strain_vectors, mStressVector, constitutive_matrices);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];
Variables.B = b_matrices[GPoint];
Variables.F = deformation_gradients[GPoint];
Variables.ConstitutiveMatrix = constitutive_matrices[GPoint];

// Compute infinitessimal strain
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

// set gauss points variables to constitutivelaw parameters
this->SetConstitutiveParameters(Variables, ConstitutiveParameters);

ConstitutiveParameters.SetStressVector(mStressVector[GPoint]);
mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters);
this->SaveGPConstitutiveTensor(ConstitutiveTensorContainer, Variables.ConstitutiveMatrix, GPoint);

// Compute DtStress
Expand Down Expand Up @@ -200,27 +194,19 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::FinalizeNonLinearIteration(const

const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = this->CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients);
const auto strain_vectors = this->CalculateStrains(
auto strain_vectors = this->CalculateStrains(
deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain);
std::vector<Matrix> constitutive_matrices;
this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters,
Variables.NContainer, Variables.DN_DXContainer,
strain_vectors, mStressVector, constitutive_matrices);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];

// Compute infinitessimal strain
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

// set gauss points variables to constitutivelaw parameters
this->SetConstitutiveParameters(Variables, ConstitutiveParameters);

// Compute ConstitutiveTensor
ConstitutiveParameters.SetStressVector(mStressVector[GPoint]);
mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters);
Variables.B = b_matrices[GPoint];
Variables.F = deformation_gradients[GPoint];
Variables.ConstitutiveMatrix = constitutive_matrices[GPoint];

// Compute DtStress
noalias(Variables.StrainVector) = prod(Variables.B, Variables.VelocityVector);
Expand Down Expand Up @@ -463,24 +449,21 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
const auto integration_coefficients =
this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer);
const auto deformation_gradients = this->CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
GeoMechanicsMathUtilities::CalculateDeterminants(deformation_gradients);
const auto strain_vectors = this->CalculateStrains(
auto strain_vectors = this->CalculateStrains(
deformation_gradients, b_matrices, Variables.DisplacementVector, Variables.UseHenckyStrain);
std::vector<Matrix> constitutive_matrices;
this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters,
Variables.NContainer, Variables.DN_DXContainer,
strain_vectors, mStressVector, constitutive_matrices);

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
// Compute Np, GradNpT, B and StrainVector
this->CalculateKinematics(Variables, GPoint);
Variables.B = b_matrices[GPoint];

// Compute infinitessimal strain
Variables.F = deformation_gradients[GPoint];
Variables.detF = determinants_of_deformation_gradients[GPoint];
Variables.StrainVector = strain_vectors[GPoint];

// set gauss points variables to constitutivelaw parameters
this->SetConstitutiveParameters(Variables, ConstitutiveParameters);
Variables.B = b_matrices[GPoint];
Variables.F = deformation_gradients[GPoint];
avdg81 marked this conversation as resolved.
Show resolved Hide resolved
Variables.StrainVector = strain_vectors[GPoint];
Variables.ConstitutiveMatrix = constitutive_matrices[GPoint];

GeoElementUtilities::CalculateNuMatrix<TDim, TNumNodes>(Variables.Nu, Variables.NContainer, GPoint);
GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
Comment on lines 468 to 469
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We checked that these functions don't change anything that the constitutive law uses (i.e. only Nu, while the constitutive law uses Np). Since we changed the order here, I'd like to double check with you: do you think this is correct @WPK4FEM?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

When there is a difference between Nu and Np ( e.g. for reason of different interpolation orders of u and Pw ) the constitutive law ( relation between strain and stress ) should be using Nu. For same order interpolation Nu and Np will be equal, still the correct one should be used.

Expand All @@ -489,10 +472,6 @@ void UPwSmallStrainFICElement<TDim, TNumNodes>::CalculateAll(MatrixType& rLeftHa
// Compute ShapeFunctionsSecondOrderGradients
this->CalculateShapeFunctionsSecondOrderGradients(FICVariables, Variables);

// Compute constitutive tensor and stresses
ConstitutiveParameters.SetStressVector(mStressVector[GPoint]);
mConstitutiveLawVector[GPoint]->CalculateMaterialResponseCauchy(ConstitutiveParameters);

this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint);

// set shear modulus from stiffness matrix
Expand Down
Loading
Loading