diff --git a/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp b/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp index 595172ec9b74..8cce462edcad 100644 --- a/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp +++ b/applications/StructuralMechanicsApplication/custom_constitutive/truss_constitutive_law.cpp @@ -156,6 +156,12 @@ double TrussConstitutiveLaw::CalculateStressElastic( CalculateValue(rParameterValues,TANGENT_MODULUS,tangent_modulus); const double current_stress = tangent_modulus*current_strain[0]; + + 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; } 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..9ac65a72ac9e --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.cpp @@ -0,0 +1,710 @@ +// 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 = 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 LinearTrussElement2D::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 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(); + 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(); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +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(); + 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); + } + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::GetShapeFunctionsValues( + 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); + } else { // 3N + rN[0] = 0.5 * xi * (xi - 1.0); + rN[2] = 0.5 * xi * (xi + 1.0); + rN[4] = (1.0 - std::pow(xi, 2)); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +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] = 0.5 * xi * (xi + 1.0); + rN[5] = (1.0 - std::pow(xi, 2)); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::GetFirstDerivativesShapeFunctionsValues( + 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; + rdN_dX[2] = inverse_l; + } else { // 3N + rdN_dX[0] = xi - 0.5; + rdN_dX[2] = xi + 0.5; + rdN_dX[4] = -2.0 * xi; + rdN_dX *= 2.0 / Length; // The Jacobian + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +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; + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +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(); + + 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); // 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]->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; + + noalias(rRHS) += N_shape * local_body_forces[0] * jacobian_weight; + noalias(rRHS) += N_shapeY * local_body_forces[1] * jacobian_weight; + } + RotateAll(rLHS, rRHS); // rotate to global + + KRATOS_CATCH("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::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); + } + 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]->CalculateMaterialResponsePK2(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("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::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); + } + 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]->CalculateMaterialResponsePK2(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("") +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::RotateLHS( + MatrixType& rLHS +) +{ + 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); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::RotateRHS( + VectorType& rRHS +) +{ + 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::BuildRotationMatrixForTruss(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); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::RotateAll( + MatrixType& rLHS, + VectorType& rRHS +) +{ + const double angle = GetAngle(); + if (std::abs(angle) > std::numeric_limits::epsilon()) { + BoundedMatrix T; + BoundedMatrix global_size_T, aux_product; + BoundedVector local_rhs; + StructuralMechanicsElementUtilities::BuildRotationMatrixForTruss(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); + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +template +void LinearTrussElement2D::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) { + 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, true); + + const double length = CalculateLength(); + + // 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); + + 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, length, xi); + + strain_vector[0] = inner_prod(B, nodal_values); + + mConstitutiveLawVector[IP]->CalculateMaterialResponsePK2(cl_values); + rOutput[IP] = cl_values.GetStressVector()[0] * area; + } + } else if (rVariable == AXIAL_STRAIN) { + const auto &r_props = GetProperties(); + const auto &r_geometry = GetGeometry(); + + ConstitutiveLaw::Parameters cl_values(r_geometry, r_props, rProcessInfo); + auto &r_cl_options = cl_values.GetOptions(); + r_cl_options.Set(ConstitutiveLaw::COMPUTE_STRESS , true); + r_cl_options.Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR, true); + + const double length = CalculateLength(); + + // 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); + + 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, length, xi); + rOutput[IP] = inner_prod(B, nodal_values); + } + } +} + +/***********************************************************************************/ +/***********************************************************************************/ + +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_ERROR_IF_NOT(GetProperties().Has(CROSS_AREA)) << "CROSS_AREA not defined in the properties" << std::endl; + + 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..40b7503e1119 --- /dev/null +++ b/applications/StructuralMechanicsApplication/custom_elements/truss_elements/linear_truss_element_2D.h @@ -0,0 +1,479 @@ +// 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 and 3 nodes. + * @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 NNodes = TNNodes; + static constexpr SizeType DofsPerNode = 2; + static constexpr SizeType SystemSize = DofsPerNode * NNodes; + 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) + { + mThisIntegrationMethod = CalculateIntegrationMethod(); + } + + // Constructor using an array of nodes with properties + LinearTrussElement2D(IndexType NewId, GeometryType::Pointer pGeometry, PropertiesType::Pointer pProperties) + : Element(NewId,pGeometry,pProperties) + { + mThisIntegrationMethod = CalculateIntegrationMethod(); + } + + // 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 This method returns the angle of the FE axis + */ + double GetAngle() const + { + 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 + */ + void GetNodalValuesVector(SystemSizeBoundedArrayType& rNodalValue) const; + + /** + * @brief Computes the length of the FE and returns it + */ + double CalculateLength() const + { + // Same implementation for 2N and 3N + 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(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 + * @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 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_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 diff --git a/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp b/applications/StructuralMechanicsApplication/structural_mechanics_application.cpp index 025b84f60eb2..e5d2772a256b 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)))), @@ -529,6 +531,8 @@ void KratosStructuralMechanicsApplication::Register() { 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) // 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 e677969efa29..ff17c6688238 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" @@ -271,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; @@ -420,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; 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..39990a9fe310 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D2N/linear_truss_2d2N_test_parameters.json @@ -0,0 +1,109 @@ +{ + "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" : "LinearTruss2D/2D2N/linear_truss_2d2N_test" + }, + "material_import_settings" : { + "materials_filename" : "LinearTruss2D/2D2N/linear_truss_2d2N_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.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/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" : [], + "vtk_output" : [] + }, + "analysis_stage" : "KratosMultiphysics.StructuralMechanicsApplication.structural_mechanics_analysis" +} 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 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..938cd7af45d7 --- /dev/null +++ b/applications/StructuralMechanicsApplication/tests/LinearTruss2D/2D3N/linear_truss_2d3N_test.mdpa @@ -0,0 +1,84 @@ +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..df81081fd538 --- /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/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/structural_mechanics_test_factory.py b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py index 2dc10830087d..e38c1a15ceea 100644 --- a/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py +++ b/applications/StructuralMechanicsApplication/tests/structural_mechanics_test_factory.py @@ -76,6 +76,12 @@ def tearDown(self): with KratosUnittest.WorkFolderScope(".", __file__): self.test.Finalize() +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 870aa1001fb8..e29dc9312a48 100644 --- a/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py +++ b/applications/StructuralMechanicsApplication/tests/test_StructuralMechanicsApplication.py @@ -96,6 +96,8 @@ ##### 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 @@ -344,6 +346,8 @@ 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'))