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

[StructuralMechanicsApplication] Make sure the structural/interface elements used in the quay wall work with the geo reset displacement process (LinearTimoshenkoBeamElement2D2N) #12948

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
+180 −45
Open

[StructuralMechanicsApplication] Make sure the structural/interface elements used in the quay wall work with the geo reset displacement process (LinearTimoshenkoBeamElement2D2N) #12948

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
69 changes: 42 additions & 27 deletions ...cturalMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -850,18 +850,11 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
const auto &r_geometry = GetGeometry();

ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
auto &r_cl_options = cl_values.GetOptions();
r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true);
r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);

const double length = CalculateLength();
const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);

// Let's initialize the cl values
VectorType strain_vector(strain_size), stress_vector(strain_size);
strain_vector.clear();
cl_values.SetStrainVector(strain_vector);
cl_values.SetStressVector(stress_vector);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);

VectorType nodal_values(mat_size);
GetNodalValuesVector(nodal_values);

Expand All @@ -880,18 +873,10 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
const auto &r_geometry = GetGeometry();

ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
auto &r_cl_options = cl_values.GetOptions();
r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true);
r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);

const double length = CalculateLength();
const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);

// Let's initialize the cl values
VectorType strain_vector(strain_size), stress_vector(strain_size);
strain_vector.clear();
cl_values.SetStrainVector(strain_vector);
cl_values.SetStressVector(stress_vector);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);
VectorType nodal_values(mat_size);
GetNodalValuesVector(nodal_values);

Expand All @@ -910,18 +895,10 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
const auto &r_geometry = GetGeometry();

ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
auto &r_cl_options = cl_values.GetOptions();
r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true);
r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);

const double length = CalculateLength();
const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);

// Let's initialize the cl values
VectorType strain_vector(strain_size), stress_vector(strain_size);
strain_vector.clear();
cl_values.SetStrainVector(strain_vector);
cl_values.SetStressVector(stress_vector);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);
VectorType nodal_values(mat_size);
GetNodalValuesVector(nodal_values);

Expand Down Expand Up @@ -964,6 +941,44 @@ void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
}
}

void LinearTimoshenkoBeamElement2D2N::CalculateOnIntegrationPoints(
const Variable<Vector>& rVariable,
std::vector<Vector>& rOutput,
const ProcessInfo& rProcessInfo
)
{
const auto& integration_points = IntegrationPoints(GetIntegrationMethod());
const SizeType strain_size = mConstitutiveLawVector[0]->GetStrainSize();
const SizeType mat_size = GetDoFsPerNode() * GetGeometry().size();
rOutput.resize(integration_points.size());
const auto &r_props = GetProperties();

if (rVariable == PK2_STRESS_VECTOR) {

const auto &r_geometry = GetGeometry();

ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo);
const double length = CalculateLength();
const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length);

VectorType strain_vector(strain_size), stress_vector(strain_size);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector);

VectorType nodal_values(mat_size);
GetNodalValuesVector(nodal_values);

// Loop over the integration points
for (SizeType integration_point = 0; integration_point < integration_points.size(); ++integration_point) {
CalculateGeneralizedStrainsVector(strain_vector, length, Phi, integration_points[integration_point].X(), nodal_values);
mConstitutiveLawVector[integration_point]->CalculateMaterialResponsePK2(cl_values);
auto stress_vector = cl_values.GetStressVector();
if ( this->GetProperties().Has(TIMOSHENKO_BEAM_PRESTRESS_PK2)) {
stress_vector += this->GetProperties()[TIMOSHENKO_BEAM_PRESTRESS_PK2];
}
rOutput[integration_point] = stress_vector;
}
}
}
/***********************************************************************************/
/***********************************************************************************/

Expand Down
12 changes: 12 additions & 0 deletions ...ructuralMechanicsApplication/custom_elements/beam_elements/timoshenko_beam_element_2D2N.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -387,6 +387,18 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTimoshenkoBeamElement2D
const ProcessInfo& rCurrentProcessInfo
) override;

/**
* @brief Calculate a double Variable on the Element Constitutive Law
* @param rVariable The variable we want to get
* @param rOutput The values obtained in the integration points
* @param rCurrentProcessInfo the current process info instance
*/
void CalculateOnIntegrationPoints(
const Variable<Vector>& rVariable,
std::vector<Vector>& rOutput,
const ProcessInfo& rCurrentProcessInfo
) override;

/**
* @brief Get on rVariable Constitutive Law from the element
* @param rVariable The variable we want to get
Expand Down
19 changes: 3 additions & 16 deletions ...ns/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -733,7 +733,7 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
VectorType strain_vector(1), stress_vector(1);
MatrixType C(1,1);
InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, C);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C);

const double length = CalculateLength();
SystemSizeBoundedArrayType nodal_values(SystemSize);
Expand All @@ -756,7 +756,7 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
VectorType strain_vector(1), stress_vector(1);
MatrixType C(1,1);
InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, C);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C);

const double length = CalculateLength();
SystemSizeBoundedArrayType nodal_values(SystemSize);
Expand All @@ -773,19 +773,6 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
}
}

template <SizeType TDimension, SizeType TNNodes>
void LinearTrussElement<TDimension, TNNodes>::InitializeConstitutiveLawValues(ConstitutiveLaw::Parameters& rValues,
VectorType& rStrainVector, VectorType& rStressVector, MatrixType rConstitutiveMatrix)
{
rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS , true);
rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);

rStrainVector.clear();
rValues.SetStrainVector(rStrainVector);
rValues.SetStressVector(rStressVector);
rValues.SetConstitutiveMatrix(rConstitutiveMatrix);
}

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

Expand All @@ -803,7 +790,7 @@ void LinearTrussElement<TDimension, TNNodes>::CalculateOnIntegrationPoints(
ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo);
VectorType strain_vector(1), stress_vector(1);
MatrixType constitutive_matrix(1,1);
InitializeConstitutiveLawValues(cl_values, strain_vector, stress_vector, constitutive_matrix);
StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, constitutive_matrix);

const double length = CalculateLength();
SystemSizeBoundedArrayType nodal_values(SystemSize);
Expand Down
2 changes: 0 additions & 2 deletions ...ions/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -477,8 +477,6 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement
///@}
///@name Private Operations
///@{
void InitializeConstitutiveLawValues(ConstitutiveLaw::Parameters& rValues,
VectorType& rStrainVector, VectorType& rStressVector, MatrixType rConstitutiveMatrix);

///@}
///@name Private Access
Expand Down
1 change: 1 addition & 0 deletions .../StructuralMechanicsApplication/custom_python/structural_mechanics_python_application.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,6 +91,7 @@ PYBIND11_MODULE(KratosStructuralMechanicsApplication,m)
KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,I33)
KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,IZ)
KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,BEAM_INITIAL_STRAIN_VECTOR)
KRATOS_REGISTER_IN_PYTHON_VARIABLE(m,TIMOSHENKO_BEAM_PRESTRESS_PK2)

// semi rigid beam variables
KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, ROTATIONAL_STIFFNESS_AXIS_2)
Expand Down
18 changes: 18 additions & 0 deletions ...tructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -635,6 +635,24 @@ double CalculatePhi(const Properties& rProperties, const double L)
/***********************************************************************************/
/***********************************************************************************/

void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
Vector& rStrainVector, Vector& rStressVector, Matrix& rConstitutiveMatrix)
{
InitializeConstitutiveLawValuesForStressCalculation(rValues, rStrainVector, rStressVector);
rValues.SetConstitutiveMatrix(rConstitutiveMatrix);
}

void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
Vector& rStrainVector, Vector& rStressVector)
{
rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_STRESS , true);
rValues.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, false);

rStrainVector.clear();
rValues.SetStrainVector(rStrainVector);
rValues.SetStressVector(rStressVector);
}

} // namespace StructuralMechanicsElementUtilities.
} // namespace Kratos.

Expand Down
6 changes: 6 additions & 0 deletions .../StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -350,5 +350,11 @@ double GetReferenceRotationAngle2D3NBeam(const GeometryType &rGeometry);
*/
KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) double CalculatePhi(const Properties& rProperties, const double L);

void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
Vector& rStrainVector, Vector& rStressVector, Matrix& rConstitutiveMatrix);

void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parameters& rValues,
Vector& rStrainVector, Vector& rStressVector);

} // namespace StructuralMechanicsElementUtilities.
} // namespace Kratos.
2 changes: 2 additions & 0 deletions applications/StructuralMechanicsApplication/structural_mechanics_application.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -347,6 +347,8 @@ void KratosStructuralMechanicsApplication::Register() {
KRATOS_REGISTER_VARIABLE(I33)
KRATOS_REGISTER_VARIABLE(LUMPED_MASS_ROTATION_COEFFICIENT)
KRATOS_REGISTER_VARIABLE(BEAM_INITIAL_STRAIN_VECTOR)
KRATOS_REGISTER_VARIABLE(TIMOSHENKO_BEAM_PRESTRESS_PK2)


// semi rigid beam variables
KRATOS_REGISTER_VARIABLE(ROTATIONAL_STIFFNESS_AXIS_2)
Expand Down
2 changes: 2 additions & 0 deletions applications/StructuralMechanicsApplication/structural_mechanics_application_variables.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,8 @@ KRATOS_CREATE_VARIABLE(double, I22)
KRATOS_CREATE_VARIABLE(double, I33)
KRATOS_CREATE_VARIABLE(double, LUMPED_MASS_ROTATION_COEFFICIENT)
KRATOS_CREATE_VARIABLE(Vector, BEAM_INITIAL_STRAIN_VECTOR)
KRATOS_CREATE_VARIABLE(Vector, TIMOSHENKO_BEAM_PRESTRESS_PK2)


// semi rigid beam variables
KRATOS_CREATE_VARIABLE(double, ROTATIONAL_STIFFNESS_AXIS_2)
Expand Down
1 change: 1 addition & 0 deletions applications/StructuralMechanicsApplication/structural_mechanics_application_variables.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@ namespace Kratos
KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, I33)
KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,double, LUMPED_MASS_ROTATION_COEFFICIENT)
KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, BEAM_INITIAL_STRAIN_VECTOR)
KRATOS_DEFINE_APPLICATION_VARIABLE( STRUCTURAL_MECHANICS_APPLICATION,Vector, TIMOSHENKO_BEAM_PRESTRESS_PK2)


// Semi rigid beam variables
Expand Down
93 changes: 93 additions & 0 deletions applications/StructuralMechanicsApplication/tests/cpp_tests/test_timoshenko_beams.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,93 @@
// KRATOS ___| | | |
// \___ \ __| __| | | __| __| | | __| _` | |
// | | | | | ( | | | | ( | |
// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS
//
// License: BSD License
// license: StructuralMechanicsApplication/license.txt
//
// Main authors: Gennady Markelov
//

// Project includes
#include "containers/model.h"
#include "structural_mechanics_fast_suite.h"
#include "structural_mechanics_application_variables.h"

#include <utility>
#include <boost/numeric/ublas/assignment.hpp>

namespace Kratos::Testing
{

void CreateBeamModel2N(std::string TimoshenkoBeamElementName)
{
Model current_model;
auto &r_model_part = current_model.CreateModelPart("ModelPart",1);
r_model_part.GetProcessInfo().SetValue(DOMAIN_SIZE, 3);
r_model_part.AddNodalSolutionStepVariable(DISPLACEMENT);
r_model_part.AddNodalSolutionStepVariable(ROTATION_Z);

// Set the element properties
auto p_elem_prop = r_model_part.CreateNewProperties(0);
constexpr auto youngs_modulus = 2.0e+06;
p_elem_prop->SetValue(YOUNG_MODULUS, youngs_modulus);
constexpr auto cross_area = 1.0;
p_elem_prop->SetValue(CROSS_AREA, cross_area);
constexpr auto i33 = 1.0;
p_elem_prop->SetValue(I33, i33);
constexpr auto area_effective_y = 1.0;
p_elem_prop->SetValue(AREA_EFFECTIVE_Y, area_effective_y);

const auto &r_clone_cl = KratosComponents<ConstitutiveLaw>::Get("TimoshenkoBeamElasticConstitutiveLaw");
p_elem_prop->SetValue(CONSTITUTIVE_LAW, r_clone_cl.Clone());

// Create the test element
constexpr double directional_length = 2.0;
auto p_node_1 = r_model_part.CreateNewNode(1, 0.0, 0.0, 0.0);
auto p_node_2 = r_model_part.CreateNewNode(2, directional_length, directional_length, directional_length);

for (auto& r_node : r_model_part.Nodes()){
r_node.AddDof(DISPLACEMENT_X);
r_node.AddDof(DISPLACEMENT_Y);
r_node.AddDof(DISPLACEMENT_Z);
r_node.AddDof(ROTATION_Z);
}

std::vector<ModelPart::IndexType> element_nodes {1,2};
auto p_element = r_model_part.CreateNewElement(std::move(TimoshenkoBeamElementName), 1, element_nodes, p_elem_prop);
const auto& r_process_info = r_model_part.GetProcessInfo();
p_element->Initialize(r_process_info); // Initialize the element to initialize the constitutive law

constexpr auto induced_strain = 0.1;
p_element->GetGeometry()[1].FastGetSolutionStepValue(DISPLACEMENT) += ScalarVector(3, induced_strain * directional_length);
p_element->GetGeometry()[1].FastGetSolutionStepValue(ROTATION_Z) += 0.1;

std::vector<Vector> stress_vectors;
p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vectors, r_process_info);

constexpr auto expected_stress = induced_strain * youngs_modulus;
constexpr auto expected_shear_stress = -37500.0;

constexpr auto tolerance = 1.0e-5;
Vector expected_stress_vector(3);
expected_stress_vector <<= expected_stress, 29631.5,expected_shear_stress;
KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[0], tolerance);
expected_stress_vector <<= expected_stress, 70710.7,expected_shear_stress;
KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[1], tolerance);
expected_stress_vector <<= expected_stress, 111790,expected_shear_stress;
KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[2], tolerance);

Vector pre_stress(3);
pre_stress <<= 1.0e5, 1.0e4, 1.0e3;
p_element->GetProperties().SetValue(TIMOSHENKO_BEAM_PRESTRESS_PK2, pre_stress);
p_element->CalculateOnIntegrationPoints(PK2_STRESS_VECTOR, stress_vectors, r_process_info);
expected_stress_vector += pre_stress;
KRATOS_EXPECT_VECTOR_RELATIVE_NEAR(expected_stress_vector, stress_vectors[2], tolerance);
}

KRATOS_TEST_CASE_IN_SUITE(LinearTimoshenkoBeam2D2N_CalculatesPK2Stress, KratosStructuralMechanicsFastSuite)
{
CreateBeamModel2N("LinearTimoshenkoBeamElement2D2N");
}
}
Loading