From 24b46e6705e09ac9bf756030ab7a531514d19f37 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Sun, 3 Nov 2024 16:44:33 +0100 Subject: [PATCH 01/32] base files added --- .../custom_constitutive/truss_constitutive_law.cpp | 3 +++ .../structural_mechanics_application.h | 1 + 2 files changed, 4 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp b/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp index 595172ec9b74..9451674b12bf 100644 --- a/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp +++ b/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp @@ -156,6 +156,9 @@ double TrussConstitutiveLaw::CalculateStressElastic( CalculateValue(rParameterValues,TANGENT_MODULUS,tangent_modulus); const double current_stress = tangent_modulus*current_strain[0]; + Matrix& r_const_matrix = rParameterValues.GetConstitutiveMatrix(); + r_const_matrix.resize(1, 1); + r_const_matrix(0, 0) = tangent_modulus; return current_stress; } diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.h b/applications/StructuralMechanicsApplication/structural_mechanics_application.h index e677969efa29..ca4b26bd05cd 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.h +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.h @@ -34,6 +34,7 @@ #include "custom_elements/truss_elements/truss_element_3D2N.hpp" #include "custom_elements/truss_elements/truss_element_linear_3D2N.hpp" #include "custom_elements/truss_elements/cable_element_3D2N.hpp" +#include "custom_elements/truss_elements/linear_truss_element_2D.h" /* Adding beam element */ #include "custom_elements/beam_elements/cr_beam_element_3D2N.hpp" From 3f409e76e8373cc32704d87f9e272968bbae967c Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Sun, 3 Nov 2024 16:45:29 +0100 Subject: [PATCH 02/32] add more --- .../linear_truss_element_2D.cpp | 569 ++++++++++++++++++ .../truss_elements/linear_truss_element_2D.h | 486 +++++++++++++++ ...enko_beam_element_2D2N.cpp:Zone.Identifier | 0 ...shenko_beam_element_2D2N.h:Zone.Identifier | 0 4 files changed, 1055 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp create mode 100644 applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h create mode 100644 applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.cpp:Zone.Identifier create mode 100644 applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp new file mode 100644 index 000000000000..859bb17a94e0 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -0,0 +1,569 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Alejandro Cornejo +// +// + +// System includes + +// External includes + +// Project includes + +// Application includes +#include "linear_truss_element_2D.h" +#include "custom_utilities/constitutive_law_utilities.h" +#include "structural_mechanics_application_variables.h" + +namespace Kratos +{ + +template +void LinearTrussElement2D::Initialize(const ProcessInfo& rCurrentProcessInfo) +{ + KRATOS_TRY + + // Initialization should not be done again in a restart! + if (!rCurrentProcessInfo[IS_RESTARTED]) { + if (this->UseGeometryIntegrationMethod()) { + if (GetProperties().Has(INTEGRATION_ORDER) ) { + mThisIntegrationMethod = static_cast(GetProperties()[INTEGRATION_ORDER] - 1); + } else { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_3; + } + } + + const auto& r_integration_points = this->IntegrationPoints(mThisIntegrationMethod); + + // Constitutive Law initialisation + if (mConstitutiveLawVector.size() != r_integration_points.size()) + mConstitutiveLawVector.resize(r_integration_points.size()); + InitializeMaterial(); + } + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::InitializeMaterial() +{ + KRATOS_TRY + + if (GetProperties()[CONSTITUTIVE_LAW] != nullptr) { + const auto& r_geometry = GetGeometry(); + const auto& r_properties = GetProperties(); + auto N_values = Vector(); + for (IndexType point_number = 0; point_number < mConstitutiveLawVector.size(); ++point_number) { + mConstitutiveLawVector[point_number] = GetProperties()[CONSTITUTIVE_LAW]->Clone(); + mConstitutiveLawVector[point_number]->InitializeMaterial(r_properties, r_geometry, N_values); + } + } else + KRATOS_ERROR << "A constitutive law needs to be specified for the element with ID " << this->Id() << std::endl; + + KRATOS_CATCH(""); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +Element::Pointer LinearTrussElement2D::Clone( + IndexType NewId, + NodesArrayType const& rThisNodes + ) const +{ + KRATOS_TRY + + LinearTrussElement2D::Pointer p_new_elem = Kratos::make_intrusive> + (NewId, GetGeometry().Create(rThisNodes), pGetProperties()); + p_new_elem->SetData(this->GetData()); + p_new_elem->Set(Flags(*this)); + + // Currently selected integration methods + p_new_elem->SetIntegrationMethod(mThisIntegrationMethod); + + // The vector containing the constitutive laws + p_new_elem->SetConstitutiveLawVector(mConstitutiveLawVector); + + return p_new_elem; + + KRATOS_CATCH(""); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::EquationIdVector( + EquationIdVectorType& rResult, + const ProcessInfo& rCurrentProcessInfo + ) const +{ + const auto& r_geometry = this->GetGeometry(); + const SizeType number_of_nodes = r_geometry.size(); + const SizeType dofs_per_node = GetDoFsPerNode(); // u, v, theta + + IndexType local_index = 0; + + if (rResult.size() != dofs_per_node * number_of_nodes) + rResult.resize(dofs_per_node * number_of_nodes, false); + + const IndexType xpos = this->GetGeometry()[0].GetDofPosition(DISPLACEMENT_X); + const IndexType rot_pos = this->GetGeometry()[0].GetDofPosition(ROTATION_Z); + + for (IndexType i = 0; i < number_of_nodes; ++i) { + rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_X, xpos ).EquationId(); + rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_Y, xpos + 1).EquationId(); + rResult[local_index++] = r_geometry[i].GetDof(ROTATION_Z , rot_pos ).EquationId(); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::GetDofList( + DofsVectorType& rElementalDofList, + const ProcessInfo& rCurrentProcessInfo + ) const +{ + KRATOS_TRY; + + const auto& r_geom = GetGeometry(); + const SizeType number_of_nodes = r_geom.size(); + const SizeType dofs_per_node = GetDoFsPerNode(); // u, v, theta + rElementalDofList.resize(dofs_per_node * number_of_nodes); + + for (IndexType i = 0; i < number_of_nodes; ++i) { + const SizeType index = i * dofs_per_node; + rElementalDofList[index] = r_geom[i].pGetDof(DISPLACEMENT_X); + rElementalDofList[index + 1] = r_geom[i].pGetDof(DISPLACEMENT_Y); + rElementalDofList[index + 2] = r_geom[i].pGetDof(ROTATION_Z ); + } + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::GetShapeFunctionsValues( + VectorType& rN, + const double Length, + const double Phi, + const double xi + ) const +{ + if (rN.size() != 4) + rN.resize(4, false); + const double one_plus_phi = 1.0 + Phi; + const double xi_square = xi * xi; + rN[0] = (xi - 1.0) * (xi + xi_square - 2.0 * one_plus_phi) / (4.0 * one_plus_phi); + rN[1] = (1.0 - xi_square) * (1.0 - xi + Phi) * Length / (8.0 * one_plus_phi); + rN[2] = (1.0 + xi) * (xi - xi_square + 2.0 * one_plus_phi) / (4.0 * one_plus_phi); + rN[3] = (xi_square - 1.0) * (1.0 + xi + Phi) * Length / (8.0 * one_plus_phi); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( + VectorType& rN, + const double Length, + const double Phi, + const double xi + ) const +{ + if (rN.size() != 4) + rN.resize(4, false); + const double one_plus_phi = 1.0 + Phi; + const double xi_square = xi * xi; + rN[0] = (-6.0 + 6.0 * xi_square - 4.0 * Phi) / (4.0 * one_plus_phi * Length); + rN[1] = (-1.0 + 3.0 * xi_square - 2.0 * xi * one_plus_phi) / (4.0 * one_plus_phi); + rN[2] = (6.0 - 6.0 * xi_square + 4.0 * Phi) / (4.0 * one_plus_phi * Length); + rN[3] = (-1.0 + 3.0 * xi_square + 2.0 * xi * one_plus_phi) / (4.0 * one_plus_phi); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::GetNodalValuesVector(VectorType& rNodalValues) const +{ + // if (rNodalValues.size() != 6) + // rNodalValues.resize(6, false); + // const auto &r_geom = GetGeometry(); + + // const double angle = GetAngle(); + + // if (std::abs(angle) > std::numeric_limits::epsilon()) { + // BoundedMatrix T; + // BoundedVector global_values; + // BoundedMatrix global_size_T; + // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); + // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); + + // const auto &r_displ_0 = r_geom[0].FastGetSolutionStepValue(DISPLACEMENT); + // global_values[0] = r_displ_0[0]; + // global_values[1] = r_displ_0[1]; + // global_values[2] = r_geom[0].FastGetSolutionStepValue(ROTATION_Z); + + // const auto& r_displ_1 = r_geom[1].FastGetSolutionStepValue(DISPLACEMENT); + // global_values[3] = r_displ_1[0]; + // global_values[4] = r_displ_1[1]; + // global_values[5] = r_geom[1].FastGetSolutionStepValue(ROTATION_Z); + + // // We rotate to local axes + // noalias(rNodalValues) = prod(trans(global_size_T), global_values); + + // } else { + // const auto &r_displ_0 = r_geom[0].FastGetSolutionStepValue(DISPLACEMENT); + // rNodalValues[0] = r_displ_0[0]; + // rNodalValues[1] = r_displ_0[1]; + // rNodalValues[2] = r_geom[0].FastGetSolutionStepValue(ROTATION_Z); + + // const auto& r_displ_1 = r_geom[1].FastGetSolutionStepValue(DISPLACEMENT); + // rNodalValues[3] = r_displ_1[0]; + // rNodalValues[4] = r_displ_1[1]; + // rNodalValues[5] = r_geom[1].FastGetSolutionStepValue(ROTATION_Z); + // } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +array_1d LinearTrussElement2D::GetLocalAxesBodyForce( + const Element &rElement, + const GeometryType::IntegrationPointsArrayType &rIntegrationPoints, + const IndexType PointNumber + ) const +{ + const double angle = GetAngle(); + const auto body_force = StructuralMechanicsElementUtilities::GetBodyForce(*this, rIntegrationPoints, PointNumber); + + const double c = std::cos(angle); + const double s = std::sin(angle); + array_1d local_body_force = ZeroVector(3); + local_body_force[0] = c * body_force[0] + s * body_force[1]; + local_body_force[1] = -s * body_force[0] + c * body_force[1]; + return local_body_force; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::CalculateLocalSystem( + MatrixType& rLHS, + VectorType& rRHS, + const ProcessInfo& rProcessInfo + ) +{ + // KRATOS_TRY; + // const auto &r_props = GetProperties(); + // const auto &r_geometry = GetGeometry(); + // const SizeType number_of_nodes = r_geometry.size(); + // const SizeType mat_size = GetDoFsPerNode() * number_of_nodes; + + // if (rLHS.size1() != mat_size || rLHS.size2() != mat_size) { + // rLHS.resize(mat_size, mat_size, false); + // } + // noalias(rLHS) = ZeroMatrix(mat_size, mat_size); + + // if (rRHS.size() != mat_size) { + // rRHS.resize(mat_size, false); + // } + // noalias(rRHS) = ZeroVector(mat_size); + + // const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + + // 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, true); + + // const double length = CalculateLength(); + // const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length); + // const double J = 0.5 * length; + // const double area = r_props[CROSS_AREA]; + + // // Let's initialize the cl values + // VectorType strain_vector(3), stress_vector(3); + // MatrixType constitutive_matrix(3, 3); + // strain_vector.clear(); + // cl_values.SetStrainVector(strain_vector); + // cl_values.SetStressVector(stress_vector); + // cl_values.SetConstitutiveMatrix(constitutive_matrix); + // VectorType nodal_values(mat_size); + // GetNodalValuesVector(nodal_values); + // VectorType global_size_N_2(mat_size), global_size_N(mat_size), N_u_derivatives(number_of_nodes), + // N_theta_derivatives(mat_size-number_of_nodes), N_theta(mat_size-number_of_nodes), N_derivatives(mat_size-number_of_nodes), + // N_u(number_of_nodes), N_shape(mat_size-number_of_nodes), N_s(mat_size-number_of_nodes); + + // // Loop over the integration points + // for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + // const auto local_body_forces = GetLocalAxesBodyForce(*this, integration_points, IP); + + // global_size_N.clear(); + // const double xi = integration_points[IP].X(); + // const double weight = integration_points[IP].Weight(); + // const double jacobian_weight = weight * J; + + // CalculateGeneralizedStrainsVector(strain_vector, length, Phi, xi, nodal_values); + + // mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); + // const Vector &r_generalized_stresses = cl_values.GetStressVector(); + // const double N = r_generalized_stresses[0]; + // const double M = r_generalized_stresses[1]; + // const double V = r_generalized_stresses[2]; + + // const MatrixType& r_constitutive_matrix = cl_values.GetConstitutiveMatrix(); + // const double dN_dEl = r_constitutive_matrix(0, 0); + // const double dM_dkappa = r_constitutive_matrix(1, 1); + // const double dV_dgamma = r_constitutive_matrix(2, 2); + + // GetFirstDerivativesNu0ShapeFunctionsValues(N_u_derivatives, length, Phi, xi); + // GetFirstDerivativesNThetaShapeFunctionsValues(N_theta_derivatives, length, Phi, xi); + // GetNThetaShapeFunctionsValues(N_theta, length, Phi, xi); + // GetFirstDerivativesShapeFunctionsValues(N_derivatives, length, Phi, xi); + // GetShapeFunctionsValues(N_shape, length, Phi, xi); + // GetNu0ShapeFunctionsValues(N_u, length, Phi, xi); + // noalias(N_s) = N_derivatives - N_theta; + + // // Axial contributions + // GlobalSizeAxialVector(global_size_N, N_u_derivatives); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N) * dN_dEl * jacobian_weight; + // noalias(rRHS) -= global_size_N * N * jacobian_weight; + + // // In here we add the cross terms + // GlobalSizeVector(global_size_N_2, N_theta_derivatives); + // const double dN_dkappa = r_constitutive_matrix(0, 1); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dN_dkappa * jacobian_weight; + + // GlobalSizeVector(global_size_N_2, N_s); + // const double dN_dgamma = r_constitutive_matrix(0, 2); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dN_dgamma * jacobian_weight; + + // // Bending contributions + // GlobalSizeVector(global_size_N, N_theta_derivatives); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N) * dM_dkappa * jacobian_weight; + // noalias(rRHS) -= global_size_N * M * jacobian_weight; + + // // In here we add the cross terms + // GlobalSizeAxialVector(global_size_N_2, N_u_derivatives); + // const double dM_dEl = r_constitutive_matrix(1, 0); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dM_dEl * jacobian_weight; + + // GlobalSizeVector(global_size_N_2, N_s); + // const double dM_dgamma = r_constitutive_matrix(1, 2); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dM_dgamma * jacobian_weight; + + // // Shear contributions + // GlobalSizeVector(global_size_N, N_s); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N) * dV_dgamma * jacobian_weight; + // noalias(rRHS) -= global_size_N * V * jacobian_weight; + + // // In here we add the cross terms + // GlobalSizeAxialVector(global_size_N_2, N_u_derivatives); + // const double dV_dEl = r_constitutive_matrix(2, 0); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dV_dEl * jacobian_weight; + + // GlobalSizeVector(global_size_N_2, N_theta_derivatives); + // const double dV_dkappa = r_constitutive_matrix(2, 1); + // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dV_dkappa * jacobian_weight; + + + // // Now we add the body forces contributions + // GlobalSizeAxialVector(global_size_N, N_u); + // noalias(rRHS) += global_size_N * local_body_forces[0] * jacobian_weight * area; + // GlobalSizeVector(global_size_N, N_shape); + // noalias(rRHS) += global_size_N * local_body_forces[1] * jacobian_weight * area; + + // } + // RotateAll(rLHS, rRHS, r_geometry); + + // KRATOS_CATCH(""); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::CalculateLeftHandSide( + MatrixType& rLHS, + const ProcessInfo& rProcessInfo + ) +{ + +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::CalculateRightHandSide( + VectorType& rRHS, + const ProcessInfo& rProcessInfo + ) +{ + +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::RotateLHS( + MatrixType& rLHS, + const GeometryType& rGeometry +) +{ + // const double angle = GetAngle(); + + // if (std::abs(angle) > std::numeric_limits::epsilon()) { + // BoundedMatrix T, Tt; + // BoundedMatrix global_size_T, aux_product; + // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); + // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); + // noalias(aux_product) = prod(rLHS, trans(global_size_T)); + // noalias(rLHS) = prod(global_size_T, aux_product); + // } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::RotateRHS( + VectorType& rRHS, + const GeometryType& rGeometry +) +{ + // const double angle = GetAngle(); + // if (std::abs(angle) > std::numeric_limits::epsilon()) { + // BoundedMatrix T; + // BoundedMatrix global_size_T; + // BoundedVector local_rhs; + // noalias(local_rhs) = rRHS; + // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); + // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); + + // noalias(rRHS) = prod(global_size_T, local_rhs); + // } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::RotateAll( + MatrixType& rLHS, + VectorType& rRHS, + const GeometryType& rGeometry +) +{ + // const double angle = GetAngle(); + // if (std::abs(angle) > std::numeric_limits::epsilon()) { + // BoundedMatrix T; + // BoundedMatrix global_size_T, aux_product; + // BoundedVector local_rhs; + // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); + // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); + + // noalias(local_rhs) = rRHS; + // noalias(rRHS) = prod(global_size_T, local_rhs); + + // noalias(aux_product) = prod(rLHS, trans(global_size_T)); + // noalias(rLHS) = prod(global_size_T, aux_product); + // } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rProcessInfo + ) +{ + +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo + ) +{ + if (rVariable == CONSTITUTIVE_LAW) { + const SizeType integration_points_number = mConstitutiveLawVector.size(); + if (rValues.size() != integration_points_number) { + rValues.resize(integration_points_number); + } + for (IndexType point_number = 0; point_number < integration_points_number; ++point_number) { + rValues[point_number] = mConstitutiveLawVector[point_number]; + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +int LinearTrussElement2D::Check(const ProcessInfo& rCurrentProcessInfo) const +{ + KRATOS_TRY; + + return mConstitutiveLawVector[0]->Check(GetProperties(), GetGeometry(), rCurrentProcessInfo); + + KRATOS_CATCH( "" ); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::save(Serializer& rSerializer) const +{ + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Element); + int IntMethod = int(this->GetIntegrationMethod()); + rSerializer.save("IntegrationMethod",IntMethod); + rSerializer.save("ConstitutiveLawVector", mConstitutiveLawVector); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::load(Serializer& rSerializer) +{ + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element); + int IntMethod; + rSerializer.load("IntegrationMethod",IntMethod); + mThisIntegrationMethod = IntegrationMethod(IntMethod); + rSerializer.load("ConstitutiveLawVector", mConstitutiveLawVector); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template class LinearTrussElement2D<2>; +template class LinearTrussElement2D<3>; + +} // Namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h new file mode 100644 index 000000000000..4f6ae346957c --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -0,0 +1,486 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Alejandro Cornejo +// +// + +#pragma once + +// System includes + +// External includes + +// Project includes +#include "includes/element.h" +#include "custom_utilities/structural_mechanics_element_utilities.h" + +namespace Kratos +{ + +///@name Kratos Globals +///@{ + +///@} +///@name Type Definitions +///@{ + +using SizeType = std::size_t; + +///@} +///@name Enum's +///@{ + +///@} +///@name Functions +///@{ + +///@} +///@name Kratos Classes +///@{ + +/** + * @class LinearTrussElement2D + * @ingroup StructuralMechanicsApplication + * @brief This is the Linear TRUSS element of 2 nodes. + * 2 and 3 noded elements + * @author Alejandro Cornejo + */ +template +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D + : public Element +{ + +public: + + ///@name Type Definitions + ///@{ + + /// The base element type + using BaseType = Element; + static constexpr SizeType DofsPerNode = 2; + static constexpr SizeType SystemSize = DofsPerNode * TNNodes; + using SystemSizeBoundedArrayType = array_1d; + + // Counted pointer of BaseSolidElement + KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(LinearTrussElement2D); + + ///@} + ///@name Life Cycle + ///@{ + + // Constructor void + LinearTrussElement2D() + { + } + + // Constructor using an array of nodes + LinearTrussElement2D(IndexType NewId, GeometryType::Pointer pGeometry) : Element(NewId, pGeometry) + { + if constexpr (TNNodes == 2) { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; + } else { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; + } + } + + // Constructor using an array of nodes with properties + LinearTrussElement2D(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) + : Element(NewId,pGeometry,pProperties) + { + // This is needed to prevent uninitialised integration method in inactive elements + if constexpr (TNNodes == 2) { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; + } else { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; + } + } + + // Copy constructor + LinearTrussElement2D(LinearTrussElement2D const& rOther) + : BaseType(rOther), + mThisIntegrationMethod(rOther.mThisIntegrationMethod), + mConstitutiveLawVector(rOther.mConstitutiveLawVector) + { + } + + // Create method + Element::Pointer Create(IndexType NewId, NodesArrayType const& ThisNodes, PropertiesType::Pointer pProperties) const override + { + return Kratos::make_intrusive(NewId, GetGeometry().Create(ThisNodes), pProperties); + } + + // Create method + Element::Pointer Create( IndexType NewId, GeometryType::Pointer pGeom, PropertiesType::Pointer pProperties ) const override + { + return Kratos::make_intrusive(NewId, pGeom, pProperties); + } + + ///@} + ///@name Operators + ///@{ + + ///@} + ///@name Operations + ///@{ + + /** + * @brief Indicates the amount of DoFs per node (u0, v, theta) + */ + IndexType GetDoFsPerNode() const + { + return 2; + } + + /** + * @brief This method returns the angle of the FE axis + */ + double GetAngle() const + { + // return StructuralMechanicsElementUtilities::GetReferenceRotationAngle2D2NBeam(GetGeometry()); + return 1.0; + } + + /** + * @brief Returns a 6 component vector including the values of the DoFs + * in LOCAL beam axes + */ + void GetNodalValuesVector(VectorType& rNodalValue) const; + + /** + * @brief Computes the length of the FE and returns it + */ + double CalculateLength() const + { + return StructuralMechanicsElementUtilities::CalculateReferenceLength2D2N(*this); + } + + /** + * @brief Called to initialize the element. + * @warning Must be called before any calculation is done + */ + void Initialize(const ProcessInfo& rCurrentProcessInfo) override; + + /** + * @brief It creates a new element pointer and clones the previous element data + * @param NewId the ID of the new element + * @param ThisNodes the nodes of the new element + * @param pProperties the properties assigned to the new element + * @return a Pointer to the new element + */ + Element::Pointer Clone( + IndexType NewId, + NodesArrayType const& rThisNodes + ) const override; + + /** + * @brief Sets on rResult the ID's of the element degrees of freedom + * @param rResult The vector containing the equation id + * @param rCurrentProcessInfo The current process info instance + */ + void EquationIdVector( + EquationIdVectorType& rResult, + const ProcessInfo& rCurrentProcessInfo + ) const override; + + /** + * @brief Sets on rElementalDofList the degrees of freedom of the considered element geometry + * @param rElementalDofList The vector containing the dof of the element + * @param rCurrentProcessInfo The current process info instance + */ + void GetDofList( + DofsVectorType& rElementalDofList, + const ProcessInfo& rCurrentProcessInfo + ) const override; + + /** + * @brief Returns the used integration method + * @return default integration method of the used Geometry + */ + IntegrationMethod GetIntegrationMethod() const override + { + return mThisIntegrationMethod; + } + + /** + * element can be integrated using the GP provided by the geometry or custom ones + * by default, the base element will use the standard integration provided by the geom + * @return bool to select if use/not use GPs given by the geometry + */ + bool UseGeometryIntegrationMethod() const + { + return true; + } + + /** + * @brief Returns the set of integration points + */ + const GeometryType::IntegrationPointsArrayType IntegrationPoints() const + { + return GetGeometry().IntegrationPoints(); + } + + /** + * @brief Returns the set of integration points + */ + const GeometryType::IntegrationPointsArrayType IntegrationPoints(IntegrationMethod ThisMethod) const + { + return GetGeometry().IntegrationPoints(ThisMethod); + } + + /** + * @brief This function returns the 4 shape functions used for interpolating the transverse displacement v. (denoted as N) + * Also its derivatives + * @param rN reference to the shape functions (or derivatives) + * @param Length The size of the beam element + * @param Phi The shear slenderness parameter + * @param xi The coordinate in the natural axes + */ + void GetShapeFunctionsValues(VectorType& rN, const double Length, const double Phi, const double xi) const; + void GetFirstDerivativesShapeFunctionsValues(VectorType& rN, const double Length, const double Phi, const double xi) const; + + /** + * @brief This function rotates the LHS from local to global coordinates + * @param rLHS the left hand side + * @param rGeometry the geometry of the FE + */ + virtual void RotateLHS( + MatrixType &rLHS, + const GeometryType &rGeometry); + + /** + * @brief This function rotates the RHS from local to global coordinates + * @param rRHS the right hand side + * @param rGeometry the geometry of the FE + */ + virtual void RotateRHS( + VectorType &rRHS, + const GeometryType &rGeometry); + + /** + * @brief This function rotates the LHS and RHS from local to global coordinates + * @param rLHS the left hand side + * @param rRHS the right hand side + * @param rGeometry the geometry of the FE + */ + virtual void RotateAll( + MatrixType &rLHS, + VectorType &rRHS, + const GeometryType &rGeometry); + + /** + * @brief This function retrieves the body forces in local axes + * @param rElement the element reference + * @param rIntegrationPoints array of IP + * @param PointNumber tthe IP to be evaluated + */ + array_1d GetLocalAxesBodyForce( + const Element &rElement, + const GeometryType::IntegrationPointsArrayType &rIntegrationPoints, + const IndexType PointNumber) const; + + /** + * @brief This function provides a more general interface to the element. + * @details It is designed so that rLHSvariables and rRHSvariables are passed to the element thus telling what is the desired output + * @param rLeftHandSideMatrix container with the output Left Hand Side matrix + * @param rRightHandSideVector container for the desired RHS output + * @param rCurrentProcessInfo the current process info instance + */ + void CalculateLocalSystem( + MatrixType& rLeftHandSideMatrix, + VectorType& rRightHandSideVector, + const ProcessInfo& rCurrentProcessInfo + ) override; + + /** + * @brief This is called during the assembling process in order to calculate the elemental left hand side matrix only + * @param rLeftHandSideMatrix the elemental left hand side matrix + * @param rCurrentProcessInfo the current process info instance + */ + void CalculateLeftHandSide( + MatrixType& rLeftHandSideMatrix, + const ProcessInfo& rCurrentProcessInfo + ) override; + + /** + * @brief This is called during the assembling process in order to calculate the elemental right hand side vector only + * @param rRightHandSideVector the elemental right hand side vector + * @param rCurrentProcessInfo the current process info instance + */ + void CalculateRightHandSide( + VectorType& rRightHandSideVector, + 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& rVariable, + std::vector& rOutput, + const ProcessInfo& rCurrentProcessInfo + ) override; + + /** + * @brief Get on rVariable Constitutive Law from the element + * @param rVariable The variable we want to get + * @param rValues The results in the integration points + * @param rCurrentProcessInfo the current process info instance + */ + void CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rValues, + const ProcessInfo& rCurrentProcessInfo + ) override; + + /** + * @brief This function provides the place to perform checks on the completeness of the input. + * @details It is designed to be called only once (or anyway, not often) typically at the beginning + * of the calculations, so to verify that nothing is missing from the input + * or that no common error is found. + * @param rCurrentProcessInfo the current process info instance + */ + int Check(const ProcessInfo &rCurrentProcessInfo) const override; + + ///@} + ///@name Access + ///@{ + + + ///@} + ///@name Inquiry + ///@{ + + + ///@} + ///@name Input and output + ///@{ + + /// Print information about this object. + void PrintInfo(std::ostream& rOStream) const override + { + rOStream << "Truss 2D Element #" << Id() << "\nConstitutive law: " << mConstitutiveLawVector[0]->Info(); + } + + /// Print object's data. + void PrintData(std::ostream& rOStream) const override + { + pGetGeometry()->PrintData(rOStream); + } + + ///@} + ///@name Friends + ///@{ + +protected: + + ///@name Protected static Member Variables + ///@{ + + ///@} + ///@name Protected member Variables + ///@{ + + IntegrationMethod mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; /// Currently selected integration methods + + std::vector mConstitutiveLawVector; /// The vector containing the constitutive laws + + ///@} + ///@name Protected Operators + ///@{ + + ///@} + ///@name Protected Operations + ///@{ + + /** + * @brief Sets the used integration method + * @param ThisIntegrationMethod Integration method used + */ + void SetIntegrationMethod(const IntegrationMethod& rThisIntegrationMethod) + { + mThisIntegrationMethod = rThisIntegrationMethod; + } + + /** + * @brief Sets the used constitutive laws + * @param ThisConstitutiveLawVector Constitutive laws used + */ + void SetConstitutiveLawVector(const std::vector& rThisConstitutiveLawVector) + { + mConstitutiveLawVector = rThisConstitutiveLawVector; + } + + /** + * @brief It initializes the material + */ + void InitializeMaterial(); + + + ///@} + ///@name Protected Access + ///@{ + + ///@} + ///@name Protected Inquiry + ///@{ + + ///@} + ///@name Protected LifeCycle + ///@{ + +private: + ///@name Static Member Variables + ///@{ + + ///@} + ///@name Member Variables + ///@{ + + ///@} + ///@name Private Operators + ///@{ + + ///@} + ///@name Private Operations + ///@{ + + + ///@} + ///@name Private Access + ///@{ + + ///@} + ///@name Private Inquiry + ///@{ + + ///@} + ///@name Serialization + ///@{ + + friend class Serializer; + + void save(Serializer &rSerializer) const override; + + void load(Serializer &rSerializer) override; + +}; // class LinearTrussElement2D. + +///@} +///@name Type Definitions +///@{ + + +///@} +///@name Input and output +///@{ + +} // namespace Kratos. diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.cpp:Zone.Identifier b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.cpp:Zone.Identifier new file mode 100644 index 000000000000..e69de29bb2d1 diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier new file mode 100644 index 000000000000..e69de29bb2d1 From d15c64b59d6c8ef248d8b091e0994cb8502e7427 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Sun, 3 Nov 2024 16:46:02 +0100 Subject: [PATCH 03/32] remove bug --- .../timoshenko_beam_element_2D2N.cpp:Zone.Identifier | 0 .../truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier | 0 2 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.cpp:Zone.Identifier delete mode 100644 applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.cpp:Zone.Identifier b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.cpp:Zone.Identifier deleted file mode 100644 index e69de29bb2d1..000000000000 diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/timoshenko_beam_element_2D2N.h:Zone.Identifier deleted file mode 100644 index e69de29bb2d1..000000000000 From d3578f7dbe7431716742173af326f6f713a09e03 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Sun, 3 Nov 2024 17:05:13 +0100 Subject: [PATCH 04/32] more code --- .../linear_truss_element_2D.cpp | 250 +----------------- .../truss_elements/linear_truss_element_2D.h | 22 +- 2 files changed, 20 insertions(+), 252 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 859bb17a94e0..4ed6d2571b6c 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -82,8 +82,7 @@ Element::Pointer LinearTrussElement2D::Clone( { KRATOS_TRY - LinearTrussElement2D::Pointer p_new_elem = Kratos::make_intrusive> - (NewId, GetGeometry().Create(rThisNodes), pGetProperties()); + LinearTrussElement2D::Pointer p_new_elem = Kratos::make_intrusive>(NewId, GetGeometry().Create(rThisNodes), pGetProperties()); p_new_elem->SetData(this->GetData()); p_new_elem->Set(Flags(*this)); @@ -109,20 +108,16 @@ void LinearTrussElement2D::EquationIdVector( { const auto& r_geometry = this->GetGeometry(); const SizeType number_of_nodes = r_geometry.size(); - const SizeType dofs_per_node = GetDoFsPerNode(); // u, v, theta - IndexType local_index = 0; - if (rResult.size() != dofs_per_node * number_of_nodes) - rResult.resize(dofs_per_node * number_of_nodes, false); + if (rResult.size() != SystemSize) + rResult.resize(SystemSize, false); const IndexType xpos = this->GetGeometry()[0].GetDofPosition(DISPLACEMENT_X); - const IndexType rot_pos = this->GetGeometry()[0].GetDofPosition(ROTATION_Z); for (IndexType i = 0; i < number_of_nodes; ++i) { rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_X, xpos ).EquationId(); rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_Y, xpos + 1).EquationId(); - rResult[local_index++] = r_geometry[i].GetDof(ROTATION_Z , rot_pos ).EquationId(); } } @@ -139,14 +134,12 @@ void LinearTrussElement2D::GetDofList( const auto& r_geom = GetGeometry(); const SizeType number_of_nodes = r_geom.size(); - const SizeType dofs_per_node = GetDoFsPerNode(); // u, v, theta - rElementalDofList.resize(dofs_per_node * number_of_nodes); + rElementalDofList.resize(SystemSize); for (IndexType i = 0; i < number_of_nodes; ++i) { - const SizeType index = i * dofs_per_node; + const SizeType index = i * DofsPerNode; rElementalDofList[index] = r_geom[i].pGetDof(DISPLACEMENT_X); rElementalDofList[index + 1] = r_geom[i].pGetDof(DISPLACEMENT_Y); - rElementalDofList[index + 2] = r_geom[i].pGetDof(ROTATION_Z ); } KRATOS_CATCH("") } @@ -158,18 +151,10 @@ template void LinearTrussElement2D::GetShapeFunctionsValues( VectorType& rN, const double Length, - const double Phi, const double xi ) const { - if (rN.size() != 4) - rN.resize(4, false); - const double one_plus_phi = 1.0 + Phi; - const double xi_square = xi * xi; - rN[0] = (xi - 1.0) * (xi + xi_square - 2.0 * one_plus_phi) / (4.0 * one_plus_phi); - rN[1] = (1.0 - xi_square) * (1.0 - xi + Phi) * Length / (8.0 * one_plus_phi); - rN[2] = (1.0 + xi) * (xi - xi_square + 2.0 * one_plus_phi) / (4.0 * one_plus_phi); - rN[3] = (xi_square - 1.0) * (1.0 + xi + Phi) * Length / (8.0 * one_plus_phi); + } /***********************************************************************************/ @@ -179,63 +164,19 @@ template void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( VectorType& rN, const double Length, - const double Phi, const double xi ) const { - if (rN.size() != 4) - rN.resize(4, false); - const double one_plus_phi = 1.0 + Phi; - const double xi_square = xi * xi; - rN[0] = (-6.0 + 6.0 * xi_square - 4.0 * Phi) / (4.0 * one_plus_phi * Length); - rN[1] = (-1.0 + 3.0 * xi_square - 2.0 * xi * one_plus_phi) / (4.0 * one_plus_phi); - rN[2] = (6.0 - 6.0 * xi_square + 4.0 * Phi) / (4.0 * one_plus_phi * Length); - rN[3] = (-1.0 + 3.0 * xi_square + 2.0 * xi * one_plus_phi) / (4.0 * one_plus_phi); + } /***********************************************************************************/ /***********************************************************************************/ template -void LinearTrussElement2D::GetNodalValuesVector(VectorType& rNodalValues) const +void LinearTrussElement2D::GetNodalValuesVector(SystemSizeBoundedArrayType& rNodalValues) const { - // if (rNodalValues.size() != 6) - // rNodalValues.resize(6, false); - // const auto &r_geom = GetGeometry(); - - // const double angle = GetAngle(); - - // if (std::abs(angle) > std::numeric_limits::epsilon()) { - // BoundedMatrix T; - // BoundedVector global_values; - // BoundedMatrix global_size_T; - // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); - // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); - - // const auto &r_displ_0 = r_geom[0].FastGetSolutionStepValue(DISPLACEMENT); - // global_values[0] = r_displ_0[0]; - // global_values[1] = r_displ_0[1]; - // global_values[2] = r_geom[0].FastGetSolutionStepValue(ROTATION_Z); - - // const auto& r_displ_1 = r_geom[1].FastGetSolutionStepValue(DISPLACEMENT); - // global_values[3] = r_displ_1[0]; - // global_values[4] = r_displ_1[1]; - // global_values[5] = r_geom[1].FastGetSolutionStepValue(ROTATION_Z); - - // // We rotate to local axes - // noalias(rNodalValues) = prod(trans(global_size_T), global_values); - - // } else { - // const auto &r_displ_0 = r_geom[0].FastGetSolutionStepValue(DISPLACEMENT); - // rNodalValues[0] = r_displ_0[0]; - // rNodalValues[1] = r_displ_0[1]; - // rNodalValues[2] = r_geom[0].FastGetSolutionStepValue(ROTATION_Z); - - // const auto& r_displ_1 = r_geom[1].FastGetSolutionStepValue(DISPLACEMENT); - // rNodalValues[3] = r_displ_1[0]; - // rNodalValues[4] = r_displ_1[1]; - // rNodalValues[5] = r_geom[1].FastGetSolutionStepValue(ROTATION_Z); - // } + } /***********************************************************************************/ @@ -248,15 +189,7 @@ array_1d LinearTrussElement2D::GetLocalAxesBodyForce( const IndexType PointNumber ) const { - const double angle = GetAngle(); - const auto body_force = StructuralMechanicsElementUtilities::GetBodyForce(*this, rIntegrationPoints, PointNumber); - - const double c = std::cos(angle); - const double s = std::sin(angle); - array_1d local_body_force = ZeroVector(3); - local_body_force[0] = c * body_force[0] + s * body_force[1]; - local_body_force[1] = -s * body_force[0] + c * body_force[1]; - return local_body_force; + } /***********************************************************************************/ @@ -269,130 +202,7 @@ void LinearTrussElement2D::CalculateLocalSystem( const ProcessInfo& rProcessInfo ) { - // KRATOS_TRY; - // const auto &r_props = GetProperties(); - // const auto &r_geometry = GetGeometry(); - // const SizeType number_of_nodes = r_geometry.size(); - // const SizeType mat_size = GetDoFsPerNode() * number_of_nodes; - - // if (rLHS.size1() != mat_size || rLHS.size2() != mat_size) { - // rLHS.resize(mat_size, mat_size, false); - // } - // noalias(rLHS) = ZeroMatrix(mat_size, mat_size); - - // if (rRHS.size() != mat_size) { - // rRHS.resize(mat_size, false); - // } - // noalias(rRHS) = ZeroVector(mat_size); - - // const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); - - // 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, true); - - // const double length = CalculateLength(); - // const double Phi = StructuralMechanicsElementUtilities::CalculatePhi(r_props, length); - // const double J = 0.5 * length; - // const double area = r_props[CROSS_AREA]; - - // // Let's initialize the cl values - // VectorType strain_vector(3), stress_vector(3); - // MatrixType constitutive_matrix(3, 3); - // strain_vector.clear(); - // cl_values.SetStrainVector(strain_vector); - // cl_values.SetStressVector(stress_vector); - // cl_values.SetConstitutiveMatrix(constitutive_matrix); - // VectorType nodal_values(mat_size); - // GetNodalValuesVector(nodal_values); - // VectorType global_size_N_2(mat_size), global_size_N(mat_size), N_u_derivatives(number_of_nodes), - // N_theta_derivatives(mat_size-number_of_nodes), N_theta(mat_size-number_of_nodes), N_derivatives(mat_size-number_of_nodes), - // N_u(number_of_nodes), N_shape(mat_size-number_of_nodes), N_s(mat_size-number_of_nodes); - - // // Loop over the integration points - // for (SizeType IP = 0; IP < integration_points.size(); ++IP) { - // const auto local_body_forces = GetLocalAxesBodyForce(*this, integration_points, IP); - - // global_size_N.clear(); - // const double xi = integration_points[IP].X(); - // const double weight = integration_points[IP].Weight(); - // const double jacobian_weight = weight * J; - - // CalculateGeneralizedStrainsVector(strain_vector, length, Phi, xi, nodal_values); - - // mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); - // const Vector &r_generalized_stresses = cl_values.GetStressVector(); - // const double N = r_generalized_stresses[0]; - // const double M = r_generalized_stresses[1]; - // const double V = r_generalized_stresses[2]; - - // const MatrixType& r_constitutive_matrix = cl_values.GetConstitutiveMatrix(); - // const double dN_dEl = r_constitutive_matrix(0, 0); - // const double dM_dkappa = r_constitutive_matrix(1, 1); - // const double dV_dgamma = r_constitutive_matrix(2, 2); - - // GetFirstDerivativesNu0ShapeFunctionsValues(N_u_derivatives, length, Phi, xi); - // GetFirstDerivativesNThetaShapeFunctionsValues(N_theta_derivatives, length, Phi, xi); - // GetNThetaShapeFunctionsValues(N_theta, length, Phi, xi); - // GetFirstDerivativesShapeFunctionsValues(N_derivatives, length, Phi, xi); - // GetShapeFunctionsValues(N_shape, length, Phi, xi); - // GetNu0ShapeFunctionsValues(N_u, length, Phi, xi); - // noalias(N_s) = N_derivatives - N_theta; - - // // Axial contributions - // GlobalSizeAxialVector(global_size_N, N_u_derivatives); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N) * dN_dEl * jacobian_weight; - // noalias(rRHS) -= global_size_N * N * jacobian_weight; - - // // In here we add the cross terms - // GlobalSizeVector(global_size_N_2, N_theta_derivatives); - // const double dN_dkappa = r_constitutive_matrix(0, 1); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dN_dkappa * jacobian_weight; - - // GlobalSizeVector(global_size_N_2, N_s); - // const double dN_dgamma = r_constitutive_matrix(0, 2); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dN_dgamma * jacobian_weight; - - // // Bending contributions - // GlobalSizeVector(global_size_N, N_theta_derivatives); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N) * dM_dkappa * jacobian_weight; - // noalias(rRHS) -= global_size_N * M * jacobian_weight; - - // // In here we add the cross terms - // GlobalSizeAxialVector(global_size_N_2, N_u_derivatives); - // const double dM_dEl = r_constitutive_matrix(1, 0); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dM_dEl * jacobian_weight; - - // GlobalSizeVector(global_size_N_2, N_s); - // const double dM_dgamma = r_constitutive_matrix(1, 2); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dM_dgamma * jacobian_weight; - - // // Shear contributions - // GlobalSizeVector(global_size_N, N_s); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N) * dV_dgamma * jacobian_weight; - // noalias(rRHS) -= global_size_N * V * jacobian_weight; - - // // In here we add the cross terms - // GlobalSizeAxialVector(global_size_N_2, N_u_derivatives); - // const double dV_dEl = r_constitutive_matrix(2, 0); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dV_dEl * jacobian_weight; - - // GlobalSizeVector(global_size_N_2, N_theta_derivatives); - // const double dV_dkappa = r_constitutive_matrix(2, 1); - // noalias(rLHS) += outer_prod(global_size_N, global_size_N_2) * dV_dkappa * jacobian_weight; - - - // // Now we add the body forces contributions - // GlobalSizeAxialVector(global_size_N, N_u); - // noalias(rRHS) += global_size_N * local_body_forces[0] * jacobian_weight * area; - // GlobalSizeVector(global_size_N, N_shape); - // noalias(rRHS) += global_size_N * local_body_forces[1] * jacobian_weight * area; - - // } - // RotateAll(rLHS, rRHS, r_geometry); - - // KRATOS_CATCH(""); + } /***********************************************************************************/ @@ -428,16 +238,7 @@ void LinearTrussElement2D::RotateLHS( const GeometryType& rGeometry ) { - // const double angle = GetAngle(); - - // if (std::abs(angle) > std::numeric_limits::epsilon()) { - // BoundedMatrix T, Tt; - // BoundedMatrix global_size_T, aux_product; - // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); - // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); - // noalias(aux_product) = prod(rLHS, trans(global_size_T)); - // noalias(rLHS) = prod(global_size_T, aux_product); - // } + } /***********************************************************************************/ @@ -449,17 +250,7 @@ void LinearTrussElement2D::RotateRHS( const GeometryType& rGeometry ) { - // const double angle = GetAngle(); - // if (std::abs(angle) > std::numeric_limits::epsilon()) { - // BoundedMatrix T; - // BoundedMatrix global_size_T; - // BoundedVector local_rhs; - // noalias(local_rhs) = rRHS; - // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); - // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); - - // noalias(rRHS) = prod(global_size_T, local_rhs); - // } + } /***********************************************************************************/ @@ -472,20 +263,7 @@ void LinearTrussElement2D::RotateAll( const GeometryType& rGeometry ) { - // const double angle = GetAngle(); - // if (std::abs(angle) > std::numeric_limits::epsilon()) { - // BoundedMatrix T; - // BoundedMatrix global_size_T, aux_product; - // BoundedVector local_rhs; - // StructuralMechanicsElementUtilities::BuildRotationMatrixForBeam(T, angle); - // StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NBeam(T, global_size_T); - - // noalias(local_rhs) = rRHS; - // noalias(rRHS) = prod(global_size_T, local_rhs); - - // noalias(aux_product) = prod(rLHS, trans(global_size_T)); - // noalias(rLHS) = prod(global_size_T, aux_product); - // } + } /***********************************************************************************/ diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h index 4f6ae346957c..b3d54e4c6a7f 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -47,8 +47,7 @@ using SizeType = std::size_t; /** * @class LinearTrussElement2D * @ingroup StructuralMechanicsApplication - * @brief This is the Linear TRUSS element of 2 nodes. - * 2 and 3 noded elements + * @brief This is the Linear TRUSS element of 2 and 3 nodes. * @author Alejandro Cornejo */ template @@ -129,28 +128,19 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D ///@name Operations ///@{ - /** - * @brief Indicates the amount of DoFs per node (u0, v, theta) - */ - IndexType GetDoFsPerNode() const - { - return 2; - } - /** * @brief This method returns the angle of the FE axis */ double GetAngle() const { - // return StructuralMechanicsElementUtilities::GetReferenceRotationAngle2D2NBeam(GetGeometry()); - return 1.0; + return StructuralMechanicsElementUtilities::GetReferenceRotationAngle2D2NBeam(GetGeometry()); } /** - * @brief Returns a 6 component vector including the values of the DoFs + * @brief Returns a n component vector including the values of the DoFs * in LOCAL beam axes */ - void GetNodalValuesVector(VectorType& rNodalValue) const; + void GetNodalValuesVector(SystemSizeBoundedArrayType& rNodalValue) const; /** * @brief Computes the length of the FE and returns it @@ -241,8 +231,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param Phi The shear slenderness parameter * @param xi The coordinate in the natural axes */ - void GetShapeFunctionsValues(VectorType& rN, const double Length, const double Phi, const double xi) const; - void GetFirstDerivativesShapeFunctionsValues(VectorType& rN, const double Length, const double Phi, const double xi) const; + void GetShapeFunctionsValues(VectorType& rN, const double Length, const double xi) const; + void GetFirstDerivativesShapeFunctionsValues(VectorType& rN, const double Length, const double xi) const; /** * @brief This function rotates the LHS from local to global coordinates From f68035fd0311c5ded947c49ea3958fa9828743a4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 11:34:58 +0100 Subject: [PATCH 05/32] adding --- .../truss_elements/linear_truss_element_2D.cpp | 17 ++++++++++++----- .../truss_elements/linear_truss_element_2D.h | 16 +++++++++------- 2 files changed, 21 insertions(+), 12 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 4ed6d2571b6c..3f636dc227d5 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -24,6 +24,9 @@ namespace Kratos { +/***********************************************************************************/ +/***********************************************************************************/ + template void LinearTrussElement2D::Initialize(const ProcessInfo& rCurrentProcessInfo) { @@ -35,7 +38,11 @@ void LinearTrussElement2D::Initialize(const ProcessInfo& rCurrentProces if (GetProperties().Has(INTEGRATION_ORDER) ) { mThisIntegrationMethod = static_cast(GetProperties()[INTEGRATION_ORDER] - 1); } else { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_3; + if constexpr (NNodes == 2) { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; + } else { + mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; + } } } @@ -56,13 +63,13 @@ template void LinearTrussElement2D::InitializeMaterial() { KRATOS_TRY + const auto &r_props = GetProperties(); - if (GetProperties()[CONSTITUTIVE_LAW] != nullptr) { + if (r_props[CONSTITUTIVE_LAW] != nullptr) { const auto& r_geometry = GetGeometry(); - const auto& r_properties = GetProperties(); auto N_values = Vector(); for (IndexType point_number = 0; point_number < mConstitutiveLawVector.size(); ++point_number) { - mConstitutiveLawVector[point_number] = GetProperties()[CONSTITUTIVE_LAW]->Clone(); + mConstitutiveLawVector[point_number] = r_props[CONSTITUTIVE_LAW]->Clone(); mConstitutiveLawVector[point_number]->InitializeMaterial(r_properties, r_geometry, N_values); } } else @@ -113,7 +120,7 @@ void LinearTrussElement2D::EquationIdVector( if (rResult.size() != SystemSize) rResult.resize(SystemSize, false); - const IndexType xpos = this->GetGeometry()[0].GetDofPosition(DISPLACEMENT_X); + const IndexType xpos = this->GetGeometry()[0].GetDofPosition(DISPLACEMENT_X); for (IndexType i = 0; i < number_of_nodes; ++i) { rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_X, xpos ).EquationId(); diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h index b3d54e4c6a7f..e189092557a9 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -62,8 +62,9 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D /// The base element type using BaseType = Element; + static constexpr SizeType NNodes = TNNodes; static constexpr SizeType DofsPerNode = 2; - static constexpr SizeType SystemSize = DofsPerNode * TNNodes; + static constexpr SizeType SystemSize = DofsPerNode * NNodes; using SystemSizeBoundedArrayType = array_1d; // Counted pointer of BaseSolidElement @@ -81,7 +82,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D // Constructor using an array of nodes LinearTrussElement2D(IndexType NewId, GeometryType::Pointer pGeometry) : Element(NewId, pGeometry) { - if constexpr (TNNodes == 2) { + if constexpr (NNodes == 2) { mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; } else { mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; @@ -93,7 +94,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D : Element(NewId,pGeometry,pProperties) { // This is needed to prevent uninitialised integration method in inactive elements - if constexpr (TNNodes == 2) { + if constexpr (NNodes == 2) { mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; } else { mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; @@ -147,6 +148,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D */ double CalculateLength() const { + // Same implementation for 2N and 3N return StructuralMechanicsElementUtilities::CalculateReferenceLength2D2N(*this); } @@ -210,7 +212,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D /** * @brief Returns the set of integration points */ - const GeometryType::IntegrationPointsArrayType IntegrationPoints() const + const GeometryType::IntegrationPointsArrayType IntegrationPoints() const { return GetGeometry().IntegrationPoints(); } @@ -239,7 +241,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param rLHS the left hand side * @param rGeometry the geometry of the FE */ - virtual void RotateLHS( + void RotateLHS( MatrixType &rLHS, const GeometryType &rGeometry); @@ -248,7 +250,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param rRHS the right hand side * @param rGeometry the geometry of the FE */ - virtual void RotateRHS( + void RotateRHS( VectorType &rRHS, const GeometryType &rGeometry); @@ -258,7 +260,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param rRHS the right hand side * @param rGeometry the geometry of the FE */ - virtual void RotateAll( + void RotateAll( MatrixType &rLHS, VectorType &rRHS, const GeometryType &rGeometry); From b803cd7241fe0cb9fd6d9afbfc839e3e1f33bab0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 11:47:22 +0100 Subject: [PATCH 06/32] shape funct --- .../truss_elements/linear_truss_element_2D.cpp | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 3f636dc227d5..4d5d621084d8 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -161,7 +161,16 @@ void LinearTrussElement2D::GetShapeFunctionsValues( const double xi ) const { - + if (rN.size() != SystemSize) + rN.resize(SystemSize, false); + if constexpr (NNodes == 2) { + rN[0] = 0.5 * (1.0 - xi); + rN[2] = 0.5 * (1.0 + xi); + } else { // 3N + rN[0] = 0.5 * xi * (xi - 1.0); + rN[2] = (1.0 - std::pow(xi, 2)); + rN[4] = 0.5 * xi * (xi + 1.0); + } } /***********************************************************************************/ From 264e19a2298b7d5a08da75d6e6b8e50277e5b88b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 11:51:36 +0100 Subject: [PATCH 07/32] derivatives of shape functions --- .../truss_elements/linear_truss_element_2D.cpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 4d5d621084d8..1a27a696d0f7 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -178,12 +178,23 @@ void LinearTrussElement2D::GetShapeFunctionsValues( template void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( - VectorType& rN, + VectorType& rdN_dX, const double Length, const double xi ) const { - + if (rdN_dX.size() != SystemSize) + rdN_dX.resize(SystemSize, false); + if constexpr (NNodes == 2) { + const double inverse_l = 1.0 / Length; + rdN_dX[0] = -inverse_l; + rdN_dX[2] = inverse_l; + } else { // 3N + rdN_dX[0] = xi - 0.5; + rdN_dX[2] = -2.0 * xi; + rdN_dX[4] = xi + 0.5; + rdN_dX *= 2.0 / Length; + } } /***********************************************************************************/ From ae5e3ecfd45db0f80909103f679b02b4b6fe9e54 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 14:34:11 +0100 Subject: [PATCH 08/32] add utils for truss rotations --- ...structural_mechanics_element_utilities.cpp | 65 +++++++++++++++++++ .../structural_mechanics_element_utilities.h | 22 +++++++ 2 files changed, 87 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp index c4ba482536cb..68a30f72edb2 100644 --- a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp +++ b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp @@ -364,6 +364,71 @@ void BuildRotationMatrixForBeam( /***********************************************************************************/ /***********************************************************************************/ +void BuildRotationMatrixForTruss( + BoundedMatrix& rRotationMatrix, + const double AlphaAngle +) +{ + rRotationMatrix.clear(); + const double s = std::sin(AlphaAngle); + const double c = std::cos(AlphaAngle); + rRotationMatrix(0, 0) = c; + rRotationMatrix(0, 1) = -s; + rRotationMatrix(1, 0) = s; + rRotationMatrix(1, 1) = c; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void BuildElementSizeRotationMatrixFor2D2NTruss( + const BoundedMatrix& rRotationMatrix, + BoundedMatrix& rElementSizeRotationMatrix +) +{ + rElementSizeRotationMatrix.clear(); + + rElementSizeRotationMatrix(0, 0) = rRotationMatrix(0, 0); + rElementSizeRotationMatrix(0, 1) = rRotationMatrix(0, 1); + rElementSizeRotationMatrix(1, 0) = rRotationMatrix(1, 0); + rElementSizeRotationMatrix(1, 1) = rRotationMatrix(1, 1); + + + rElementSizeRotationMatrix(2, 2) = rRotationMatrix(0, 0); + rElementSizeRotationMatrix(2, 3) = rRotationMatrix(0, 1); + rElementSizeRotationMatrix(3, 2) = rRotationMatrix(1, 0); + rElementSizeRotationMatrix(3, 3) = rRotationMatrix(1, 1); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +void BuildElementSizeRotationMatrixFor2D3NTruss( + const BoundedMatrix& rRotationMatrix, + BoundedMatrix& rElementSizeRotationMatrix +) +{ + rElementSizeRotationMatrix.clear(); + + rElementSizeRotationMatrix(0, 0) = rRotationMatrix(0, 0); + rElementSizeRotationMatrix(0, 1) = rRotationMatrix(0, 1); + rElementSizeRotationMatrix(1, 0) = rRotationMatrix(1, 0); + rElementSizeRotationMatrix(1, 1) = rRotationMatrix(1, 1); + + rElementSizeRotationMatrix(2, 2) = rRotationMatrix(0, 0); + rElementSizeRotationMatrix(2, 3) = rRotationMatrix(0, 1); + rElementSizeRotationMatrix(3, 2) = rRotationMatrix(1, 0); + rElementSizeRotationMatrix(3, 3) = rRotationMatrix(1, 1); + + rElementSizeRotationMatrix(4, 4) = rRotationMatrix(0, 0); + rElementSizeRotationMatrix(4, 5) = rRotationMatrix(0, 1); + rElementSizeRotationMatrix(5, 4) = rRotationMatrix(1, 0); + rElementSizeRotationMatrix(5, 5) = rRotationMatrix(1, 1); +} + +/***********************************************************************************/ +/***********************************************************************************/ + double GetReferenceRotationAngle2D2NBeam(const GeometryType& rGeometry) { const auto &r_node_1 = rGeometry[0]; diff --git a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h index 89fbd9a6affc..349d69682d2b 100644 --- a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h +++ b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h @@ -286,6 +286,28 @@ void BuildElementSizeRotationMatrixFor2D3NBeam( const BoundedMatrix& rRotationMatrix, BoundedMatrix& rElementSizeRotationMatrix); +void BuildRotationMatrixForTruss( + BoundedMatrix& rRotationMatrix, + const double AlphaAngle); + +/** + * @brief This function fills an element size rotation matrix a local rotation matrix + * @param rRotationMatrix The rotation matrix from local to global axes + * It assumes 3 dofs per node: u,v,theta + */ +void BuildElementSizeRotationMatrixFor2D2NTruss( + const BoundedMatrix& rRotationMatrix, + BoundedMatrix& rElementSizeRotationMatrix); + +/** + * @brief This function fills an element size rotation matrix a local rotation matrix + * @param rRotationMatrix The rotation matrix from local to global axes + * It assumes 3 dofs per node: u,v,theta + */ +void BuildElementSizeRotationMatrixFor2D3NTruss( + const BoundedMatrix& rRotationMatrix, + BoundedMatrix& rElementSizeRotationMatrix); + /** * @brief This function computes the inclination angle of a 2 noded beam * @param rGeometry The geometry of the beam From 68a66938e2e5a292cefad6babd553729ae7e3493 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 14:47:54 +0100 Subject: [PATCH 09/32] more code --- .../linear_truss_element_2D.cpp | 42 +++++++++++++++++-- 1 file changed, 39 insertions(+), 3 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 1a27a696d0f7..ba3136fdcc02 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -70,7 +70,7 @@ void LinearTrussElement2D::InitializeMaterial() auto N_values = Vector(); for (IndexType point_number = 0; point_number < mConstitutiveLawVector.size(); ++point_number) { mConstitutiveLawVector[point_number] = r_props[CONSTITUTIVE_LAW]->Clone(); - mConstitutiveLawVector[point_number]->InitializeMaterial(r_properties, r_geometry, N_values); + mConstitutiveLawVector[point_number]->InitializeMaterial(r_props, r_geometry, N_values); } } else KRATOS_ERROR << "A constitutive law needs to be specified for the element with ID " << this->Id() << std::endl; @@ -193,7 +193,7 @@ void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( rdN_dX[0] = xi - 0.5; rdN_dX[2] = -2.0 * xi; rdN_dX[4] = xi + 0.5; - rdN_dX *= 2.0 / Length; + rdN_dX *= 2.0 / Length; // The Jacobian } } @@ -203,7 +203,35 @@ void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( template void LinearTrussElement2D::GetNodalValuesVector(SystemSizeBoundedArrayType& rNodalValues) const { + if (rNodalValues.size() != SystemSize) + rNodalValues.resize(SystemSize, false); + const auto &r_geom = GetGeometry(); + const double angle = GetAngle(); + + BoundedVector global_values; + + // We fill the vector with global values + for (SizeType i = 0; i < NNodes; ++i) { + const auto& r_displ = r_geom[i].FastGetSolutionStepValue(DISPLACEMENT); + global_values[i * DofsPerNode] = r_displ[0]; + global_values[i * DofsPerNode + 1] = r_displ[1]; + } + + if (std::abs(angle) > std::numeric_limits::epsilon()) { + BoundedMatrix T; + BoundedMatrix global_size_T; + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); + if constexpr (NNodes == 2) { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + } else { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D3NTruss(T, global_size_T); + } + // We rotate to local axes + noalias(rNodalValues) = prod(trans(global_size_T), global_values); + } else { + noalias(rNodalValues) = global_values; + } } /***********************************************************************************/ @@ -216,7 +244,15 @@ array_1d LinearTrussElement2D::GetLocalAxesBodyForce( const IndexType PointNumber ) const { - + const double angle = GetAngle(); + const auto body_force = StructuralMechanicsElementUtilities::GetBodyForce(*this, rIntegrationPoints, PointNumber); + + const double c = std::cos(angle); + const double s = std::sin(angle); + array_1d local_body_force = ZeroVector(3); + local_body_force[0] = c * body_force[0] + s * body_force[1]; + local_body_force[1] = -s * body_force[0] + c * body_force[1]; + return local_body_force; } /***********************************************************************************/ From fade7603a36fef9186e3d8ee0cfd164a393bc167 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 15:44:59 +0100 Subject: [PATCH 10/32] rotations --- .../linear_truss_element_2D.cpp | 140 +++++++++++++++++- .../truss_elements/linear_truss_element_2D.h | 5 +- 2 files changed, 141 insertions(+), 4 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index ba3136fdcc02..237d4596b7f4 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -156,13 +156,15 @@ void LinearTrussElement2D::GetDofList( template void LinearTrussElement2D::GetShapeFunctionsValues( - VectorType& rN, + SystemSizeBoundedArrayType& rN, const double Length, const double xi ) const { if (rN.size() != SystemSize) rN.resize(SystemSize, false); + + rN.clear(); if constexpr (NNodes == 2) { rN[0] = 0.5 * (1.0 - xi); rN[2] = 0.5 * (1.0 + xi); @@ -176,15 +178,41 @@ void LinearTrussElement2D::GetShapeFunctionsValues( /***********************************************************************************/ /***********************************************************************************/ +template +void LinearTrussElement2D::GetShapeFunctionsValuesY( + SystemSizeBoundedArrayType& rN, + const double Length, + const double xi + ) const +{ + if (rN.size() != SystemSize) + rN.resize(SystemSize, false); + + rN.clear(); + if constexpr (NNodes == 2) { + rN[1] = 0.5 * (1.0 - xi); + rN[3] = 0.5 * (1.0 + xi); + } else { // 3N + rN[1] = 0.5 * xi * (xi - 1.0); + rN[3] = (1.0 - std::pow(xi, 2)); + rN[5] = 0.5 * xi * (xi + 1.0); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + template void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( - VectorType& rdN_dX, + SystemSizeBoundedArrayType& rdN_dX, const double Length, const double xi ) const { if (rdN_dX.size() != SystemSize) rdN_dX.resize(SystemSize, false); + + rdN_dX.clear(); if constexpr (NNodes == 2) { const double inverse_l = 1.0 / Length; rdN_dX[0] = -inverse_l; @@ -265,6 +293,68 @@ void LinearTrussElement2D::CalculateLocalSystem( const ProcessInfo& rProcessInfo ) { + KRATOS_TRY; + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + + if (rLHS.size1() != SystemSize || rLHS.size2() != SystemSize) { + rLHS.resize(SystemSize, SystemSize, false); + } + noalias(rLHS) = ZeroMatrix(SystemSize, SystemSize); + + if (rRHS.size() != SystemSize) { + rRHS.resize(SystemSize, false); + } + noalias(rRHS) = ZeroVector(SystemSize); + + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + + 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, true); + + const double length = CalculateLength(); + const double J = 0.5 * length; + const double area = r_props[CROSS_AREA]; + + // Let's initialize the cl values + VectorType strain_vector(1), stress_vector(1); + MatrixType constitutive_matrix(1, 1); // Young modulus + + strain_vector.clear(); + cl_values.SetStrainVector(strain_vector); + cl_values.SetStressVector(stress_vector); + cl_values.SetConstitutiveMatrix(constitutive_matrix); + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); + + SystemSizeBoundedArrayType B, N_shape, N_shapeY; + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const auto local_body_forces = GetLocalAxesBodyForce(*this, integration_points, IP); + + GetShapeFunctionsValues(N_shape, length, xi); + GetShapeFunctionsValuesY(N_shapeY, length, xi); + GetFirstDerivativesShapeFunctionsValues(B, length, xi); + + const double xi = integration_points[IP].X(); + const double weight = integration_points[IP].Weight(); + const double jacobian_weight = weight * J * area; + + strain_vector[0] = inner_prod(B, nodal_values); + mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix + + noalias(rLHS) += outer_prod(B, B) * constitutive_matrix(0, 0) * jacobian_weight; + noalias(rRHS) -= B * stress_vector[0] * jacobian_weight; + + noalias(rRHS) += N_shape * local_body_forces[0] * jacobian_weight; + noalias(rRHS) += N_shapeY * local_body_forces[1] * jacobian_weight; + } + RotateAll(rLHS, rRHS, r_geometry); + + KRATOS_CATCH("") } @@ -301,7 +391,21 @@ void LinearTrussElement2D::RotateLHS( const GeometryType& rGeometry ) { + const double angle = GetAngle(); + + if (std::abs(angle) > std::numeric_limits::epsilon()) { + BoundedMatrix T, Tt; + BoundedMatrix global_size_T, aux_product; + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); + if constexpr (NNodes == 2) { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + } else { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D3NTruss(T, global_size_T); + } + noalias(aux_product) = prod(rLHS, trans(global_size_T)); + noalias(rLHS) = prod(global_size_T, aux_product); + } } /***********************************************************************************/ @@ -313,7 +417,21 @@ void LinearTrussElement2D::RotateRHS( const GeometryType& rGeometry ) { + const double angle = GetAngle(); + if (std::abs(angle) > std::numeric_limits::epsilon()) { + BoundedMatrix T; + BoundedMatrix global_size_T; + BoundedVector local_rhs; + noalias(local_rhs) = rRHS; + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, angle); + if constexpr (NNodes == 2) { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + } else { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D3NTruss(T, global_size_T); + } + noalias(rRHS) = prod(global_size_T, local_rhs); + } } /***********************************************************************************/ @@ -326,7 +444,25 @@ void LinearTrussElement2D::RotateAll( const GeometryType& rGeometry ) { + const double angle = GetAngle(); + if (std::abs(angle) > std::numeric_limits::epsilon()) { + BoundedMatrix T; + BoundedMatrix global_size_T, aux_product; + BoundedVector local_rhs; + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, angle); + if constexpr (NNodes == 2) { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + } else { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D3NTruss(T, global_size_T); + } + + noalias(local_rhs) = rRHS; + noalias(rRHS) = prod(global_size_T, local_rhs); + + noalias(aux_product) = prod(rLHS, trans(global_size_T)); + noalias(rLHS) = prod(global_size_T, aux_product); + } } /***********************************************************************************/ diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h index e189092557a9..dae49d49a3a3 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -233,8 +233,9 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param Phi The shear slenderness parameter * @param xi The coordinate in the natural axes */ - void GetShapeFunctionsValues(VectorType& rN, const double Length, const double xi) const; - void GetFirstDerivativesShapeFunctionsValues(VectorType& rN, const double Length, const double xi) const; + void GetShapeFunctionsValues(SystemSizeBoundedArrayType& rN, const double Length, const double xi) const; + void GetShapeFunctionsValuesY(SystemSizeBoundedArrayType& rN, const double Length, const double xi) const; + void GetFirstDerivativesShapeFunctionsValues(SystemSizeBoundedArrayType& rN, const double Length, const double xi) const; /** * @brief This function rotates the LHS from local to global coordinates From 3773baa60a29ab038e396642981fd4856754203f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 4 Nov 2024 15:57:19 +0100 Subject: [PATCH 11/32] calculate LHS --- .../linear_truss_element_2D.cpp | 23 ++++++++----------- .../truss_elements/linear_truss_element_2D.h | 9 +++----- 2 files changed, 13 insertions(+), 19 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 237d4596b7f4..188259f01907 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -327,7 +327,7 @@ void LinearTrussElement2D::CalculateLocalSystem( cl_values.SetStressVector(stress_vector); cl_values.SetConstitutiveMatrix(constitutive_matrix); SystemSizeBoundedArrayType nodal_values(SystemSize); - GetNodalValuesVector(nodal_values); + GetNodalValuesVector(nodal_values); // In local axes SystemSizeBoundedArrayType B, N_shape, N_shapeY; @@ -335,13 +335,13 @@ void LinearTrussElement2D::CalculateLocalSystem( for (SizeType IP = 0; IP < integration_points.size(); ++IP) { const auto local_body_forces = GetLocalAxesBodyForce(*this, integration_points, IP); + const double xi = integration_points[IP].X(); + const double weight = integration_points[IP].Weight(); + const double jacobian_weight = weight * J * area; GetShapeFunctionsValues(N_shape, length, xi); GetShapeFunctionsValuesY(N_shapeY, length, xi); GetFirstDerivativesShapeFunctionsValues(B, length, xi); - const double xi = integration_points[IP].X(); - const double weight = integration_points[IP].Weight(); - const double jacobian_weight = weight * J * area; strain_vector[0] = inner_prod(B, nodal_values); mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix @@ -352,7 +352,7 @@ void LinearTrussElement2D::CalculateLocalSystem( noalias(rRHS) += N_shape * local_body_forces[0] * jacobian_weight; noalias(rRHS) += N_shapeY * local_body_forces[1] * jacobian_weight; } - RotateAll(rLHS, rRHS, r_geometry); + RotateAll(rLHS, rRHS); // rotate to global KRATOS_CATCH("") @@ -387,8 +387,7 @@ void LinearTrussElement2D::CalculateRightHandSide( template void LinearTrussElement2D::RotateLHS( - MatrixType& rLHS, - const GeometryType& rGeometry + MatrixType& rLHS ) { const double angle = GetAngle(); @@ -413,8 +412,7 @@ void LinearTrussElement2D::RotateLHS( template void LinearTrussElement2D::RotateRHS( - VectorType& rRHS, - const GeometryType& rGeometry + VectorType& rRHS ) { const double angle = GetAngle(); @@ -423,7 +421,7 @@ void LinearTrussElement2D::RotateRHS( BoundedMatrix global_size_T; BoundedVector local_rhs; noalias(local_rhs) = rRHS; - StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, angle); + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); if constexpr (NNodes == 2) { StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); } else { @@ -440,8 +438,7 @@ void LinearTrussElement2D::RotateRHS( template void LinearTrussElement2D::RotateAll( MatrixType& rLHS, - VectorType& rRHS, - const GeometryType& rGeometry + VectorType& rRHS ) { const double angle = GetAngle(); @@ -449,7 +446,7 @@ void LinearTrussElement2D::RotateAll( BoundedMatrix T; BoundedMatrix global_size_T, aux_product; BoundedVector local_rhs; - StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, angle); + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); if constexpr (NNodes == 2) { StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h index dae49d49a3a3..e070dd6ae2f1 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -243,8 +243,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param rGeometry the geometry of the FE */ void RotateLHS( - MatrixType &rLHS, - const GeometryType &rGeometry); + MatrixType &rLHS); /** * @brief This function rotates the RHS from local to global coordinates @@ -252,8 +251,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D * @param rGeometry the geometry of the FE */ void RotateRHS( - VectorType &rRHS, - const GeometryType &rGeometry); + VectorType &rRHS); /** * @brief This function rotates the LHS and RHS from local to global coordinates @@ -263,8 +261,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D */ void RotateAll( MatrixType &rLHS, - VectorType &rRHS, - const GeometryType &rGeometry); + VectorType &rRHS); /** * @brief This function retrieves the body forces in local axes From 4d63e9e2b561c4fe3de7d8e85b5a94b8596e49eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Tue, 5 Nov 2024 08:18:14 +0100 Subject: [PATCH 12/32] calc LHS --- .../linear_truss_element_2D.cpp | 49 +++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 188259f01907..6118a0bad487 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -367,7 +367,56 @@ void LinearTrussElement2D::CalculateLeftHandSide( const ProcessInfo& rProcessInfo ) { + KRATOS_TRY; + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + + if (rLHS.size1() != SystemSize || rLHS.size2() != SystemSize) { + rLHS.resize(SystemSize, SystemSize, false); + } + noalias(rLHS) = ZeroMatrix(SystemSize, SystemSize); + + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + + 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, true); + + const double length = CalculateLength(); + const double J = 0.5 * length; + const double area = r_props[CROSS_AREA]; + + // Let's initialize the cl values + VectorType strain_vector(1), stress_vector(1); + MatrixType constitutive_matrix(1, 1); // Young modulus + + strain_vector.clear(); + cl_values.SetStrainVector(strain_vector); + cl_values.SetStressVector(stress_vector); + cl_values.SetConstitutiveMatrix(constitutive_matrix); + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); // In local axes + SystemSizeBoundedArrayType B; + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const auto local_body_forces = GetLocalAxesBodyForce(*this, integration_points, IP); + + const double xi = integration_points[IP].X(); + const double weight = integration_points[IP].Weight(); + const double jacobian_weight = weight * J * area; + GetFirstDerivativesShapeFunctionsValues(B, length, xi); + + strain_vector[0] = inner_prod(B, nodal_values); + mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix + + noalias(rLHS) += outer_prod(B, B) * constitutive_matrix(0, 0) * jacobian_weight; + } + RotateLHS(rLHS); // rotate to global + + KRATOS_CATCH("") } /***********************************************************************************/ From c09682fd3b8005c4a7f98d0a26a1da063a03f970 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Tue, 5 Nov 2024 08:30:17 +0100 Subject: [PATCH 13/32] RHS --- .../linear_truss_element_2D.cpp | 56 ++++++++++++++++++- 1 file changed, 55 insertions(+), 1 deletion(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 6118a0bad487..bc38d25248a3 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -355,7 +355,6 @@ void LinearTrussElement2D::CalculateLocalSystem( RotateAll(rLHS, rRHS); // rotate to global KRATOS_CATCH("") - } /***********************************************************************************/ @@ -428,7 +427,62 @@ void LinearTrussElement2D::CalculateRightHandSide( const ProcessInfo& rProcessInfo ) { + KRATOS_TRY; + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + + if (rRHS.size() != SystemSize) { + rRHS.resize(SystemSize, false); + } + noalias(rRHS) = ZeroVector(SystemSize); + + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + + 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, true); + + const double length = CalculateLength(); + const double J = 0.5 * length; + const double area = r_props[CROSS_AREA]; + + // Let's initialize the cl values + VectorType strain_vector(1), stress_vector(1); + MatrixType constitutive_matrix(1, 1); // Young modulus + strain_vector.clear(); + cl_values.SetStrainVector(strain_vector); + cl_values.SetStressVector(stress_vector); + cl_values.SetConstitutiveMatrix(constitutive_matrix); + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); // In local axes + + SystemSizeBoundedArrayType B, N_shape, N_shapeY; + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const auto local_body_forces = GetLocalAxesBodyForce(*this, integration_points, IP); + + const double xi = integration_points[IP].X(); + const double weight = integration_points[IP].Weight(); + const double jacobian_weight = weight * J * area; + GetShapeFunctionsValues(N_shape, length, xi); + GetShapeFunctionsValuesY(N_shapeY, length, xi); + GetFirstDerivativesShapeFunctionsValues(B, length, xi); + + + strain_vector[0] = inner_prod(B, nodal_values); + mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix + + noalias(rRHS) -= B * stress_vector[0] * jacobian_weight; + + noalias(rRHS) += N_shape * local_body_forces[0] * jacobian_weight; + noalias(rRHS) += N_shapeY * local_body_forces[1] * jacobian_weight; + } + RotateRHS(rRHS); // rotate to global + + KRATOS_CATCH("") } /***********************************************************************************/ From e0d5b4f93290e2b0a1fa49db32e988ad4d83c3c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Tue, 5 Nov 2024 14:50:28 +0100 Subject: [PATCH 14/32] registering elements --- .../structural_mechanics_application.cpp | 4 ++++ .../structural_mechanics_application.h | 4 +++- 2 files changed, 7 insertions(+), 1 deletion(-) diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp index fae3cabda680..faee6fe9d4d1 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp @@ -65,6 +65,8 @@ KratosStructuralMechanicsApplication::KratosStructuralMechanicsApplication() mTrussElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), mTrussLinearElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), mCableElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), + mLinearTrussElement2D2N(0, Element::GeometryType::Pointer(new Line2D2(Element::GeometryType::PointsArrayType(2)))), + mLinearTrussElement2D3N(0, Element::GeometryType::Pointer(new Line2D3(Element::GeometryType::PointsArrayType(3)))), // Adding the beam elements mCrBeamElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), mCrLinearBeamElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), @@ -525,6 +527,8 @@ void KratosStructuralMechanicsApplication::Register() { //Register the truss element KRATOS_REGISTER_ELEMENT("TrussElement3D2N", mTrussElement3D2N) KRATOS_REGISTER_ELEMENT("TrussLinearElement3D2N", mTrussLinearElement3D2N) + KRATOS_REGISTER_ELEMENT("LinearTrussElement2D2N", mLinearTrussElement2D2N) + KRATOS_REGISTER_ELEMENT("LinearTrussElement2D3N", mLinearTrussElement2D3N) KRATOS_REGISTER_ELEMENT("CableElement3D2N", mCableElement3D2N) // Register the beam element diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.h b/applications/StructuralMechanicsApplication/structural_mechanics_application.h index ca4b26bd05cd..1ab192b764d9 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.h +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.h @@ -281,6 +281,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) KratosStructuralMechanicsAppl const LinearTimoshenkoBeamElement2D2N mLinearTimoshenkoBeamElement2D2N; const LinearTimoshenkoBeamElement2D3N mLinearTimoshenkoBeamElement2D3N; const LinearTimoshenkoCurvedBeamElement2D3N mLinearTimoshenkoCurvedBeamElement2D3N; + const LinearTrussElement2D<2> mLinearTrussElement2D2N; + const LinearTrussElement2D<3> mLinearTrussElement2D3N; // Adding the shells elements @@ -421,7 +423,7 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) KratosStructuralMechanicsAppl const AxisymUpdatedLagrangian mAxisymUpdatedLagrangian2D9N; // Adding the bushing element - const BushingElement mBushingElement3D2N; + const BushingElement mBushingElement3D2N; // Adding the spring damper element const SpringDamperElement<2> mSpringDamperElement2D; From cfa8179631e89c34c9b786a0749c24df30f95fbf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Tue, 5 Nov 2024 16:24:13 +0100 Subject: [PATCH 15/32] minor --- .../truss_elements/linear_truss_element_2D.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index bc38d25248a3..095d36a519ce 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -53,6 +53,7 @@ void LinearTrussElement2D::Initialize(const ProcessInfo& rCurrentProces mConstitutiveLawVector.resize(r_integration_points.size()); InitializeMaterial(); } + KRATOS_CATCH("") } @@ -294,6 +295,7 @@ void LinearTrussElement2D::CalculateLocalSystem( ) { KRATOS_TRY; + const auto &r_props = GetProperties(); const auto &r_geometry = GetGeometry(); @@ -344,7 +346,7 @@ void LinearTrussElement2D::CalculateLocalSystem( strain_vector[0] = inner_prod(B, nodal_values); - mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix noalias(rLHS) += outer_prod(B, B) * constitutive_matrix(0, 0) * jacobian_weight; noalias(rRHS) -= B * stress_vector[0] * jacobian_weight; @@ -409,7 +411,7 @@ void LinearTrussElement2D::CalculateLeftHandSide( GetFirstDerivativesShapeFunctionsValues(B, length, xi); strain_vector[0] = inner_prod(B, nodal_values); - mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix noalias(rLHS) += outer_prod(B, B) * constitutive_matrix(0, 0) * jacobian_weight; } @@ -473,7 +475,7 @@ void LinearTrussElement2D::CalculateRightHandSide( strain_vector[0] = inner_prod(B, nodal_values); - mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); // fills stress and const. matrix + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix noalias(rRHS) -= B * stress_vector[0] * jacobian_weight; From 5a8f3a546819080b431a735e0d6536489876a214 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Wed, 6 Nov 2024 08:46:41 +0100 Subject: [PATCH 16/32] minor --- .../truss_elements/linear_truss_element_2D.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 095d36a519ce..735b6fea0b80 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -171,8 +171,8 @@ void LinearTrussElement2D::GetShapeFunctionsValues( rN[2] = 0.5 * (1.0 + xi); } else { // 3N rN[0] = 0.5 * xi * (xi - 1.0); - rN[2] = (1.0 - std::pow(xi, 2)); - rN[4] = 0.5 * xi * (xi + 1.0); + rN[2] = 0.5 * xi * (xi + 1.0); + rN[4] = (1.0 - std::pow(xi, 2)); } } @@ -195,8 +195,8 @@ void LinearTrussElement2D::GetShapeFunctionsValuesY( rN[3] = 0.5 * (1.0 + xi); } else { // 3N rN[1] = 0.5 * xi * (xi - 1.0); - rN[3] = (1.0 - std::pow(xi, 2)); - rN[5] = 0.5 * xi * (xi + 1.0); + rN[3] = 0.5 * xi * (xi + 1.0); + rN[5] = (1.0 - std::pow(xi, 2)); } } @@ -220,8 +220,8 @@ void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( rdN_dX[2] = inverse_l; } else { // 3N rdN_dX[0] = xi - 0.5; - rdN_dX[2] = -2.0 * xi; - rdN_dX[4] = xi + 0.5; + rdN_dX[2] = xi + 0.5; + rdN_dX[4] = -2.0 * xi; rdN_dX *= 2.0 / Length; // The Jacobian } } From 24f27731ba2b67a38d5bf1efc236d702f7e65d21 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Wed, 6 Nov 2024 10:46:54 +0100 Subject: [PATCH 17/32] remove space --- .../custom_elements/truss_elements/linear_truss_element_2D.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 735b6fea0b80..d12736a3ceb8 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -473,7 +473,6 @@ void LinearTrussElement2D::CalculateRightHandSide( GetShapeFunctionsValuesY(N_shapeY, length, xi); GetFirstDerivativesShapeFunctionsValues(B, length, xi); - strain_vector[0] = inner_prod(B, nodal_values); mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix From 10232b388582faac472a14fbb18257f386435469 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Wed, 6 Nov 2024 12:50:08 +0100 Subject: [PATCH 18/32] prints --- .../linear_truss_element_2D.cpp | 62 +++++++++++++++++++ 1 file changed, 62 insertions(+) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index d12736a3ceb8..9713fb0c005a 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -576,7 +576,69 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( const ProcessInfo& rProcessInfo ) { + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + rOutput.resize(integration_points.size()); + + if (rVariable == AXIAL_FORCE) { + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + const double area = r_props[CROSS_AREA]; + + 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(); + + // Let's initialize the cl values + VectorType strain_vector(1), stress_vector(1); + strain_vector.clear(); + cl_values.SetStrainVector(strain_vector); + cl_values.SetStressVector(stress_vector); + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); + SystemSizeBoundedArrayType B; + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const double xi = integration_points[IP].X(); + + strain_vector[0] = inner_prod(B, nodal_values); + + mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); + rOutput[IP] = cl_values.GetStressVector()[0] * area; + } + } else if (rVariable == AXIAL_STRAIN) { + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + const double area = r_props[CROSS_AREA]; + + 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(); + + // Let's initialize the cl values + VectorType strain_vector(1), stress_vector(1); + strain_vector.clear(); + cl_values.SetStrainVector(strain_vector); + cl_values.SetStressVector(stress_vector); + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); + + SystemSizeBoundedArrayType B; + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const double xi = integration_points[IP].X(); + + rOutput[IP] = inner_prod(B, nodal_values); + } + } } /***********************************************************************************/ From 69d859b62f0a2b5cfd022a2822fc2e46d0bbee4c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Wed, 6 Nov 2024 13:02:51 +0100 Subject: [PATCH 19/32] missing C --- .../truss_elements/linear_truss_element_2D.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 9713fb0c005a..1037f95f1605 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -593,9 +593,11 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( // Let's initialize the cl values VectorType strain_vector(1), stress_vector(1); + MatrixType C(1,1); strain_vector.clear(); cl_values.SetStrainVector(strain_vector); cl_values.SetStressVector(stress_vector); + cl_values.SetConstitutiveMatrix(C); SystemSizeBoundedArrayType nodal_values(SystemSize); GetNodalValuesVector(nodal_values); @@ -607,7 +609,7 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( strain_vector[0] = inner_prod(B, nodal_values); - mConstitutiveLawVector[IP]->CalculateMaterialResponseCauchy(cl_values); + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); rOutput[IP] = cl_values.GetStressVector()[0] * area; } } else if (rVariable == AXIAL_STRAIN) { @@ -624,9 +626,11 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( // Let's initialize the cl values VectorType strain_vector(1), stress_vector(1); + MatrixType C(1,1); strain_vector.clear(); cl_values.SetStrainVector(strain_vector); cl_values.SetStressVector(stress_vector); + cl_values.SetConstitutiveMatrix(C); SystemSizeBoundedArrayType nodal_values(SystemSize); GetNodalValuesVector(nodal_values); @@ -635,7 +639,6 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( // Loop over the integration points for (SizeType IP = 0; IP < integration_points.size(); ++IP) { const double xi = integration_points[IP].X(); - rOutput[IP] = inner_prod(B, nodal_values); } } From f45077bdfa7ebafe0ec3d5c404be4d6866614eaf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Wed, 6 Nov 2024 13:11:26 +0100 Subject: [PATCH 20/32] good print now --- .../truss_elements/linear_truss_element_2D.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 1037f95f1605..dd4ff0947011 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -587,7 +587,7 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( 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); + r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true); const double length = CalculateLength(); @@ -606,6 +606,7 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( // Loop over the integration points for (SizeType IP = 0; IP < integration_points.size(); ++IP) { const double xi = integration_points[IP].X(); + GetFirstDerivativesShapeFunctionsValues(B, length, xi); strain_vector[0] = inner_prod(B, nodal_values); @@ -620,7 +621,7 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( 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); + r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true); const double length = CalculateLength(); @@ -639,6 +640,7 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( // Loop over the integration points for (SizeType IP = 0; IP < integration_points.size(); ++IP) { const double xi = integration_points[IP].X(); + GetFirstDerivativesShapeFunctionsValues(B, length, xi); rOutput[IP] = inner_prod(B, nodal_values); } } From acf22ab33839e43bc06304dd62023bc71047d709 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Wed, 6 Nov 2024 17:10:52 +0100 Subject: [PATCH 21/32] adding 2d2N test input --- .../2D2N/linear_truss_2d2N_test.mdpa | 54 +++++++++++ .../linear_truss_2d2N_test_materials.json | 18 ++++ .../linear_truss_2d2N_test_parameters.json | 95 +++++++++++++++++++ .../structural_mechanics_test_factory.py | 3 + .../test_StructuralMechanicsApplication.py | 2 + 5 files changed, 172 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test.mdpa create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_materials.json create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test.mdpa b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test.mdpa new file mode 100644 index 000000000000..a8b44179f2a8 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test.mdpa @@ -0,0 +1,54 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties +Begin Nodes + 1 0.0000000000 0.0000000000 0.0000000000 + 2 1.0000000000 1.0000000000 0.0000000000 + 3 2.0000000000 0.0000000000 0.0000000000 +End Nodes + + +Begin Elements LinearTrussElement2D2N// GUI group identifier: Truss Auto1 + 1 0 1 2 + 2 0 2 3 +End Elements + +Begin Conditions PointLoadCondition2D1N// GUI group identifier: Load on points Auto1 + 2 0 2 +End Conditions + +Begin SubModelPart Parts_Truss_Truss_Auto1 // Group Truss Auto1 // Subtree Parts_Truss + Begin SubModelPartNodes + 1 + 2 + 3 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + End SubModelPartElements +End SubModelPart +Begin SubModelPart DISPLACEMENT_Displacement_Auto1 // Group Displacement Auto1 // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 3 + End SubModelPartNodes +End SubModelPart +Begin SubModelPart SelfWeight2D_Self_weight_Auto1 // Group Self weight Auto1 // Subtree SelfWeight2D + Begin SubModelPartNodes + 1 + 2 + 3 + End SubModelPartNodes +End SubModelPart +Begin SubModelPart PointLoad2D_Load_on_points_Auto1 // Group Load on points Auto1 // Subtree PointLoad2D + Begin SubModelPartNodes + 2 + End SubModelPartNodes + Begin SubModelPartConditions + 2 + End SubModelPartConditions +End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_materials.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_materials.json new file mode 100644 index 000000000000..2e6c0d37a4ce --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_materials.json @@ -0,0 +1,18 @@ +{ + "properties" : [{ + "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + "properties_id" : 1, + "Material" : { + "constitutive_law" : { + "name" : "TrussConstitutiveLaw" + }, + "Variables" : { + "DENSITY" : 7850.0, + "YOUNG_MODULUS" : 200000000000.0, + "CROSS_AREA" : 0.001, + "TRUSS_PRESTRESS_PK2" : 0.0 + }, + "Tables" : null + } + }] +} diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json new file mode 100644 index 000000000000..83f1bb871615 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json @@ -0,0 +1,95 @@ +{ + "problem_data" : { + "problem_name" : "2d2N", + "parallel_type" : "OpenMP", + "echo_level" : 1, + "start_time" : 0.0, + "end_time" : 1.0 + }, + "solver_settings" : { + "time_stepping" : { + "time_step" : 1.1 + }, + "solver_type" : "Static", + "model_part_name" : "Structure", + "domain_size" : 2, + "echo_level" : 1, + "analysis_type" : "non_linear", + "model_import_settings" : { + "input_type" : "mdpa", + "input_filename" : "2d2N" + }, + "material_import_settings" : { + "materials_filename" : "StructuralMaterials.json" + }, + "line_search" : false, + "convergence_criterion" : "residual_criterion", + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1e-9, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1e-9, + "max_iteration" : 10, + "use_old_stiffness_in_first_iteration" : false, + "rotation_dofs" : false, + "volumetric_strain_dofs" : false + }, + "processes" : { + "constraints_process_list" : [{ + "python_module" : "assign_vector_variable_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "AssignVectorVariableProcess", + "Parameters" : { + "model_part_name" : "Structure.DISPLACEMENT_Displacement_Auto1", + "variable_name" : "DISPLACEMENT", + "interval" : [0.0,"End"], + "constrained" : [true,true,true], + "value" : [0.0,0.0,0.0] + } + }], + "loads_process_list" : [{ + "python_module" : "assign_vector_by_direction_process", + "kratos_module" : "KratosMultiphysics", + "check" : "DirectorVectorNonZero direction", + "process_name" : "AssignVectorByDirectionProcess", + "Parameters" : { + "model_part_name" : "Structure.SelfWeight2D_Self_weight_Auto1", + "variable_name" : "VOLUME_ACCELERATION", + "interval" : [0.0,"End"], + "constrained" : false, + "modulus" : 10.0, + "direction" : [0.0,-1.0,0.0] + } + },{ + "python_module" : "assign_vector_by_direction_to_condition_process", + "kratos_module" : "KratosMultiphysics", + "check" : "DirectorVectorNonZero direction", + "process_name" : "AssignVectorByDirectionToConditionProcess", + "Parameters" : { + "model_part_name" : "Structure.PointLoad2D_Load_on_points_Auto1", + "variable_name" : "POINT_LOAD", + "interval" : [0.0,"End"], + "modulus" : 10.0, + "direction" : [0.0,-1.0,0.0] + } + }], + "list_other_processes" : [], + "_json_output_process" : [ + { + "python_module" : "json_output_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "JsonOutputProcess", + "Parameters" : { + "output_variables" : ["DISPLACEMENT"], + "gauss_points_output_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], + "output_file_name" : "LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json", + "model_part_name" : "Structure.test", + "time_frequency" : 0.1 + } + }] + }, + "output_processes" : { + "gid_output" : [], + "vtk_output" : [] + }, + "analysis_stage" : "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis" +} diff --git a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py index 2dc10830087d..652458a6d978 100644 --- a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py +++ b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py @@ -76,6 +76,9 @@ def tearDown(self): with KratosUnittest.WorkFolderScope(".", __file__): self.test.Finalize() +class LinearTruss2D2NTest(StructuralMechanicsTestFactory): + file_name = "LinearTruss2D/2D2N/linear_truss_2d2N_test" + class TimoshenkoBeam2D2NTest(StructuralMechanicsTestFactory): file_name = "TimoshenkoBeams/2D2N/timoshenko_beam_2d2N_test" diff --git a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py index 870aa1001fb8..3fd59a0a2560 100644 --- a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py +++ b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py @@ -96,6 +96,7 @@ ##### NIGHTLY TESTS ##### # Patch test Small Displacements +from structural_mechanics_test_factory import LinearTruss2D2NTest as TLinearTruss2D2NTest from structural_mechanics_test_factory import TimoshenkoBeam2D2NTest as TTimoshenkoBeam2D2NTest from structural_mechanics_test_factory import TimoshenkoBeam2D3NTest as TTimoshenkoBeam2D3NTest from structural_mechanics_test_factory import TimoshenkoCurvedBeam2D3NTest as TTimoshenkoCurvedBeam2D3NTest @@ -344,6 +345,7 @@ def AssembleTestSuites(): ### Adding Nightly Tests # Patch test Small Displacements + smallSuite.addTest(TLinearTruss2D2NTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam2D2NTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam2D3NTest('test_execution')) smallSuite.addTest(TTimoshenkoCurvedBeam2D3NTest('test_execution')) From c4a2f632f4d6b1caebff2628ef03bb9b89c662f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Fri, 8 Nov 2024 10:22:04 +0100 Subject: [PATCH 22/32] adding test --- .../2D2N/linear_truss_2d2N_test_parameters.json | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json index 83f1bb871615..05984aa1ba35 100644 --- a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json @@ -17,10 +17,10 @@ "analysis_type" : "non_linear", "model_import_settings" : { "input_type" : "mdpa", - "input_filename" : "2d2N" + "input_filename" : "LinearTruss2D/2D2N/linear_truss_2d2N_test" }, "material_import_settings" : { - "materials_filename" : "StructuralMaterials.json" + "materials_filename" : "LinearTruss2D/2D2N/linear_truss_2d2N_test_materials.json" }, "line_search" : false, "convergence_criterion" : "residual_criterion", @@ -82,8 +82,8 @@ "output_variables" : ["DISPLACEMENT"], "gauss_points_output_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], "output_file_name" : "LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json", - "model_part_name" : "Structure.test", - "time_frequency" : 0.1 + "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + "time_frequency" : 1.0 } }] }, From 0b3eefd6189ab021326282079b1bf0ac230e8d71 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Fri, 8 Nov 2024 17:37:58 +0100 Subject: [PATCH 23/32] results --- .../2D2N/linear_truss_2d2N_test_results.json | 62 +++++++++++++++++++ 1 file changed, 62 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json new file mode 100644 index 000000000000..20084809d889 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json @@ -0,0 +1,62 @@ +{ + "TIME": [ + 1.1 + ], + "NODE_1": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_2": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -8.557106781186548e-07 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_3": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "ELEMENT_1": { + "AXIAL_FORCE": { + "0": [ + -85.57106781186548 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -4.278553390593274e-07 + ] + } + }, + "ELEMENT_2": { + "AXIAL_FORCE": { + "0": [ + -85.57106781186548 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -4.278553390593274e-07 + ] + } + } +} \ No newline at end of file From 1a13a56c5d29f6b40985b97ac776e7e430ab79a9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Fri, 8 Nov 2024 17:38:09 +0100 Subject: [PATCH 24/32] improve truss cl --- .../custom_constitutive/truss_constitutive_law.cpp | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp b/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp index 9451674b12bf..8cce462edcad 100644 --- a/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp +++ b/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp @@ -156,9 +156,12 @@ double TrussConstitutiveLaw::CalculateStressElastic( CalculateValue(rParameterValues,TANGENT_MODULUS,tangent_modulus); const double current_stress = tangent_modulus*current_strain[0]; - Matrix& r_const_matrix = rParameterValues.GetConstitutiveMatrix(); - r_const_matrix.resize(1, 1); - r_const_matrix(0, 0) = tangent_modulus; + + if (rParameterValues.GetOptions().Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)) { + Matrix& r_const_matrix = rParameterValues.GetConstitutiveMatrix(); + r_const_matrix.resize(1, 1); + r_const_matrix(0, 0) = tangent_modulus; + } return current_stress; } From 051816410c5d32deccf755c9d95ab7a8c08ca633 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Fri, 8 Nov 2024 17:42:54 +0100 Subject: [PATCH 25/32] add check --- .../linear_truss_2d2N_test_parameters.json | 40 +++++++++++++------ 1 file changed, 27 insertions(+), 13 deletions(-) diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json index 05984aa1ba35..39990a9fe310 100644 --- a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json @@ -73,19 +73,33 @@ } }], "list_other_processes" : [], - "_json_output_process" : [ - { - "python_module" : "json_output_process", - "kratos_module" : "KratosMultiphysics", - "process_name" : "JsonOutputProcess", - "Parameters" : { - "output_variables" : ["DISPLACEMENT"], - "gauss_points_output_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], - "output_file_name" : "LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json", - "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", - "time_frequency" : 1.0 - } - }] + "json_check_process" : [ + { + "python_module" : "from_json_check_result_process", + "kratos_module" : "KratosMultiphysics", + "help" : "", + "process_name" : "FromJsonCheckResultProcess", + "Parameters" : { + "check_variables" : ["DISPLACEMENT"], + "gauss_points_check_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], + "input_file_name" : "LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json", + "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + "time_frequency" : 1.0 + } + }] + //"_json_output_process" : [ + //{ + // "python_module" : "json_output_process", + // "kratos_module" : "KratosMultiphysics", + // "process_name" : "JsonOutputProcess", + // "Parameters" : { + // "output_variables" : ["DISPLACEMENT"], + // "gauss_points_output_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], + // "output_file_name" : "LinearTruss2D/2D2N/linear_truss_2d2N_test_results.json", + // "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + // "time_frequency" : 1.0 + // } + //}] }, "output_processes" : { "gid_output" : [], From dcd20550376a331d97ecf9d3def68a5e1ce7dfb5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Fri, 8 Nov 2024 18:06:06 +0100 Subject: [PATCH 26/32] adding 2D3N test --- .../2D3N/linear_truss_2D3N_test_results.json | 108 ++++++++++++++++ .../2D3N/linear_truss_2d3N_test.mdpa | 87 +++++++++++++ .../linear_truss_2d3N_test_materials.json | 33 +++++ .../linear_truss_2d3N_test_parameters.json | 120 ++++++++++++++++++ .../structural_mechanics_test_factory.py | 3 + .../test_StructuralMechanicsApplication.py | 2 + 6 files changed, 353 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_materials.json create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json new file mode 100644 index 000000000000..fd16e6c18ec1 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json @@ -0,0 +1,108 @@ +{ + "TIME": [ + 1.1 + ], + "NODE_1": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_2": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -1.1332500897343743e-06 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_3": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_4": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -8.66952324223082e-07 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_5": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -8.66952324223082e-07 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "ELEMENT_1": { + "AXIAL_FORCE": { + "0": [ + -182.68262320588295 + ], + "1": [ + -43.96739474099189 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -9.134131160294146e-07 + ], + "1": [ + -2.1983697370495944e-07 + ] + } + }, + "ELEMENT_2": { + "AXIAL_FORCE": { + "0": [ + -43.96739474099189 + ], + "1": [ + -182.68262320588295 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -2.1983697370495944e-07 + ], + "1": [ + -9.134131160294146e-07 + ] + } + }, + "ELEMENT_3": { + "AXIAL_FORCE": { + "0": [ + 0.0 + ] + }, + "AXIAL_STRAIN": { + "0": [ + 0.0 + ] + } + } +} \ No newline at end of file diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa new file mode 100644 index 000000000000..55574f9a64eb --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa @@ -0,0 +1,87 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties +Begin Nodes + 1 0.0 0.0 0.0 + 2 1.0 1.0 0.0 + 3 2.0 0.0 0.0 + 4 0.5 0.5 0.0 + 5 1.5 0.5 0.0 +End Nodes + + +Begin Elements LinearTrussElement2D3N// GUI group identifier: Truss Auto1 + 1 0 1 2 4 + 2 0 2 3 5 +End Elements + +Begin Elements LinearTrussElement2D2N// GUI group identifier: Truss Auto2 + 3 0 4 5 +End Elements + +Begin Conditions PointLoadCondition2D1N// GUI group identifier: Load on points Auto1 + 2 0 2 +End Conditions + +Begin SubModelPart Parts_Truss_Truss_Auto1 // Group Truss Auto1 // Subtree Parts_Truss + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart truss2N // Group Truss Auto1 // Subtree Parts_Truss + Begin SubModelPartNodes + 4 + 5 + End SubModelPartNodes + Begin SubModelPartElements + 3 + End SubModelPartElements +End SubModelPart + + + + + +Begin SubModelPart DISPLACEMENT_Displacement_Auto1 + Begin SubModelPartNodes + 1 + 3 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Displacement_Auto2 + Begin SubModelPartNodes + 4 + 5 + End SubModelPartNodes +End SubModelPart + +Begin SubModelPart SelfWeight2D_Self_weight_Auto1 // Group Self weight Auto1 // Subtree SelfWeight2D + Begin SubModelPartNodes + 1 + 2 + 3 + 4 + 5 + End SubModelPartNodes +End SubModelPart +Begin SubModelPart PointLoad2D_Load_on_points_Auto1 // Group Load on points Auto1 // Subtree PointLoad2D + Begin SubModelPartNodes + 2 + End SubModelPartNodes + Begin SubModelPartConditions + 2 + End SubModelPartConditions +End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_materials.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_materials.json new file mode 100644 index 000000000000..9dccc3859cd0 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_materials.json @@ -0,0 +1,33 @@ +{ + "properties" : [{ + "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + "properties_id" : 1, + "Material" : { + "constitutive_law" : { + "name" : "TrussConstitutiveLaw" + }, + "Variables" : { + "DENSITY" : 7850.0, + "YOUNG_MODULUS" : 200000000000.0, + "CROSS_AREA" : 0.001, + "TRUSS_PRESTRESS_PK2" : 0.0 + }, + "Tables" : null + } + },{ + "model_part_name" : "Structure.truss2N", + "properties_id" : 1, + "Material" : { + "constitutive_law" : { + "name" : "TrussConstitutiveLaw" + }, + "Variables" : { + "DENSITY" : 7850.0, + "YOUNG_MODULUS" : 200000000000.0, + "CROSS_AREA" : 0.001, + "TRUSS_PRESTRESS_PK2" : 0.0 + }, + "Tables" : null + } + }] +} diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json new file mode 100644 index 000000000000..bd6f4118fe67 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json @@ -0,0 +1,120 @@ +{ + "problem_data" : { + "problem_name" : "2D3N", + "parallel_type" : "OpenMP", + "echo_level" : 1, + "start_time" : 0.0, + "end_time" : 1.0 + }, + "solver_settings" : { + "time_stepping" : { + "time_step" : 1.1 + }, + "solver_type" : "Static", + "model_part_name" : "Structure", + "domain_size" : 2, + "echo_level" : 1, + "analysis_type" : "non_linear", + "model_import_settings" : { + "input_type" : "mdpa", + "input_filename" : "LinearTruss2D/2D3N/linear_truss_2D3N_test" + }, + "material_import_settings" : { + "materials_filename" : "LinearTruss2D/2D3N/linear_truss_2D3N_test_materials.json" + }, + "line_search" : false, + "convergence_criterion" : "residual_criterion", + "displacement_relative_tolerance" : 0.0001, + "displacement_absolute_tolerance" : 1e-9, + "residual_relative_tolerance" : 0.0001, + "residual_absolute_tolerance" : 1e-9, + "max_iteration" : 10, + "use_old_stiffness_in_first_iteration" : false, + "rotation_dofs" : false, + "volumetric_strain_dofs" : false + }, + "processes" : { + "constraints_process_list" : [{ + "python_module" : "assign_vector_variable_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "AssignVectorVariableProcess", + "Parameters" : { + "model_part_name" : "Structure.DISPLACEMENT_Displacement_Auto1", + "variable_name" : "DISPLACEMENT", + "interval" : [0.0,"End"], + "constrained" : [true,true,true], + "value" : [0.0,0.0,0.0] + } + },{ + "python_module" : "assign_vector_variable_process", + "kratos_module" : "KratosMultiphysics", + "process_name" : "AssignVectorVariableProcess", + "Parameters" : { + "model_part_name" : "Structure.DISPLACEMENT_Displacement_Auto2", + "variable_name" : "DISPLACEMENT", + "interval" : [0.0,"End"], + "constrained" : [true,false,true], + "value" : [0.0,null,0.0] + } + }], + "loads_process_list" : [{ + "python_module" : "assign_vector_by_direction_process", + "kratos_module" : "KratosMultiphysics", + "check" : "DirectorVectorNonZero direction", + "process_name" : "AssignVectorByDirectionProcess", + "Parameters" : { + "model_part_name" : "Structure.SelfWeight2D_Self_weight_Auto1", + "variable_name" : "VOLUME_ACCELERATION", + "interval" : [0.0,"End"], + "constrained" : false, + "modulus" : 10.0, + "direction" : [0.0,-1.0,0.0] + } + },{ + "python_module" : "assign_vector_by_direction_to_condition_process", + "kratos_module" : "KratosMultiphysics", + "check" : "DirectorVectorNonZero direction", + "process_name" : "AssignVectorByDirectionToConditionProcess", + "Parameters" : { + "model_part_name" : "Structure.PointLoad2D_Load_on_points_Auto1", + "variable_name" : "POINT_LOAD", + "interval" : [0.0,"End"], + "modulus" : 10.0, + "direction" : [0.0,-1.0,0.0] + } + }], + "list_other_processes" : [], + "json_check_process" : [ + { + "python_module" : "from_json_check_result_process", + "kratos_module" : "KratosMultiphysics", + "help" : "", + "process_name" : "FromJsonCheckResultProcess", + "Parameters" : { + "check_variables" : ["DISPLACEMENT"], + "gauss_points_check_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], + "input_file_name" : "LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json", + "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + "time_frequency" : 1.0 + } + }] + //"_json_output_process" : [ + //{ + // "python_module" : "json_output_process", + // "kratos_module" : "KratosMultiphysics", + // "process_name" : "JsonOutputProcess", + // "Parameters" : { + // "output_variables" : ["DISPLACEMENT"], + // "gauss_points_output_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], + // "output_file_name" : "LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json", + // "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + // "time_frequency" : 1.0 + // } + //}] + }, + "output_processes" : { + "gid_output" : [], + "vtk_output" : [] + }, + "analysis_stage" : "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis" +} diff --git a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py index 652458a6d978..e38c1a15ceea 100644 --- a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py +++ b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py @@ -79,6 +79,9 @@ def tearDown(self): class LinearTruss2D2NTest(StructuralMechanicsTestFactory): file_name = "LinearTruss2D/2D2N/linear_truss_2d2N_test" +class LinearTruss2D3NTest(StructuralMechanicsTestFactory): + file_name = "LinearTruss2D/2D3N/linear_truss_2d3N_test" + class TimoshenkoBeam2D2NTest(StructuralMechanicsTestFactory): file_name = "TimoshenkoBeams/2D2N/timoshenko_beam_2d2N_test" diff --git a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py index 3fd59a0a2560..e29dc9312a48 100644 --- a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py +++ b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py @@ -97,6 +97,7 @@ ##### NIGHTLY TESTS ##### # Patch test Small Displacements from structural_mechanics_test_factory import LinearTruss2D2NTest as TLinearTruss2D2NTest +from structural_mechanics_test_factory import LinearTruss2D3NTest as TLinearTruss2D3NTest from structural_mechanics_test_factory import TimoshenkoBeam2D2NTest as TTimoshenkoBeam2D2NTest from structural_mechanics_test_factory import TimoshenkoBeam2D3NTest as TTimoshenkoBeam2D3NTest from structural_mechanics_test_factory import TimoshenkoCurvedBeam2D3NTest as TTimoshenkoCurvedBeam2D3NTest @@ -346,6 +347,7 @@ def AssembleTestSuites(): ### Adding Nightly Tests # Patch test Small Displacements smallSuite.addTest(TLinearTruss2D2NTest('test_execution')) + smallSuite.addTest(TLinearTruss2D3NTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam2D2NTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam2D3NTest('test_execution')) smallSuite.addTest(TTimoshenkoCurvedBeam2D3NTest('test_execution')) From 3c49879bd21339a5a1bf866cbbe037642b905a98 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Sat, 9 Nov 2024 17:44:15 +0100 Subject: [PATCH 27/32] unused var --- .../custom_elements/truss_elements/linear_truss_element_2D.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index dd4ff0947011..2b6d6da4b054 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -616,7 +616,6 @@ void LinearTrussElement2D::CalculateOnIntegrationPoints( } else if (rVariable == AXIAL_STRAIN) { const auto &r_props = GetProperties(); const auto &r_geometry = GetGeometry(); - const double area = r_props[CROSS_AREA]; ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo); auto &r_cl_options = cl_values.GetOptions(); From 21be5ba1c17db5d8ea5b69249fd4bbe8c4b20972 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alejandro=20Cornejo=20Vel=C3=A1zquez?= Date: Mon, 11 Nov 2024 10:08:48 +0100 Subject: [PATCH 28/32] avoid repetition --- .../truss_elements/linear_truss_element_2D.h | 25 +++++++++++-------- 1 file changed, 14 insertions(+), 11 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h index e070dd6ae2f1..40b7503e1119 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -82,23 +82,14 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D // Constructor using an array of nodes LinearTrussElement2D(IndexType NewId, GeometryType::Pointer pGeometry) : Element(NewId, pGeometry) { - if constexpr (NNodes == 2) { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; - } else { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; - } + mThisIntegrationMethod = CalculateIntegrationMethod(); } // Constructor using an array of nodes with properties LinearTrussElement2D(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) : Element(NewId,pGeometry,pProperties) { - // This is needed to prevent uninitialised integration method in inactive elements - if constexpr (NNodes == 2) { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; - } else { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; - } + mThisIntegrationMethod = CalculateIntegrationMethod(); } // Copy constructor @@ -137,6 +128,18 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) LinearTrussElement2D return StructuralMechanicsElementUtilities::GetReferenceRotationAngle2D2NBeam(GetGeometry()); } + /** + * @brief This method returns the integration method depending on the Number of Nodes + */ + IntegrationMethod CalculateIntegrationMethod() + { + if constexpr (NNodes == 2) { + return GeometryData::IntegrationMethod::GI_GAUSS_1; + } else { + return GeometryData::IntegrationMethod::GI_GAUSS_2; + } + } + /** * @brief Returns a n component vector including the values of the DoFs * in LOCAL beam axes From 19bae45b6b6b5c02e20df7201a7fbb334e0e4a28 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Wed, 13 Nov 2024 08:16:53 +0100 Subject: [PATCH 29/32] reordering register of FEs --- .../structural_mechanics_application.cpp | 2 +- .../structural_mechanics_application.h | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp index 27b4a9f27ac3..e5d2772a256b 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp @@ -530,9 +530,9 @@ void KratosStructuralMechanicsApplication::Register() { //Register the truss element KRATOS_REGISTER_ELEMENT("TrussElement3D2N", mTrussElement3D2N) KRATOS_REGISTER_ELEMENT("TrussLinearElement3D2N", mTrussLinearElement3D2N) + KRATOS_REGISTER_ELEMENT("CableElement3D2N", mCableElement3D2N) KRATOS_REGISTER_ELEMENT("LinearTrussElement2D2N", mLinearTrussElement2D2N) KRATOS_REGISTER_ELEMENT("LinearTrussElement2D3N", mLinearTrussElement2D3N) - KRATOS_REGISTER_ELEMENT("CableElement3D2N", mCableElement3D2N) // Register the beam element KRATOS_REGISTER_ELEMENT("CrBeamElement3D2N", mCrBeamElement3D2N) diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.h b/applications/StructuralMechanicsApplication/structural_mechanics_application.h index 1ab192b764d9..ff17c6688238 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.h +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.h @@ -272,6 +272,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) KratosStructuralMechanicsAppl const TrussElement3D2N mTrussElement3D2N; const TrussElementLinear3D2N mTrussLinearElement3D2N; const CableElement3D2N mCableElement3D2N; + const LinearTrussElement2D<2> mLinearTrussElement2D2N; + const LinearTrussElement2D<3> mLinearTrussElement2D3N; // Adding the beam element const CrBeamElement3D2N mCrBeamElement3D2N; @@ -281,8 +283,6 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) KratosStructuralMechanicsAppl const LinearTimoshenkoBeamElement2D2N mLinearTimoshenkoBeamElement2D2N; const LinearTimoshenkoBeamElement2D3N mLinearTimoshenkoBeamElement2D3N; const LinearTimoshenkoCurvedBeamElement2D3N mLinearTimoshenkoCurvedBeamElement2D3N; - const LinearTrussElement2D<2> mLinearTrussElement2D2N; - const LinearTrussElement2D<3> mLinearTrussElement2D3N; // Adding the shells elements From 9acbfe58c26d9eb847494bdea0209f26aeb6951e Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Wed, 13 Nov 2024 10:16:43 +0100 Subject: [PATCH 30/32] minor in test --- .../2D3N/linear_truss_2D3N_test_results.json | 108 ------------------ .../2D3N/linear_truss_2d3N_test.mdpa | 27 ++--- .../linear_truss_2d3N_test_parameters.json | 8 +- 3 files changed, 16 insertions(+), 127 deletions(-) delete mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json deleted file mode 100644 index fd16e6c18ec1..000000000000 --- a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json +++ /dev/null @@ -1,108 +0,0 @@ -{ - "TIME": [ - 1.1 - ], - "NODE_1": { - "DISPLACEMENT_X": [ - 0.0 - ], - "DISPLACEMENT_Y": [ - 0.0 - ], - "DISPLACEMENT_Z": [ - 0.0 - ] - }, - "NODE_2": { - "DISPLACEMENT_X": [ - 0.0 - ], - "DISPLACEMENT_Y": [ - -1.1332500897343743e-06 - ], - "DISPLACEMENT_Z": [ - 0.0 - ] - }, - "NODE_3": { - "DISPLACEMENT_X": [ - 0.0 - ], - "DISPLACEMENT_Y": [ - 0.0 - ], - "DISPLACEMENT_Z": [ - 0.0 - ] - }, - "NODE_4": { - "DISPLACEMENT_X": [ - 0.0 - ], - "DISPLACEMENT_Y": [ - -8.66952324223082e-07 - ], - "DISPLACEMENT_Z": [ - 0.0 - ] - }, - "NODE_5": { - "DISPLACEMENT_X": [ - 0.0 - ], - "DISPLACEMENT_Y": [ - -8.66952324223082e-07 - ], - "DISPLACEMENT_Z": [ - 0.0 - ] - }, - "ELEMENT_1": { - "AXIAL_FORCE": { - "0": [ - -182.68262320588295 - ], - "1": [ - -43.96739474099189 - ] - }, - "AXIAL_STRAIN": { - "0": [ - -9.134131160294146e-07 - ], - "1": [ - -2.1983697370495944e-07 - ] - } - }, - "ELEMENT_2": { - "AXIAL_FORCE": { - "0": [ - -43.96739474099189 - ], - "1": [ - -182.68262320588295 - ] - }, - "AXIAL_STRAIN": { - "0": [ - -2.1983697370495944e-07 - ], - "1": [ - -9.134131160294146e-07 - ] - } - }, - "ELEMENT_3": { - "AXIAL_FORCE": { - "0": [ - 0.0 - ] - }, - "AXIAL_STRAIN": { - "0": [ - 0.0 - ] - } - } -} \ No newline at end of file diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa index 55574f9a64eb..938cd7af45d7 100644 --- a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa @@ -8,8 +8,8 @@ Begin Nodes 1 0.0 0.0 0.0 2 1.0 1.0 0.0 3 2.0 0.0 0.0 - 4 0.5 0.5 0.0 - 5 1.5 0.5 0.0 + 4 0.5 0.5 0.0 + 5 1.5 0.5 0.0 End Nodes @@ -19,7 +19,7 @@ Begin Elements LinearTrussElement2D3N// GUI group identifier: Truss Auto1 End Elements Begin Elements LinearTrussElement2D2N// GUI group identifier: Truss Auto2 - 3 0 4 5 + 3 0 4 5 End Elements Begin Conditions PointLoadCondition2D1N// GUI group identifier: Load on points Auto1 @@ -31,8 +31,8 @@ Begin SubModelPart Parts_Truss_Truss_Auto1 // Group Truss Auto1 // Subtree Parts 1 2 3 - 4 - 5 + 4 + 5 End SubModelPartNodes Begin SubModelPartElements 1 @@ -42,18 +42,15 @@ End SubModelPart Begin SubModelPart truss2N // Group Truss Auto1 // Subtree Parts_Truss Begin SubModelPartNodes - 4 - 5 + 4 + 5 End SubModelPartNodes Begin SubModelPartElements - 3 + 3 End SubModelPartElements End SubModelPart - - - Begin SubModelPart DISPLACEMENT_Displacement_Auto1 Begin SubModelPartNodes 1 @@ -63,8 +60,8 @@ End SubModelPart Begin SubModelPart DISPLACEMENT_Displacement_Auto2 Begin SubModelPartNodes - 4 - 5 + 4 + 5 End SubModelPartNodes End SubModelPart @@ -73,8 +70,8 @@ Begin SubModelPart SelfWeight2D_Self_weight_Auto1 // Group Self weight Auto1 // 1 2 3 - 4 - 5 + 4 + 5 End SubModelPartNodes End SubModelPart Begin SubModelPart PointLoad2D_Load_on_points_Auto1 // Group Load on points Auto1 // Subtree PointLoad2D diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json index bd6f4118fe67..df81081fd538 100644 --- a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_parameters.json @@ -17,10 +17,10 @@ "analysis_type" : "non_linear", "model_import_settings" : { "input_type" : "mdpa", - "input_filename" : "LinearTruss2D/2D3N/linear_truss_2D3N_test" + "input_filename" : "LinearTruss2D/2D3N/linear_truss_2d3N_test" }, "material_import_settings" : { - "materials_filename" : "LinearTruss2D/2D3N/linear_truss_2D3N_test_materials.json" + "materials_filename" : "LinearTruss2D/2D3N/linear_truss_2d3N_test_materials.json" }, "line_search" : false, "convergence_criterion" : "residual_criterion", @@ -93,7 +93,7 @@ "Parameters" : { "check_variables" : ["DISPLACEMENT"], "gauss_points_check_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], - "input_file_name" : "LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json", + "input_file_name" : "LinearTruss2D/2D3N/linear_truss_2d3N_test_results.json", "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", "time_frequency" : 1.0 } @@ -106,7 +106,7 @@ // "Parameters" : { // "output_variables" : ["DISPLACEMENT"], // "gauss_points_output_variables" : ["AXIAL_FORCE", "AXIAL_STRAIN"], - // "output_file_name" : "LinearTruss2D/2D3N/linear_truss_2D3N_test_results.json", + // "output_file_name" : "LinearTruss2D/2D3N/linear_truss_2d3N_test_results.json", // "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", // "time_frequency" : 1.0 // } From 345a423346bcf41527e6327bc468c474e6796243 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Wed, 13 Nov 2024 10:35:03 +0100 Subject: [PATCH 31/32] add chack and avoid a repet --- .../truss_elements/linear_truss_element_2D.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp index 2b6d6da4b054..9ac65a72ac9e 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -38,11 +38,7 @@ void LinearTrussElement2D::Initialize(const ProcessInfo& rCurrentProces if (GetProperties().Has(INTEGRATION_ORDER) ) { mThisIntegrationMethod = static_cast(GetProperties()[INTEGRATION_ORDER] - 1); } else { - if constexpr (NNodes == 2) { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_1; - } else { - mThisIntegrationMethod = GeometryData::IntegrationMethod::GI_GAUSS_2; - } + mThisIntegrationMethod = CalculateIntegrationMethod(); } } @@ -675,6 +671,7 @@ int LinearTrussElement2D::Check(const ProcessInfo& rCurrentProcessInfo) KRATOS_TRY; return mConstitutiveLawVector[0]->Check(GetProperties(), GetGeometry(), rCurrentProcessInfo); + KRATOS_ERROR_IF_NOT(GetProperties().Has(CROSS_AREA)) << "CROSS_AREA not defined in the properties" << std::endl; KRATOS_CATCH( "" ); } From 31cf8cb0833759d5c719c364daac3bd2929b4c47 Mon Sep 17 00:00:00 2001 From: AlejandroCornejo Date: Wed, 13 Nov 2024 11:29:41 +0100 Subject: [PATCH 32/32] add new res file 2D3N --- .../2D3N/linear_truss_2d3N_test_results.json | 108 ++++++++++++++++++ 1 file changed, 108 insertions(+) create mode 100644 applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_results.json diff --git a/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_results.json b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_results.json new file mode 100644 index 000000000000..fd16e6c18ec1 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test_results.json @@ -0,0 +1,108 @@ +{ + "TIME": [ + 1.1 + ], + "NODE_1": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_2": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -1.1332500897343743e-06 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_3": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_4": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -8.66952324223082e-07 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_5": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + -8.66952324223082e-07 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "ELEMENT_1": { + "AXIAL_FORCE": { + "0": [ + -182.68262320588295 + ], + "1": [ + -43.96739474099189 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -9.134131160294146e-07 + ], + "1": [ + -2.1983697370495944e-07 + ] + } + }, + "ELEMENT_2": { + "AXIAL_FORCE": { + "0": [ + -43.96739474099189 + ], + "1": [ + -182.68262320588295 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -2.1983697370495944e-07 + ], + "1": [ + -9.134131160294146e-07 + ] + } + }, + "ELEMENT_3": { + "AXIAL_FORCE": { + "0": [ + 0.0 + ] + }, + "AXIAL_STRAIN": { + "0": [ + 0.0 + ] + } + } +} \ No newline at end of file