Skip to content

[Core][StructuralMechanicsApplication] Refactor solid element methods to ensure consistent Jacobian computation for Updated and Total Lagrangian elements #13367

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

Closed
wants to merge 2 commits into from
Closed

[Core][StructuralMechanicsApplication] Refactor solid element methods to ensure consistent Jacobian computation for Updated and Total Lagrangian elements #13367

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
0412861
[StructuralMechanicsApplication] Consistent initial configuration sol…
loumalouomega Apr 23, 2025
ae440c4
Merge branch 'master' into structural/consistent-initial-configuratio…
loumalouomega Apr 24, 2025
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
17 changes: 3 additions & 14 deletions ...ions/StructuralMechanicsApplication/custom_elements/solid_elements/base_solid_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1695,24 +1695,13 @@ void BaseSolidElement::SetConstitutiveVariables(
/***********************************************************************************/
/***********************************************************************************/

Matrix& BaseSolidElement::CalculateDeltaDisplacement(Matrix& DeltaDisplacement) const
Matrix& BaseSolidElement::CalculateDeltaDisplacement(Matrix& rDeltaDisplacement) const
{
KRATOS_TRY

const SizeType number_of_nodes = GetGeometry().PointsNumber();
const SizeType dimension = GetGeometry().WorkingSpaceDimension();

DeltaDisplacement.resize(number_of_nodes , dimension, false);

for ( IndexType i_node = 0; i_node < number_of_nodes; i_node++ ) {
const array_1d<double, 3 >& current_displacement = GetGeometry()[i_node].FastGetSolutionStepValue(DISPLACEMENT);
const array_1d<double, 3 >& previous_displacement = GetGeometry()[i_node].FastGetSolutionStepValue(DISPLACEMENT,1);

for ( IndexType j_dim = 0; j_dim < dimension; ++j_dim )
DeltaDisplacement(i_node, j_dim) = current_displacement[j_dim] - previous_displacement[j_dim];
}
GeometryUtils::CalculateDeltaDisplacement(this->GetGeometry(), rDeltaDisplacement);

return DeltaDisplacement;
return rDeltaDisplacement;

KRATOS_CATCH( "" )
}
Expand Down
29 changes: 13 additions & 16 deletions ...ations/StructuralMechanicsApplication/custom_elements/solid_elements/base_solid_element.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,6 @@
// External includes

// Project includes
#include "includes/define.h"
#include "includes/element.h"
#include "utilities/integration_utilities.h"
#include "structural_mechanics_application_variables.h"
Expand All @@ -36,10 +35,10 @@ namespace Kratos
///@{

/// The definition of the index type
typedef std::size_t IndexType;
using IndexType = std::size_t;

/// The definition of the sizetype
typedef std::size_t SizeType;
using SizeType = std::size_t;

///@}
///@name Enum's
Expand Down Expand Up @@ -139,26 +138,24 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) BaseSolidElement
///@{

///Reference type definition for constitutive laws
typedef ConstitutiveLaw ConstitutiveLawType;
using ConstitutiveLawType = ConstitutiveLaw;

///Pointer type for constitutive laws
typedef ConstitutiveLawType::Pointer ConstitutiveLawPointerType;
using ConstitutiveLawPointerType = ConstitutiveLawType::Pointer;

///StressMeasure from constitutive laws
typedef ConstitutiveLawType::StressMeasure StressMeasureType;
using StressMeasureType = ConstitutiveLawType::StressMeasure;

///Type definition for integration methods
typedef GeometryData::IntegrationMethod IntegrationMethod;

/// This is the definition of the node.
typedef Node NodeType;
using IntegrationMethod = GeometryData::IntegrationMethod;

/// The base element type
typedef Element BaseType;
using BaseType = Element;

// Counted pointer of BaseSolidElement
/// Counted pointer of BaseSolidElement
KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION( BaseSolidElement );

/// The definition of the rotation flag
KRATOS_DEFINE_LOCAL_FLAG(ROTATED);

///@}
Expand Down Expand Up @@ -779,10 +776,10 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) BaseSolidElement

/**
* @brief This methods gives us a matrix with the increment of displacement
* @param DeltaDisplacement The matrix containing the increment of displacements
* @return DeltaDisplacement: The matrix containing the increment of displacements
* @param rDeltaDisplacement The matrix containing the increment of displacements
* @return rDeltaDisplacement: The matrix containing the increment of displacements
*/
Matrix& CalculateDeltaDisplacement(Matrix& DeltaDisplacement) const;
Matrix& CalculateDeltaDisplacement(Matrix& rDeltaDisplacement) const;

/**
* @brief This functions calculate the derivatives in the reference frame
Expand All @@ -793,7 +790,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) BaseSolidElement
* @param ThisIntegrationMethod The integration method considered
* @return The determinant of the jacobian in the reference configuration
*/
virtual double CalculateDerivativesOnReferenceConfiguration(
double CalculateDerivativesOnReferenceConfiguration(
Matrix& rJ0,
Matrix& rInvJ0,
Matrix& rDN_DX,
Expand Down
29 changes: 0 additions & 29 deletions ...ions/StructuralMechanicsApplication/custom_elements/solid_elements/updated_lagrangian.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -339,35 +339,6 @@ void UpdatedLagrangian::CalculateKinematicVariables(
/***********************************************************************************/
/***********************************************************************************/

double UpdatedLagrangian::CalculateDerivativesOnReferenceConfiguration(
Matrix& J0,
Matrix& InvJ0,
Matrix& DN_DX,
const IndexType PointNumber,
IntegrationMethod ThisIntegrationMethod
) const
{
J0.clear();

double detJ0;

Matrix delta_displacement;
delta_displacement = this->CalculateDeltaDisplacement(delta_displacement);

J0 = this->GetGeometry().Jacobian( J0, PointNumber, ThisIntegrationMethod, delta_displacement);

const Matrix& DN_De = this->GetGeometry().ShapeFunctionsLocalGradients(ThisIntegrationMethod)[PointNumber];

MathUtils<double>::InvertMatrix( J0, InvJ0, detJ0 );

noalias( DN_DX ) = prod( DN_De, InvJ0);

return detJ0;
}

/***********************************************************************************/
/***********************************************************************************/

void UpdatedLagrangian::CalculateB(
Matrix& rB,
const Matrix& rDN_DX,
Expand Down
17 changes: 0 additions & 17 deletions ...ations/StructuralMechanicsApplication/custom_elements/solid_elements/updated_lagrangian.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -384,23 +384,6 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) UpdatedLagrangian
const GeometryType::IntegrationMethod& rIntegrationMethod
) override;

/**
* @brief This functions calculate the derivatives in the reference frame
* @param J0 The jacobian in the reference configuration
* @param InvJ0 The inverse of the jacobian in the reference configuration
* @param DN_DX The gradient derivative of the shape function
* @param PointNumber The id of the integration point considered
* @param ThisIntegrationMethod The integration method considered
* @return The determinant of the jacobian in the reference configuration
*/
double CalculateDerivativesOnReferenceConfiguration(
Matrix& J0,
Matrix& InvJ0,
Matrix& DN_DX,
const IndexType PointNumber,
IntegrationMethod ThisIntegrationMethod
) const override;

///@}
///@name Protected Operations
///@{
Expand Down
35 changes: 30 additions & 5 deletions kratos/utilities/geometry_utilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,8 +18,9 @@
// External includes

// Project includes
#include "geometries/geometry.h"
#include "includes/node.h"
#include "geometries/geometry.h"
#include "includes/variables.h"

namespace Kratos
{
Expand Down Expand Up @@ -712,13 +713,37 @@ class KRATOS_API(KRATOS_CORE) GeometryUtils
)
{
Matrix delta_position(rGeom.PointsNumber(), rGeom.WorkingSpaceDimension());
for (std::size_t i = 0; i < rGeom.PointsNumber(); ++i)
for (std::size_t j = 0; j < rGeom.WorkingSpaceDimension(); ++j)
delta_position(i, j) = rGeom[i].Coordinates()[j] -
rGeom[i].GetInitialPosition().Coordinates()[j];
CalculateDeltaDisplacement(rGeom, delta_position);
rGeom.Jacobian(rJ0, rCoords, delta_position);
}

/**
* @brief This function computes the delta displacement of the element nodes.
* @details The delta displacement is defined as the difference between the current and previous displacements of the nodes.
* @param rGeom element geometry.
* @param rDeltaDisplacement delta displacement matrix.
*/
static void CalculateDeltaDisplacement(
const GeometryType& rGeom,
Matrix& rDeltaDisplacement
)
{
const SizeType number_of_nodes = rGeom.PointsNumber();
const SizeType dimension = rGeom.WorkingSpaceDimension();
rDeltaDisplacement.resize(number_of_nodes , dimension, false);

for ( IndexType i_node = 0; i_node < number_of_nodes; i_node++ ) {
auto& r_node = rGeom[i_node];
const array_1d<double, 3>& r_current_displacement = r_node.FastGetSolutionStepValue(DISPLACEMENT);
const array_1d<double, 3>& r_previous_displacement = r_node.FastGetSolutionStepValue(DISPLACEMENT,1);

for ( IndexType j_dim = 0; j_dim < dimension; ++j_dim ) {
rDeltaDisplacement(i_node, j_dim) = r_current_displacement[j_dim] - r_previous_displacement[j_dim];
// rDeltaDisplacement(i_node, j_dim) = r_node.Coordinates()[j_dim] - r_node.GetInitialPosition().Coordinates()[j_dim]; // NOTE: This is problematic when mesh is not moved, for example during NL iterations. Mesh is moved usually at the end of a step.
}
}
}

/**
* @brief Calculate the Jacobian on the initial configuration.
* @param rGeom element geometry.
Expand Down
Loading