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.

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

[SMApp] Adding a mixed Displacement/Strain element #12897

Open
wants to merge 155 commits into
base: master
Choose a base branch
from

Conversation

AlejandroCornejo
Copy link
Member

@AlejandroCornejo AlejandroCornejo commented Dec 4, 2024

📝 Description
Adding a mixed u/E high accuracy element in small strains. This element is based on:

  • Finite element modeling of quasi-brittle cracks in 2D and 3D with enhanced strain accuracy, M. Cervera, G. Barbat and M. Chiumenti, Computational Mechanics, (60) 767-796, 2017. DOI: https://doi.org/10.1007/s00466-017-1438-8

with some improvements such as the usageof tangent constitutive matrix in the balance of linear mom equation and the elastic constitutive matrix in the strain compatibility equation.

// Compute the symmetric gradient of the displacements
noalias(rKinVariables.SymmGradientDispl) = prod(rKinVariables.B, rKinVariables.NodalDisplacements);

noalias(rKinVariables.NodalStrain) = prod(rKinVariables.N_epsilon, rKinVariables.NodalStrains);
Copy link
Member

Choose a reason for hiding this comment

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

Are we sure of this? Maybe I'm wrong, but it seems that we are using the same container for the nodal values and the integration point interpolation.

Copy link
Member Author

Choose a reason for hiding this comment

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

noalias(rKinVariables.NodalStrain) = prod(rKinVariables.N_epsilon, rKinVariables.NodalStrains); The left term is a interpolated strain Vector and the right one is the shape function matrix and a element size strain vector

Copy link
Member

Choose a reason for hiding this comment

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

Ufff... now I see the s at the end... This is very misleading IMO.

Taking into account that NodalStrain is only used to calculate the equivalent strain, I'd avoid saving this in the kinematic variables container and directly do the interpolation it in the CalculateEquivalentStrain.

Comment on lines 74 to 141
protected:

/**
* Internal variables used in the kinematic calculations
*/
struct KinematicVariables
{
Vector N;
Matrix B;
Matrix N_epsilon; // Used to interpolate nodal strains
double detJ0;
Matrix J0;
Matrix InvJ0;
Matrix DN_DX;
Vector NodalDisplacements; // Displacement DoFs -> U
Vector NodalStrains; // Strains stored at the nodes (strain DoFs) -> E
Vector EquivalentStrain; // Stabilized strain field E = (1-tau) N_e · E + tau Bu · U
Vector SymmGradientDispl; // Symmetric gradient of the nodal displacements: Bu·U
Vector NodalStrain; // N_e · E

/**
* The default constructor
* @param StrainSize The size of the strain vector in Voigt notation
* @param Dimension The problem dimension: 2D or 3D
* @param NumberOfNodes The number of nodes of the element
*/
KinematicVariables(
const SizeType StrainSize,
const SizeType Dimension,
const SizeType NumberOfNodes
)
{
detJ0 = 1.0;
N = ZeroVector(NumberOfNodes);
B = ZeroMatrix(StrainSize, Dimension * NumberOfNodes);
DN_DX = ZeroMatrix(NumberOfNodes, Dimension);
J0 = ZeroMatrix(Dimension, Dimension);
InvJ0 = ZeroMatrix(Dimension, Dimension);
NodalDisplacements = ZeroVector(Dimension * NumberOfNodes);
NodalStrains = ZeroVector(NumberOfNodes * StrainSize);
EquivalentStrain = ZeroVector(StrainSize);
SymmGradientDispl = ZeroVector(StrainSize);
NodalStrain = ZeroVector(StrainSize);
N_epsilon = ZeroMatrix(StrainSize, StrainSize * NumberOfNodes);
}
};

/**
* Internal variables used in the kinematic calculations
*/
struct ConstitutiveVariables
{
Vector StrainVector;
Vector StressVector;
Matrix D;

/**
* The default constructor
* @param StrainSize The size of the strain vector in Voigt notation
*/
ConstitutiveVariables(const SizeType StrainSize)
{
StrainVector = ZeroVector(StrainSize);
StressVector = ZeroVector(StrainSize);
D = ZeroMatrix(StrainSize, StrainSize);
}
};

Copy link
Member

Choose a reason for hiding this comment

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

I think that at the moment we can make all these private.

Copy link
Member

Choose a reason for hiding this comment

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

This is not done yet.

Copy link
Member

@rubenzorrilla rubenzorrilla left a comment

Choose a reason for hiding this comment

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

I finished the review. Overall it's a great contribution, thanks. Nevertheless, there are still some minors that we can easily address.

@AlejandroCornejo
Copy link
Member Author

AlejandroCornejo commented Jan 17, 2025

I finished the review. Overall it's a great contribution, thanks. Nevertheless, there are still some minors that we can easily address.

Let's see if I've addressed them all :)

cover2-3833907614

Comment on lines +201 to +220
switch (integration_order)
{
case 1:
mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1;
break;
case 2:
mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2;
break;
case 3:
mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_3;
break;
case 4:
mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_4;
break;
case 5:
mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_5;
break;
default:
KRATOS_WARNING("SmallDisplacementMixedStrainDisplacementElement") << "Integration order " << integration_order << " is not available, using GI_LOBATTO_1" << std::endl;
mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_LOBATTO_1;
Copy link
Member

Choose a reason for hiding this comment

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

I think it is possible to directly specify the quadrature. I think that it is more convenient in this case to do so rather than the order.

Copy link
Member Author

Choose a reason for hiding this comment

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

How can we do this?

@rubenzorrilla
Copy link
Member

There are still things to be addressed, specially the issue with the variable.

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

Successfully merging this pull request may close these issues.

3 participants