diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp index 927b8d0806e7..4b5da0b37f12 100644 --- a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element.cpp @@ -711,6 +711,11 @@ void LinearTrussElement::CalculateOnIntegrationPoints( SystemSizeBoundedArrayType B; const double area = GetProperties()[CROSS_AREA]; + double pre_stress = 0.0; + if (GetProperties().Has(TRUSS_PRESTRESS_PK2)) { + pre_stress = GetProperties()[TRUSS_PRESTRESS_PK2]; + } + // Loop over the integration points for (SizeType IP = 0; IP < integration_points.size(); ++IP) { const double xi = integration_points[IP].X(); @@ -720,11 +725,6 @@ void LinearTrussElement::CalculateOnIntegrationPoints( mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); - double pre_stress = 0.0; - if (GetProperties().Has(TRUSS_PRESTRESS_PK2)) { - pre_stress = GetProperties()[TRUSS_PRESTRESS_PK2]; - } - rOutput[IP] = (cl_values.GetStressVector()[0] + pre_stress) * area; } } else if (rVariable == AXIAL_STRAIN) { diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/total_lagrangian_truss_element.cpp b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/total_lagrangian_truss_element.cpp new file mode 100644 index 000000000000..942561fcf990 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/total_lagrangian_truss_element.cpp @@ -0,0 +1,1065 @@ +// KRATOS ___| | | | +// \___ \ __| __| | | __| __| | | __| _` | | +// | | | | | ( | | | | ( | | +// _____/ \__|_| \__,_|\___|\__|\__,_|_| \__,_|_| MECHANICS +// +// License: BSD License +// license: StructuralMechanicsApplication/license.txt +// +// Main authors: Alejandro Cornejo +// +// + +// System includes + +// External includes + +// Project includes + +// Application includes +#include "total_lagrangian_truss_element.h" +#include "structural_mechanics_application_variables.h" + +namespace Kratos +{ + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::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 = CalculateIntegrationMethod(); + } + } + + 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 TotalLagrangianTrussElement::InitializeMaterial() +{ + KRATOS_TRY + const auto &r_props = GetProperties(); + + if (r_props[CONSTITUTIVE_LAW] != nullptr) { + const auto& r_geometry = GetGeometry(); + 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_props, 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 TotalLagrangianTrussElement::Clone( + IndexType NewId, + NodesArrayType const& rThisNodes + ) const +{ + KRATOS_TRY + + TotalLagrangianTrussElement::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 TotalLagrangianTrussElement::EquationIdVector( + EquationIdVectorType& rResult, + const ProcessInfo& rCurrentProcessInfo + ) const +{ + const auto& r_geometry = this->GetGeometry(); + const SizeType number_of_nodes = r_geometry.size(); + IndexType local_index = 0; + + if (rResult.size() != SystemSize) + rResult.resize(SystemSize, false); + + 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(); + rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_Y, xpos + 1).EquationId(); + if constexpr (Dimension == 3) { + rResult[local_index++] = r_geometry[i].GetDof(DISPLACEMENT_Z, xpos + 2).EquationId(); + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetDofList( + DofsVectorType& rElementalDofList, + const ProcessInfo& rCurrentProcessInfo + ) const +{ + KRATOS_TRY; + + const auto& r_geom = GetGeometry(); + const SizeType number_of_nodes = r_geom.size(); + rElementalDofList.resize(SystemSize); + + for (IndexType i = 0; i < number_of_nodes; ++i) { + const SizeType index = i * DofsPerNode; + rElementalDofList[index] = r_geom[i].pGetDof(DISPLACEMENT_X); + rElementalDofList[index + 1] = r_geom[i].pGetDof(DISPLACEMENT_Y); + if constexpr (Dimension == 3) { + rElementalDofList[index + 2] = r_geom[i].pGetDof(DISPLACEMENT_Z); + } + } + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetShapeFunctionsValues( + SystemSizeBoundedArrayType& rN, + const double Length, + const double xi + ) const +{ + if (rN.size() != SystemSize) + rN.resize(SystemSize, false); + + rN.clear(); + array_1d base_N; + noalias(base_N) = GetBaseShapeFunctions(xi); + + if constexpr (Dimension == 2) { + rN[0] = base_N[0]; + rN[2] = base_N[1]; + } else { + rN[0] = base_N[0]; + rN[3] = base_N[1]; + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetShapeFunctionsValuesY( + SystemSizeBoundedArrayType& rN, + const double Length, + const double xi + ) const +{ + if (rN.size() != SystemSize) + rN.resize(SystemSize, false); + + rN.clear(); + array_1d base_N; + noalias(base_N) = GetBaseShapeFunctions(xi); + + if constexpr (Dimension == 2) { + rN[1] = base_N[0]; + rN[3] = base_N[1]; + } else { + rN[1] = base_N[0]; + rN[4] = base_N[1]; + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetShapeFunctionsValuesZ( + SystemSizeBoundedArrayType& rN, + const double Length, + const double xi + ) const +{ + if (rN.size() != SystemSize) + rN.resize(SystemSize, false); + + rN.clear(); + + if constexpr (Dimension == 3) { + array_1d base_N; + noalias(base_N) = GetBaseShapeFunctions(xi); + rN[2] = base_N[0]; + rN[5] = base_N[1]; + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetFirstDerivativesShapeFunctionsValues( + SystemSizeBoundedArrayType& rdN_dX, + const double RefLength, + const double xi + ) const +{ + if (rdN_dX.size() != SystemSize) + rdN_dX.resize(SystemSize, false); + + rdN_dX.clear(); + + Vector coord(3); + Matrix dN_de(NNodes, 1); + coord.clear(); + coord[0] = xi; + GetGeometry().ShapeFunctionsLocalGradients(dN_de, coord); + + if constexpr (Dimension == 2) { + rdN_dX[0] = dN_de(0, 0); + rdN_dX[2] = dN_de(1, 0); + } else { + rdN_dX[0] = dN_de(0, 0); + rdN_dX[3] = dN_de(1, 0); + } + rdN_dX *= 2.0 / RefLength; // The Jacobian +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +BoundedMatrix TotalLagrangianTrussElement::GetFrenetSerretMatrix() const +{ + // The Frenet-Serret in the CURRENT configuration + return StructuralMechanicsElementUtilities::GetFrenetSerretMatrix3D(GetGeometry(), true); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetNodalValuesVector( + SystemSizeBoundedArrayType& rNodalValues +) const +{ + if (rNodalValues.size() != SystemSize) + rNodalValues.resize(SystemSize, false); + + const auto &r_geom = GetGeometry(); + BoundedVector global_values; + BoundedMatrix T; + BoundedMatrix global_size_T; + + if constexpr (Dimension == 2) { + const double angle = GetCurrentAngle(); + // 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]; + } + + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + noalias(rNodalValues) = prod(trans(global_size_T), global_values); + } else { + // 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]; + global_values[i * DofsPerNode + 2] = r_displ[2]; + } + + noalias(T) = GetFrenetSerretMatrix(); // global to local + if constexpr (NNodes == 2) { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor3D2NTruss(T, global_size_T); + } else { + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor3D3NTruss(T, global_size_T); + } + noalias(rNodalValues) = prod(global_size_T, global_values); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +array_1d TotalLagrangianTrussElement::GetLocalAxesBodyForce( + const Element &rElement, + const GeometryType::IntegrationPointsArrayType &rIntegrationPoints, + const IndexType PointNumber + ) const +{ + const auto body_force = StructuralMechanicsElementUtilities::GetBodyForce(*this, rIntegrationPoints, PointNumber); + array_1d local_body_force; + + if constexpr (Dimension == 2) { + const double angle = GetCurrentAngle(); + + const double c = std::cos(angle); + const double s = std::sin(angle); + local_body_force[0] = c * body_force[0] + s * body_force[1]; + local_body_force[1] = -s * body_force[0] + c * body_force[1]; + local_body_force[2] = 0.0; + return local_body_force; + } else { + BoundedMatrix T; + noalias(T) = GetFrenetSerretMatrix(); // global to local + noalias(local_body_force) = prod(T, body_force); + return local_body_force; + } + + return array_1d(); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +typename TotalLagrangianTrussElement::MatrixType TotalLagrangianTrussElement:: +CalculateGeometricStiffnessMatrix( + const double Stress_x, + const double Ref_Length +) +{ + KRATOS_TRY + MatrixType K_geometric(SystemSize, SystemSize); + K_geometric.clear(); + + for (IndexType i = 0; i < SystemSize; ++i) { + for (IndexType j = 0; j < SystemSize; ++j) { + if (i == j) { + K_geometric(i, i) = 1.0; + } else if (i == j + Dimension) { + K_geometric(i, j) = -1.0; + } else if (i + Dimension == j) { + K_geometric(i, j) = -1.0; + } + } + } + return K_geometric * Stress_x / std::pow(Ref_Length, 2); + + KRATOS_CATCH("CalculateGeometricStiffnessMatrix") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::CalculateLocalSystem( + MatrixType& rLHS, + VectorType& rRHS, + 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); + } + rLHS.clear(); + + if (rRHS.size() != SystemSize) { + rRHS.resize(SystemSize, false); + } + rRHS.clear(); + + const auto& r_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 ref_length = CalculateReferenceLength(); + const double curr_length = CalculateCurrentLength(); + const double J = 0.5 * ref_length; // Jacobian to the natural coordinates + const double ref_area = r_props[CROSS_AREA]; + const double Fxx = curr_length / ref_length; // Deformation gradient in the axial direction + const auto pre_stress = r_props.Has(TRUSS_PRESTRESS_PK2) ? r_props[TRUSS_PRESTRESS_PK2] : 0.0; + + // 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, N_shapeZ; + + // Loop over the integration points + for (SizeType IP = 0; IP < r_integration_points.size(); ++IP) { + const auto local_body_forces = GetLocalAxesBodyForce(*this, r_integration_points, IP); + + const double xi = r_integration_points[IP].X(); + const double weight = r_integration_points[IP].Weight(); + const double jacobian_weight = weight * J * ref_area; + GetShapeFunctionsValues(N_shape, ref_length, xi); + GetShapeFunctionsValuesY(N_shapeY, ref_length, xi); + GetShapeFunctionsValuesZ(N_shapeZ, ref_length, xi); + GetFirstDerivativesShapeFunctionsValues(B, ref_length, xi); + + strain_vector[0] = CalculateGreenLagrangeStrain(curr_length, ref_length); + + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix + + noalias(rRHS) -= Fxx * B * (stress_vector[0] + pre_stress) * jacobian_weight; + + noalias(rLHS) += Fxx * Fxx * outer_prod(B, B) * constitutive_matrix(0, 0) * jacobian_weight; // Material stiffness + noalias(rLHS) += CalculateGeometricStiffnessMatrix(stress_vector[0] + pre_stress, ref_length) * jacobian_weight; // Geometric stiffness + + noalias(rRHS) += N_shape * local_body_forces[0] * jacobian_weight; + noalias(rRHS) += N_shapeY * local_body_forces[1] * jacobian_weight; + noalias(rRHS) += N_shapeZ * local_body_forces[2] * jacobian_weight; + } + RotateAll(rLHS, rRHS); // rotate to global + + KRATOS_CATCH("CalculateLocalSystem") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::CalculateLeftHandSide( + MatrixType& rLHS, + 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); + } + rLHS.clear(); + + const auto& r_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 , false); + r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true); + + const double ref_length = CalculateReferenceLength(); + const double curr_length = CalculateCurrentLength(); + const double J = 0.5 * ref_length; // Jacobian to the natural coordinates + const double ref_area = r_props[CROSS_AREA]; + const double Fxx = curr_length / ref_length; // Deformation gradient in the axial direction + const auto pre_stress = r_props.Has(TRUSS_PRESTRESS_PK2) ? r_props[TRUSS_PRESTRESS_PK2] : 0.0; + + // 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; // Derivatives of the shape functions with respect to reference coordinates + + // Loop over the integration points + for (SizeType IP = 0; IP < r_integration_points.size(); ++IP) { + const auto local_body_forces = GetLocalAxesBodyForce(*this, r_integration_points, IP); + + const double xi = r_integration_points[IP].X(); + const double weight = r_integration_points[IP].Weight(); + const double jacobian_weight = weight * J * ref_area; + GetFirstDerivativesShapeFunctionsValues(B, ref_length, xi); + + strain_vector[0] = CalculateGreenLagrangeStrain(curr_length, ref_length); + + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix + + noalias(rLHS) += Fxx * Fxx * outer_prod(B, B) * constitutive_matrix(0, 0) * jacobian_weight; // Material stiffness + noalias(rLHS) += CalculateGeometricStiffnessMatrix(stress_vector[0] + pre_stress, ref_length) * jacobian_weight; // Geometric stiffness + + } + RotateLHS(rLHS); // rotate to global + + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::CalculateRightHandSide( + VectorType& rRHS, + const ProcessInfo& rProcessInfo + ) +{ + KRATOS_TRY; + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + + if (rRHS.size() != SystemSize) { + rRHS.resize(SystemSize, false); + } + rRHS.clear(); + + const auto& r_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, false); + + const double ref_length = CalculateReferenceLength(); + const double curr_length = CalculateCurrentLength(); + const double J = 0.5 * ref_length; + const double ref_area = r_props[CROSS_AREA]; + const double Fxx = curr_length / ref_length; + const auto pre_stress = r_props.Has(TRUSS_PRESTRESS_PK2) ? r_props[TRUSS_PRESTRESS_PK2] : 0.0; + + // 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, N_shapeZ; + + // Loop over the integration points + for (SizeType IP = 0; IP < r_integration_points.size(); ++IP) { + const auto local_body_forces = GetLocalAxesBodyForce(*this, r_integration_points, IP); + + const double xi = r_integration_points[IP].X(); + const double weight = r_integration_points[IP].Weight(); + const double jacobian_weight = weight * J * ref_area; + GetShapeFunctionsValues(N_shape, ref_length, xi); + GetShapeFunctionsValuesY(N_shapeY, ref_length, xi); + GetShapeFunctionsValuesZ(N_shapeZ, ref_length, xi); + GetFirstDerivativesShapeFunctionsValues(B, ref_length, xi); + + strain_vector[0] = CalculateGreenLagrangeStrain(curr_length, ref_length); + + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); // fills stress and const. matrix + + noalias(rRHS) -= Fxx * B * (stress_vector[0] + pre_stress) * jacobian_weight; + + noalias(rRHS) += N_shape * local_body_forces[0] * jacobian_weight; + noalias(rRHS) += N_shapeY * local_body_forces[1] * jacobian_weight; + noalias(rRHS) += N_shapeZ * local_body_forces[2] * jacobian_weight; + } + RotateRHS(rRHS); // rotate to global + + KRATOS_CATCH("CalculateRightHandSide") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::RotateLHS( + MatrixType& rLHS +) +{ + BoundedMatrix global_size_T, aux_product; + BoundedMatrix T; + + if constexpr (Dimension == 2) { + const double angle = GetCurrentAngle(); + + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + + } else { + noalias(T) = trans(GetFrenetSerretMatrix()); // global to local + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor3D2NTruss(T, global_size_T); + } + + noalias(aux_product) = prod(rLHS, trans(global_size_T)); + noalias(rLHS) = prod(global_size_T, aux_product); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::RotateRHS( + VectorType& rRHS +) +{ + BoundedMatrix T; + BoundedMatrix global_size_T; + BoundedVector local_rhs; + + if constexpr (Dimension == 2) { + const double angle = GetCurrentAngle(); + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + + } else { + noalias(T) = trans(GetFrenetSerretMatrix()); // global to local + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor3D2NTruss(T, global_size_T); + } + noalias(local_rhs) = rRHS; + noalias(rRHS) = prod(global_size_T, local_rhs); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::RotateAll( + MatrixType& rLHS, + VectorType& rRHS +) +{ + BoundedMatrix T; + BoundedMatrix global_size_T, aux_product; + BoundedVector local_rhs; + + if constexpr (Dimension == 2) { + const double angle = GetCurrentAngle(); + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(T, angle); + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor2D2NTruss(T, global_size_T); + + } else { + noalias(T) = trans(GetFrenetSerretMatrix()); // global to local + StructuralMechanicsElementUtilities::BuildElementSizeRotationMatrixFor3D2NTruss(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 TotalLagrangianTrussElement::CalculateMassMatrix( + TotalLagrangianTrussElement::MatrixType& rMassMatrix, + const ProcessInfo& rCurrentProcessInfo +) +{ + KRATOS_TRY + + const auto &r_props = GetProperties(); + + if (rMassMatrix.size1() != SystemSize || rMassMatrix.size2() != SystemSize) { + rMassMatrix.resize(SystemSize, SystemSize, false); + } + rMassMatrix.clear(); + + const double A = r_props[CROSS_AREA]; + const double L0 = CalculateReferenceLength(); + const double rho = r_props[DENSITY]; + const double total_mass = rho * A * L0; + + if (StructuralMechanicsElementUtilities::ComputeLumpedMassMatrix(r_props, rCurrentProcessInfo)) { + // Compute lumped mass matrix + noalias(rMassMatrix) = 0.5 * total_mass * IdentityMatrix(SystemSize, SystemSize); + } else { + // Compute consistent mass matrix + for (IndexType i = 0; i < SystemSize; ++i) { + for (IndexType j = 0; j < SystemSize; ++j) { + if (i == j) { + rMassMatrix(i, i) = 2.0; + } else if (i == j + Dimension) { + rMassMatrix(i, j) = 1.0; + } else if (i + Dimension == j) { + rMassMatrix(i, j) = 1.0; + } + } + } + rMassMatrix *= total_mass / 6.0; + } + + KRATOS_CATCH("CalculateMassMatrix") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::CalculateDampingMatrix( + TotalLagrangianTrussElement::MatrixType& rDampingMatrix, + const ProcessInfo& rCurrentProcessInfo +) +{ + KRATOS_TRY + StructuralMechanicsElementUtilities::CalculateRayleighDampingMatrix( + *this, + rDampingMatrix, + rCurrentProcessInfo, + SystemSize); + KRATOS_CATCH("CalculateDampingMatrix") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rProcessInfo + ) +{ + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + rOutput.resize(integration_points.size()); + + if (rVariable == AXIAL_FORCE) { + ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo); + VectorType strain_vector(1), stress_vector(1); + MatrixType C(1,1); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C); + + const double ref_length = CalculateReferenceLength(); + const double curr_length = CalculateCurrentLength(); + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); + + SystemSizeBoundedArrayType B; + const double ref_area = GetProperties()[CROSS_AREA]; + + double pre_stress = 0.0; + if (GetProperties().Has(TRUSS_PRESTRESS_PK2)) { + pre_stress = GetProperties()[TRUSS_PRESTRESS_PK2]; + } + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const double xi = integration_points[IP].X(); + GetFirstDerivativesShapeFunctionsValues(B, ref_length, xi); + + strain_vector[0] = CalculateGreenLagrangeStrain(curr_length, ref_length); + + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); + + rOutput[IP] = (cl_values.GetStressVector()[0] + pre_stress) * ref_area; + } + } else if (rVariable == AXIAL_STRAIN) { + ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo); + VectorType strain_vector(1), stress_vector(1); + MatrixType C(1,1); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, C); + + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); + + const double ref_length = CalculateReferenceLength(); + const double curr_length = CalculateCurrentLength(); + + SystemSizeBoundedArrayType B; + + // Loop over the integration points + for (SizeType IP = 0; IP < integration_points.size(); ++IP) { + const double xi = integration_points[IP].X(); + GetFirstDerivativesShapeFunctionsValues(B, ref_length, xi); + rOutput[IP] = CalculateGreenLagrangeStrain(curr_length, ref_length); + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::CalculateOnIntegrationPoints( + const Variable& rVariable, + std::vector& rOutput, + const ProcessInfo& rProcessInfo + ) +{ + const auto& integration_points = IntegrationPoints(GetIntegrationMethod()); + rOutput.resize(integration_points.size()); + + if (rVariable == PK2_STRESS_VECTOR) { + ConstitutiveLaw::Parameters cl_values(GetGeometry(), GetProperties(), rProcessInfo); + VectorType strain_vector(1), stress_vector(1); + MatrixType constitutive_matrix(1,1); + StructuralMechanicsElementUtilities::InitializeConstitutiveLawValuesForStressCalculation(cl_values, strain_vector, stress_vector, constitutive_matrix); + + const double ref_length = CalculateReferenceLength(); + const double curr_length = CalculateCurrentLength(); + + SystemSizeBoundedArrayType nodal_values(SystemSize); + GetNodalValuesVector(nodal_values); + + SystemSizeBoundedArrayType dN_dX; + + for (SizeType integration_point = 0; integration_point < integration_points.size(); ++integration_point) { + GetFirstDerivativesShapeFunctionsValues(dN_dX, ref_length, integration_points[integration_point].X()); + strain_vector[0] = CalculateGreenLagrangeStrain(curr_length, ref_length); + mConstitutiveLawVector[integration_point]->CalculateMaterialResponsePK2(cl_values); + + auto& r_stress = cl_values.GetStressVector()[0]; + if (GetProperties().Has(TRUSS_PRESTRESS_PK2)) { + r_stress += GetProperties()[TRUSS_PRESTRESS_PK2]; + } + rOutput[integration_point] = ScalarVector(1, r_stress); + } + } +} +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::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 TotalLagrangianTrussElement::Check(const ProcessInfo& rCurrentProcessInfo) const +{ + 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(""); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::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 TotalLagrangianTrussElement::load(Serializer& rSerializer) +{ + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Element); + int IntMethod; + rSerializer.load("IntegrationMethod",IntMethod); + mThisIntegrationMethod = IntegrationMethod(IntMethod); + rSerializer.load("ConstitutiveLawVector", mConstitutiveLawVector); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +Vector TotalLagrangianTrussElement::GetBaseShapeFunctions(const double xi) const +{ + Vector coord(3), N(NNodes); + coord.clear(); + coord[0] = xi; + GetGeometry().ShapeFunctionsValues(N, coord); + return N; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +typename TotalLagrangianTrussElement::MatrixType TotalLagrangianTrussElement:: +CalculateClosedFormK( + const double Sxx +) +{ + MatrixType K = ZeroMatrix(SystemSize, SystemSize); + const auto &r_props = GetProperties(); + const double E = r_props[YOUNG_MODULUS]; + const double A = r_props[CROSS_AREA]; + const double L0 = CalculateReferenceLength(); + const double L = CalculateCurrentLength(); + const double Fxx = L / L0; + const double theta = GetCurrentAngle(); + + const double k = Fxx * Fxx * E * A / L0; + + // Material stiffness matrix + K(0,0) = k * std::cos(theta) * std::cos(theta); + K(0,1) = k * std::cos(theta) * std::sin(theta); + K(1,0) = K(0,1); + K(1,1) = k * std::sin(theta) * std::sin(theta); + + K(0,2) = -k * std::cos(theta) * std::cos(theta); + K(1,2) = -k * std::sin(theta) * std::cos(theta); + K(2,0) = K(0,2); + K(2,1) = K(1,2); + K(2,2) = k * std::cos(theta) * std::cos(theta); + + K(3,3) = k * std::sin(theta) * std::sin(theta); + K(2,3) = k * std::cos(theta) * std::sin(theta); + K(3,2) = K(2,3); + K(1,3) = -k * std::sin(theta) * std::sin(theta); + K(3,1) = K(1,3); + K(0,3) = -k * std::cos(theta) * std::sin(theta); + K(3,0) = K(0,3); + + // Geometric stiffness matrix + const double k_fact_g = Sxx * A / L0; + + for (IndexType i = 0; i < SystemSize; ++i) { + for (IndexType j = 0; j < SystemSize; ++j) { + if (i == j) { + K(i, i) += k_fact_g; + } else if (i == j + Dimension) { + K(i, j) += -k_fact_g; + } else if (i + Dimension == j) { + K(i, j) += -k_fact_g; + } + } + } + + return K; +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetValuesVector(Vector& rValues, int Step) const +{ + + KRATOS_TRY + if (rValues.size() != SystemSize) { + rValues.resize(SystemSize, false); + } + + for (IndexType i = 0; i < NNodes; ++i) { + int index = i * Dimension; + const auto& r_disp = + GetGeometry()[i].FastGetSolutionStepValue(DISPLACEMENT, Step); + + rValues[index] = r_disp[0]; + rValues[index + 1] = r_disp[1]; + if constexpr (Dimension == 3) + rValues[index + 2] = r_disp[2]; + } + KRATOS_CATCH("GetValuesVector") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetFirstDerivativesVector(Vector& rValues, int Step) const +{ + + KRATOS_TRY + if (rValues.size() != SystemSize) { + rValues.resize(SystemSize, false); + } + + for (IndexType i = 0; i < NNodes; ++i) { + int index = i * Dimension; + const auto& r_vel = + GetGeometry()[i].FastGetSolutionStepValue(VELOCITY, Step); + + rValues[index] = r_vel[0]; + rValues[index + 1] = r_vel[1]; + if constexpr (Dimension == 3) + rValues[index + 2] = r_vel[2]; + } + KRATOS_CATCH("GetFirstDerivativesVector") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void TotalLagrangianTrussElement::GetSecondDerivativesVector(Vector& rValues, int Step) const +{ + + KRATOS_TRY + if (rValues.size() != SystemSize) { + rValues.resize(SystemSize, false); + } + + for (IndexType i = 0; i < NNodes; ++i) { + int index = i * Dimension; + const auto& r_acc = GetGeometry()[i].FastGetSolutionStepValue(ACCELERATION, Step); + + rValues[index] = r_acc[0]; + rValues[index + 1] = r_acc[1]; + if constexpr (Dimension == 3) + rValues[index + 2] = r_acc[2]; + } + + KRATOS_CATCH("GetSecondDerivativesVector") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template class TotalLagrangianTrussElement<2>; +template class TotalLagrangianTrussElement<3>; + +} // Namespace Kratos diff --git a/applications/StructuralMechanicsApplication/custom_elements/truss_elements/total_lagrangian_truss_element.h b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/total_lagrangian_truss_element.h new file mode 100644 index 000000000000..97efa2a54448 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/total_lagrangian_truss_element.h @@ -0,0 +1,578 @@ +// 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 TotalLagrangianTrussElement + * @ingroup StructuralMechanicsApplication + * @brief This is the Total Lagrangian TRUSS element of 2 nodes for 2D and 3D. + * O---------O --> x' + * 0 1 + * @author Alejandro Cornejo + * Reference: "Non-linear Finite Element Analysis of Solids and Structures" by Belytschko, Liu, Moran. + */ +template +class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) TotalLagrangianTrussElement + : public Element +{ + +public: + + ///@name Type Definitions + ///@{ + + /// The base element type + using BaseType = Element; + using MatrixType = typename BaseType::MatrixType; + static constexpr SizeType NNodes = 2; + static constexpr SizeType Dimension = TDimension; + static constexpr SizeType DofsPerNode = TDimension; + static constexpr SizeType SystemSize = TDimension * 2; + using SystemSizeBoundedArrayType = array_1d; + + // Counted pointer of BaseSolidElement + KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(TotalLagrangianTrussElement); + + ///@} + ///@name Life Cycle + ///@{ + + // Constructor void + TotalLagrangianTrussElement() + { + } + + // Constructor using an array of nodes + TotalLagrangianTrussElement(IndexType NewId, GeometryType::Pointer pGeometry) : Element(NewId, pGeometry) + { + mThisIntegrationMethod = CalculateIntegrationMethod(); + } + + // Constructor using an array of nodes with properties + TotalLagrangianTrussElement(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) + : Element(NewId,pGeometry,pProperties) + { + mThisIntegrationMethod = CalculateIntegrationMethod(); + } + + // Copy constructor + TotalLagrangianTrussElement(TotalLagrangianTrussElement 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 This method returns the angle of the FE axis in the reference configuration + * 2D calculations + */ + double GetReferenceAngle() const + { + return StructuralMechanicsElementUtilities::GetReferenceRotationAngle2D2NBeam(GetGeometry()); + } + + /** + * @brief This method returns the angle of the FE axis in the current configuration + * 2D calculations + */ + double GetCurrentAngle() const + { + return StructuralMechanicsElementUtilities::GetCurrentRotationAngle2D2NBeam(GetGeometry()); + } + + /** + * @brief This function builds the Frenet Serret matrix that rotates from global to local axes in the reference configuration + * T = | <- t -> | x, local + * | <- n -> | y, local + * | <- m -> | z, local + * 3D calculations + */ + BoundedMatrix GetFrenetSerretMatrix() const; + + /** + * @brief This method returns the integration method depending on the Number of Nodes + */ + IntegrationMethod CalculateIntegrationMethod() + { + return GeometryData::IntegrationMethod::GI_GAUSS_1; + } + + /** + * @brief This method returns the base shape functions in a reduced size vector + */ + Vector GetBaseShapeFunctions(const double xi) const; + + /** + * @brief Returns a n component vector including the values of the DoFs + * in LOCAL truss axes + */ + void GetNodalValuesVector(SystemSizeBoundedArrayType& rNodalValue) const; + + /** + * @brief Computes the reference (initial) length of the FE and returns it + */ + double CalculateReferenceLength() const + { + if constexpr (Dimension == 2) { + return StructuralMechanicsElementUtilities::CalculateReferenceLength2D2N(*this); + } else { + return StructuralMechanicsElementUtilities::CalculateReferenceLength3D2N(*this); + } + } + + /** + * @brief Computes the current length (updated) of the FE and returns it + */ + double CalculateCurrentLength() const + { + if constexpr (Dimension == 2) { + return StructuralMechanicsElementUtilities::CalculateCurrentLength2D2N(*this); + } else { + return StructuralMechanicsElementUtilities::CalculateCurrentLength3D2N(*this); + } + } + + /** + * @brief Computes the Green-Lagrange strain of the element + * @param CurrentLength The current length of the element + * @param ReferenceLength The reference length of the element + * @return The Green-Lagrange strain + */ + double CalculateGreenLagrangeStrain(const double CurrentLength, const double ReferenceLength) const + { + return 0.5 * (std::pow(CurrentLength / ReferenceLength, 2) - 1.0); + } + + /** + * @brief Computes the Geometric stiffness matrix of the element + * @param Stress_x The axial Cauchy stress of the element + * @param Ref_Length The reference length of the element + */ + MatrixType CalculateGeometricStiffnessMatrix( + const double Stress_x, + const double Ref_Length + ); + + MatrixType CalculateClosedFormK(const double Sxx); + + void GetValuesVector( + Vector& rValues, + int Step = 0) const override; + + void GetSecondDerivativesVector( + Vector& rValues, + int Step = 0) const override; + + void GetFirstDerivativesVector( + Vector& rValues, + int Step = 0) const override; + + /** + * @brief Called to initialize the element. + * @warning Must be called before any calculation is done + */ + void Initialize(const ProcessInfo& rCurrentProcessInfo) override; + + /** + * @brief This method is called to calculate the mass matrix + * This method discerns wether to calculate a lumped or consistent mass matrix + */ + void CalculateMassMatrix( + MatrixType& rMassMatrix, + const ProcessInfo& rCurrentProcessInfo) override; + + /** + * @brief This method calculates the damping matrix + */ + void CalculateDampingMatrix( + MatrixType& rDampingMatrix, + 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 2 shape functions used for interpolating the displacements in x,y,z + * Also the derivatives of u to compute the longitudinal strain + * @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(SystemSizeBoundedArrayType& rN, const double Length, const double xi) const; + void GetShapeFunctionsValuesY(SystemSizeBoundedArrayType& rN, const double Length, const double xi) const; + void GetShapeFunctionsValuesZ(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 + * @param rLHS the left hand side + * @param rGeometry the geometry of the FE + */ + void RotateLHS( + MatrixType &rLHS); + + /** + * @brief This function rotates the RHS from local to global coordinates + * @param rRHS the right hand side + * @param rGeometry the geometry of the FE + */ + void RotateRHS( + VectorType &rRHS); + + /** + * @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 + */ + void RotateAll( + MatrixType &rLHS, + VectorType &rRHS); + + /** + * @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 Calculate a Vector 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 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_1; /// 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 TotalLagrangianTrussElement. + +///@} +///@name Type Definitions +///@{ + + +///@} +///@name Input and output +///@{ + +} // namespace Kratos. diff --git a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp index 0e3cb0519e7e..be68c6716827 100644 --- a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp +++ b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.cpp @@ -253,13 +253,15 @@ double CalculateReferenceLength2D2N(const Element& rElement) double CalculateCurrentLength2D2N(const Element& rElement) { - KRATOS_TRY; + KRATOS_TRY + + const auto& r_geometry = rElement.GetGeometry(); const array_1d delta_pos = - rElement.GetGeometry()[1].GetInitialPosition().Coordinates() - - rElement.GetGeometry()[0].GetInitialPosition().Coordinates() + - rElement.GetGeometry()[1].FastGetSolutionStepValue(DISPLACEMENT) - - rElement.GetGeometry()[0].FastGetSolutionStepValue(DISPLACEMENT); + r_geometry[1].GetInitialPosition().Coordinates() - + r_geometry[0].GetInitialPosition().Coordinates() + + r_geometry[1].FastGetSolutionStepValue(DISPLACEMENT) - + r_geometry[0].FastGetSolutionStepValue(DISPLACEMENT); const double l = std::sqrt((delta_pos[0] * delta_pos[0]) + (delta_pos[1] * delta_pos[1])); @@ -293,13 +295,15 @@ double CalculateReferenceLength3D2N(const Element& rElement) double CalculateCurrentLength3D2N(const Element& rElement) { - KRATOS_TRY; + KRATOS_TRY + + const auto& r_geometry = rElement.GetGeometry(); const array_1d delta_pos = - rElement.GetGeometry()[1].GetInitialPosition().Coordinates() - - rElement.GetGeometry()[0].GetInitialPosition().Coordinates() + - rElement.GetGeometry()[1].FastGetSolutionStepValue(DISPLACEMENT) - - rElement.GetGeometry()[0].FastGetSolutionStepValue(DISPLACEMENT); + r_geometry[1].GetInitialPosition().Coordinates() - + r_geometry[0].GetInitialPosition().Coordinates() + + r_geometry[1].FastGetSolutionStepValue(DISPLACEMENT) - + r_geometry[0].FastGetSolutionStepValue(DISPLACEMENT); const double l = MathUtils::Norm3(delta_pos); @@ -511,11 +515,24 @@ void BuildElementSizeRotationMatrixFor3D3NTruss( double GetReferenceRotationAngle2D2NBeam(const GeometryType& rGeometry) { - const auto &r_node_1 = rGeometry[0]; - const auto &r_node_2 = rGeometry[1]; + const auto &r_coords_1 = rGeometry[0].GetInitialPosition().Coordinates(); + const auto &r_coords_2 = rGeometry[1].GetInitialPosition().Coordinates(); - const double delta_x = r_node_2.X0() - r_node_1.X0(); - const double delta_y = r_node_2.Y0() - r_node_1.Y0(); + const double delta_x = r_coords_2[0] - r_coords_1[0]; + const double delta_y = r_coords_2[1] - r_coords_1[1]; + return std::atan2(delta_y, delta_x); +} + +/***********************************************************************************/ +/***********************************************************************************/ + +double GetCurrentRotationAngle2D2NBeam(const GeometryType& rGeometry) +{ + const auto &r_coords_1 = rGeometry[0].Coordinates(); + const auto &r_coords_2 = rGeometry[1].Coordinates(); + + const double delta_x = r_coords_2[0] - r_coords_1[0]; + const double delta_y = r_coords_2[1] - r_coords_1[1]; return std::atan2(delta_y, delta_x); } @@ -621,7 +638,7 @@ void BuildElementSizeRotationMatrixFor2D3NBeam( double CalculatePhi(const Properties& rProperties, const double L, const SizeType Plane) { - const double E = rProperties[YOUNG_MODULUS]; + const double E = rProperties[YOUNG_MODULUS]; double I, G, A_s; if (Plane == 0) { @@ -667,7 +684,10 @@ void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parame /***********************************************************************************/ /***********************************************************************************/ -BoundedMatrix GetFrenetSerretMatrix3D(const GeometryType& rGeometry) +BoundedMatrix GetFrenetSerretMatrix3D( + const GeometryType& rGeometry, + const bool UseCurrentConfiguration + ) { BoundedMatrix T; T.clear(); // global to local @@ -677,7 +697,11 @@ BoundedMatrix GetFrenetSerretMatrix3D(const GeometryType& rGeometr array_1d m; // t is the axis of the truss - noalias(t) = rGeometry[1].GetInitialPosition() - rGeometry[0].GetInitialPosition(); + if (UseCurrentConfiguration) { + noalias(t) = rGeometry[1].Coordinates() - rGeometry[0].Coordinates(); + } else { + noalias(t) = rGeometry[1].GetInitialPosition() - rGeometry[0].GetInitialPosition(); + } t /= norm_2(t); n.clear(); diff --git a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h index d9b1c6587c8a..9955c2aead02 100644 --- a/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h +++ b/applications/StructuralMechanicsApplication/custom_utilities/structural_mechanics_element_utilities.h @@ -331,12 +331,17 @@ void BuildElementSizeRotationMatrixFor3D3NTruss( BoundedMatrix& rElementSizeRotationMatrix); /** - * @brief This function computes the inclination angle of a 2 noded beam + * @brief This function computes the reference inclination angle of a 2 noded beam * @param rGeometry The geometry of the beam - * It assumes 3 dofs per node: u,v,theta */ double GetReferenceRotationAngle2D2NBeam(const GeometryType &rGeometry); +/** + * @brief This function computes the current inclination angle of a 2 noded beam + * @param rGeometry The geometry of the beam + */ +double GetCurrentRotationAngle2D2NBeam(const GeometryType& rGeometry); + /** * @brief This function computes the inclination angle of a 3 noded beam * @param rGeometry The geometry of the beam @@ -361,7 +366,7 @@ void InitializeConstitutiveLawValuesForStressCalculation(ConstitutiveLaw::Parame * This is valid for STRAIGHT lines * @param rGeometry: The geometry of the line element */ -BoundedMatrix GetFrenetSerretMatrix3D(const GeometryType& rGeometry); +BoundedMatrix GetFrenetSerretMatrix3D(const GeometryType& rGeometry, const bool UseCurrentConfiguration = false); } // namespace StructuralMechanicsElementUtilities. } // namespace Kratos. diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp index 7dfdf08375f9..868f9600e4e4 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp @@ -69,6 +69,8 @@ KratosStructuralMechanicsApplication::KratosStructuralMechanicsApplication() mLinearTrussElement2D3N(0, Element::GeometryType::Pointer(new Line2D3(Element::GeometryType::PointsArrayType(3)))), mLinearTrussElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), mLinearTrussElement3D3N(0, Element::GeometryType::Pointer(new Line3D3(Element::GeometryType::PointsArrayType(3)))), + mTotalLagrangianTrussElement2D2N(0, Element::GeometryType::Pointer(new Line2D2(Element::GeometryType::PointsArrayType(2)))), + mTotalLagrangianTrussElement3D2N(0, Element::GeometryType::Pointer(new Line3D2(Element::GeometryType::PointsArrayType(2)))), // 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)))), @@ -564,6 +566,8 @@ void KratosStructuralMechanicsApplication::Register() { KRATOS_REGISTER_ELEMENT("LinearTrussElement2D3N", mLinearTrussElement2D3N) KRATOS_REGISTER_ELEMENT("LinearTrussElement3D2N", mLinearTrussElement3D2N) KRATOS_REGISTER_ELEMENT("LinearTrussElement3D3N", mLinearTrussElement3D3N) + KRATOS_REGISTER_ELEMENT("TotalLagrangianTrussElement2D2N", mTotalLagrangianTrussElement2D2N) + KRATOS_REGISTER_ELEMENT("TotalLagrangianTrussElement3D2N", mTotalLagrangianTrussElement3D2N) // 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 5a29fb7bd0e6..79c374cd9584 100644 --- a/applications/StructuralMechanicsApplication/structural_mechanics_application.h +++ b/applications/StructuralMechanicsApplication/structural_mechanics_application.h @@ -35,6 +35,7 @@ #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.h" +#include "custom_elements/truss_elements/total_lagrangian_truss_element.h" /* Adding beam element */ #include "custom_elements/beam_elements/cr_beam_element_3D2N.hpp" @@ -284,6 +285,8 @@ class KRATOS_API(STRUCTURAL_MECHANICS_APPLICATION) KratosStructuralMechanicsAppl const LinearTrussElement<2, 3> mLinearTrussElement2D3N; const LinearTrussElement<3, 2> mLinearTrussElement3D2N; const LinearTrussElement<3, 3> mLinearTrussElement3D3N; + const TotalLagrangianTrussElement<2> mTotalLagrangianTrussElement2D2N; + const TotalLagrangianTrussElement<3> mTotalLagrangianTrussElement3D2N; // Adding the beam element const CrBeamElement3D2N mCrBeamElement3D2N; diff --git a/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test.mdpa b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test.mdpa new file mode 100644 index 000000000000..10231ba888dc --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test.mdpa @@ -0,0 +1,55 @@ +Begin ModelPartData +// VARIABLE_NAME value +End ModelPartData + +Begin Properties 0 +End Properties +Begin Nodes + 1 10.1703000000 -4.8240600000 0.0000000000 + 2 -3.1328040000 0.3632200000 0.0000000000 + 3 0.1816120000 5.5732100000 0.0000000000 + 5 0.1816120000 5.5732100000 10.0000000000 +End Nodes + + +Begin Elements TotalLagrangianTrussElement3D2N// GUI group identifier: Truss Auto1 + 1 0 3 1 + 2 0 3 5 + 3 0 1 5 + 4 0 2 3 + 5 0 2 5 + 6 0 2 1 +End Elements + +Begin SubModelPart Parts_Truss_Truss_Auto1 // Group Truss Auto1 // Subtree Parts_Truss + Begin SubModelPartNodes + 1 + 2 + 3 + 5 + End SubModelPartNodes + Begin SubModelPartElements + 1 + 2 + 3 + 4 + 5 + 6 + End SubModelPartElements +End SubModelPart + +Begin SubModelPart DISPLACEMENT_Displacement_Auto1 // Group Displacement Auto1 // Subtree DISPLACEMENT + Begin SubModelPartNodes + 1 + 5 + 2 + End SubModelPartNodes +End SubModelPart +Begin SubModelPart SelfWeight3D_Self_weight_Auto1 // Group Self weight Auto1 // Subtree SelfWeight3D + Begin SubModelPartNodes + 1 + 2 + 3 + 5 + End SubModelPartNodes +End SubModelPart diff --git a/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_materials.json b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_materials.json new file mode 100644 index 000000000000..ead77e1c1200 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_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" : 206900000000.0, + "CROSS_AREA" : 1.0, + "TRUSS_PRESTRESS_PK2" : 0.0 + }, + "Tables" : null + } + }] +} diff --git a/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_parameters.json b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_parameters.json new file mode 100644 index 000000000000..d123d0eabd21 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_parameters.json @@ -0,0 +1,100 @@ +{ + "problem_data" : { + "problem_name" : "test_truss_3d", + "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" : 3, + "echo_level" : 1, + "analysis_type" : "non_linear", + "model_import_settings" : { + "input_type" : "mdpa", + "input_filename" : "TLTruss3D/tl_3d_truss_test" + }, + "material_import_settings" : { + "materials_filename" : "TLTruss3D/tl_3d_truss_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] + } + }], + "loads_process_list" : [{ + "python_module" : "assign_vector_by_direction_process", + "kratos_module" : "KratosMultiphysics", + "check" : "DirectorVectorNonZero direction", + "process_name" : "AssignVectorByDirectionProcess", + "Parameters" : { + "model_part_name" : "Structure.SelfWeight3D_Self_weight_Auto1", + "variable_name" : "VOLUME_ACCELERATION", + "interval" : [0.0,"End"], + "constrained" : false, + "modulus" : 9.81, + "direction" : [0.0,-1.0,-1.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" : "TLTruss3D/tl_3d_truss_test_results.json", + "model_part_name" : "Structure.Parts_Truss_Truss_Auto1", + "tolerance" : 1.0e-8, + "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" : "TLTruss3D/tl_3d_truss_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/TLTruss3D/tl_3d_truss_test_results.json b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_results.json new file mode 100644 index 000000000000..e0e9ce08cfe8 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/TLTruss3D/tl_3d_truss_test_results.json @@ -0,0 +1,121 @@ +{ + "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": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "NODE_3": { + "DISPLACEMENT_X": [ + 1.4689718166419459e-05 + ], + "DISPLACEMENT_Y": [ + -3.035326472476075e-05 + ], + "DISPLACEMENT_Z": [ + -4.0258073304280336e-05 + ] + }, + "NODE_5": { + "DISPLACEMENT_X": [ + 0.0 + ], + "DISPLACEMENT_Y": [ + 0.0 + ], + "DISPLACEMENT_Z": [ + 0.0 + ] + }, + "ELEMENT_1": { + "AXIAL_FORCE": { + "0": [ + -460147.6160865814 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -2.2240097442560725e-06 + ] + } + }, + "ELEMENT_2": { + "AXIAL_FORCE": { + "0": [ + 832942.3896180965 + ] + }, + "AXIAL_STRAIN": { + "0": [ + 4.025821119468809e-06 + ] + } + }, + "ELEMENT_3": { + "AXIAL_FORCE": { + "0": [ + 0.0 + ] + }, + "AXIAL_STRAIN": { + "0": [ + 0.0 + ] + } + }, + "ELEMENT_4": { + "AXIAL_FORCE": { + "0": [ + -593910.209888815 + ] + }, + "AXIAL_STRAIN": { + "0": [ + -2.870518172493064e-06 + ] + } + }, + "ELEMENT_5": { + "AXIAL_FORCE": { + "0": [ + 0.0 + ] + }, + "AXIAL_STRAIN": { + "0": [ + 0.0 + ] + } + }, + "ELEMENT_6": { + "AXIAL_FORCE": { + "0": [ + 0.0 + ] + }, + "AXIAL_STRAIN": { + "0": [ + 0.0 + ] + } + } +} \ No newline at end of file diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_structural_mechanics_element_utilities.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_structural_mechanics_element_utilities.cpp index b428f582b8c4..5f87386975fa 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_structural_mechanics_element_utilities.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_structural_mechanics_element_utilities.cpp @@ -108,9 +108,9 @@ KRATOS_TEST_CASE_IN_SUITE(CalculateElementLength, KratosStructuralMechanicsFastS model_part.CreateNewNode(3, 0.0, 0.0, 0.0); // Set the updated node-position (due to MOVE-MESH) - p_node_2->X() = 0.5; - p_node_2->Y() = 0.5; - p_node_2->Z() = 0.5; + p_node_2->X() = 2.0; + p_node_2->Y() = 2.0; + p_node_2->Z() = 2.0; // Set the DISPLACEMENT p_node_2->FastGetSolutionStepValue(DISPLACEMENT_X) = 1.0; diff --git a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp index 75c18de5afd7..80212d031bfd 100644 --- a/applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp +++ b/applications/StructuralMechanicsApplication/tests/cpp_tests/test_truss.cpp @@ -333,6 +333,7 @@ void CreateTrussModel2N_and_CheckPK2Stress(std::string TrussElementName) p_bottom_node->FastGetSolutionStepValue(DISPLACEMENT) = array_1d{0.0, 0.0, 0.0}; p_top_node->FastGetSolutionStepValue(DISPLACEMENT) = array_1d{0.0, 0.0, elongation}; + p_top_node->Coordinates()[2] += elongation; // to account for Green-Lagrange strain auto p_truss_element = dynamic_cast(p_element.get()); KRATOS_EXPECT_NE(p_truss_element, nullptr); @@ -360,6 +361,7 @@ void CreateTrussModel2N_and_CheckPK2Stress(std::string TrussElementName) p_bottom_node->FastGetSolutionStepValue(DISPLACEMENT) = array_1d{0.0, 0.0, 0.0}; p_top_node->FastGetSolutionStepValue(DISPLACEMENT) = array_1d{0.0, 0.0, elongation}; + p_top_node->Coordinates()[2] += elongation; auto p_truss_element = dynamic_cast(p_element.get()); KRATOS_EXPECT_NE(p_truss_element, nullptr); diff --git a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py index 0248f0dfc621..b7e97fe0b845 100644 --- a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py +++ b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py @@ -94,6 +94,9 @@ class LinearTruss2D3NTest(StructuralMechanicsTestFactory): class LinearTruss3DTest(StructuralMechanicsTestFactory): file_name = "LinearTruss3D/linear_3d_truss_test" +class TLTruss3DTest(StructuralMechanicsTestFactory): + file_name = "TLTruss3D/tl_3d_truss_test" + class TimoshenkoBeam3D2NTest(StructuralMechanicsTestFactory): file_name = "TimoshenkoBeams/3D2N_straight/timoshenko_beam_3d2N_test" diff --git a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py index 41a1554d0b03..faae135ed443 100644 --- a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py +++ b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py @@ -108,6 +108,7 @@ from structural_mechanics_test_factory import LinearTruss2D2NTest as TLinearTruss2D2NTest from structural_mechanics_test_factory import LinearTruss2D3NTest as TLinearTruss2D3NTest from structural_mechanics_test_factory import LinearTruss3DTest as TLinearTruss3DTest +from structural_mechanics_test_factory import TLTruss3DTest as TTLTruss3DTest from structural_mechanics_test_factory import TimoshenkoBeam3D2NTest as TTimoshenkoBeam3D2NTest from structural_mechanics_test_factory import TimoshenkoBeam2D2NTest as TTimoshenkoBeam2D2NTest from structural_mechanics_test_factory import TimoshenkoBeam2D3NTest as TTimoshenkoBeam2D3NTest @@ -369,6 +370,7 @@ def AssembleTestSuites(): smallSuite.addTest(TLinearTruss2D2NTest('test_execution')) smallSuite.addTest(TLinearTruss2D3NTest('test_execution')) smallSuite.addTest(TLinearTruss3DTest('test_execution')) + smallSuite.addTest(TTLTruss3DTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam3D2NTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam2D2NTest('test_execution')) smallSuite.addTest(TTimoshenkoBeam2D3NTest('test_execution'))