diff --git a/src/for_3D_build/particle_dynamics/slender_structure_math.cpp b/src/for_3D_build/particle_dynamics/slender_structure_math.cpp index c1aef91895..ce30e394df 100644 --- a/src/for_3D_build/particle_dynamics/slender_structure_math.cpp +++ b/src/for_3D_build/particle_dynamics/slender_structure_math.cpp @@ -6,20 +6,7 @@ namespace SPH namespace slender_structure_dynamics { //=================================================================================================// - Vec2d getVectorAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles) - { - /**The rotation matrix. */ - Real sin_angle = sin(rotation_angles[0]); - Real cos_angle = cos(rotation_angles[0]); - - Mat2d rotation_matrix{ - {cos_angle, sin_angle}, // First row - {-sin_angle,cos_angle}, //Second row - }; - - return rotation_matrix * initial_vector; - } - //=================================================================================================// + Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles) { Real sin_angle_x = sin(rotation_angles[0]); @@ -44,11 +31,7 @@ namespace SPH return rotation_matrix * initial_vector; } - //=================================================================================================// - Vec2d getVectorChangeRateAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles, const Vec2d &angular_vel) - { - return Vec2d(cos(rotation_angles[0]) * angular_vel[0], -sin(rotation_angles[0]) * angular_vel[0]); - } + //=================================================================================================// Vec3d getVectorChangeRateAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles, const Vec3d &angular_vel) { @@ -65,17 +48,7 @@ namespace SPH return Vec3d(dpseudo_n_dt_0, dpseudo_n_dt_1, dpseudo_n_dt_2); } //=================================================================================================// - Vec2d getRotationFromPseudoNormalForFiniteDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt) - { - Real cos_rotation_0 = cos(rotation[0]); - Real sin_rotation_0 = sin(rotation[0]); - - Real angle_vel_dt_0 = cos_rotation_0 * (dpseudo_n_d2t[0] + sin_rotation_0 * angular_vel[0] * angular_vel[0]) - - sin_rotation_0 * (dpseudo_n_d2t[1] + cos_rotation_0 * angular_vel[0] * angular_vel[0]); - - return Vec2d(angle_vel_dt_0, 0.0); - } - //=================================================================================================// + Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt) { Real sin_rotation_0 = sin(rotation[0]); @@ -107,11 +80,7 @@ namespace SPH return Vec3d(angle_vel_dt_0, angle_vel_dt_1, 0.0); } //=================================================================================================// - Vec2d getRotationFromPseudoNormalForSmallDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt) - { - return Vec2d(dpseudo_n_d2t[0], 0); - } - //=================================================================================================// + Vec3d getRotationFromPseudoNormalForSmallDeformation(const Vec3d &dpseudo_b_n_d2t, const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt) { return Vec3d(0.0, dpseudo_n_d2t[0], 0.0); @@ -124,11 +93,7 @@ namespace SPH } //=================================================================================================// - Vec2d getNormalFromDeformationGradientTensor(const Mat2d &F) - { - return Vec2d(-F.col(0)[1], F.col(0)[0]).normalized(); - } - //=================================================================================================// + Vec3d getNormalFromDeformationGradientTensor(const Mat3d &F) { return F.col(0).cross(F.col(1)).normalized(); @@ -205,14 +170,7 @@ namespace SPH return getWENOStateWithStencilPoints(v1, v2, v3, v4); } //=================================================================================================// - Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix) - { - Vec2d local_rotation_jump = Vec2d::Zero(); - Vec2d local_pseuodo_n_jump = transformation_matrix * pseudo_n_jump; - local_rotation_jump[0] = local_pseuodo_n_jump[0]; - return transformation_matrix.transpose() * local_rotation_jump; - } - //=================================================================================================// + Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix) { Vec3d local_rotation_jump = Vec3d::Zero(); @@ -222,13 +180,7 @@ namespace SPH return transformation_matrix.transpose() * local_rotation_jump; } //=================================================================================================// - Mat2d getCorrectedAlmansiStrain(const Mat2d ¤t_local_almansi_strain, const Real &nu_) - { - Mat2d corrected_almansi_strain = current_local_almansi_strain; - corrected_almansi_strain(1,1) = -nu_ * current_local_almansi_strain(0,0) / (1.0 - nu_); - return corrected_almansi_strain; - } - //=================================================================================================// + Mat3d getCorrectedAlmansiStrain(const Mat3d ¤t_local_almansi_strain, const Real &nu_) { Mat3d corrected_almansi_strain = current_local_almansi_strain; @@ -237,15 +189,7 @@ namespace SPH return corrected_almansi_strain; } //=================================================================================================// - Mat2d getCorrectionMatrix(const Mat2d &local_deformation_part_one) - { - Real one_over_local_deformation = 1.0 / local_deformation_part_one(0,0); - return Mat2d{ - {one_over_local_deformation, 0}, - {0,0}, - }; - } - //=================================================================================================// + Mat3d getCorrectionMatrix(const Mat3d &local_deformation_part_one) { Mat3d correction_matrix = Mat3d::Zero(); diff --git a/src/for_3D_build/particle_dynamics/slender_structure_math.h b/src/for_3D_build/particle_dynamics/slender_structure_math.h index 9e65aaeda6..41305c6c72 100644 --- a/src/for_3D_build/particle_dynamics/slender_structure_math.h +++ b/src/for_3D_build/particle_dynamics/slender_structure_math.h @@ -44,24 +44,19 @@ namespace SPH * when the axis about which they occur points toward the observer, * and the coordinate system is right-handed. */ - Vec2d getVectorAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles); Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles); /** Vector change rate after rotation. */ - Vec2d getVectorChangeRateAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles, const Vec2d &angular_vel); Vec3d getVectorChangeRateAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles, const Vec3d &angular_vel); /** get the rotation from pseudo-normal for finite deformation. */ - Vec2d getRotationFromPseudoNormalForFiniteDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt); Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt); /** get the rotation from pseudo-normal for small deformation. */ - Vec2d getRotationFromPseudoNormalForSmallDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt); Vec3d getRotationFromPseudoNormalForSmallDeformation(const Vec3d &dpseudo_b_n_d2t, const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt); Vec3d getRotationFromPseudoNormalForSmallDeformation_b(const Vec3d &dpseudo_b_n_d2t, const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt); /** get the current normal direction from deformation gradient tensor. */ - Vec2d getNormalFromDeformationGradientTensor(const Mat2d &F); Vec3d getNormalFromDeformationGradientTensor(const Mat3d &F); Vec3d getBinormalFromDeformationGradientTensor(const Mat3d &F); @@ -79,15 +74,12 @@ namespace SPH const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value); /** get the artificial rotation from the pseudo-normal jump. */ - Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix); Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix); /** get the corrected Eulerian Almansi strain tensor according to plane stress problem. */ - Mat2d getCorrectedAlmansiStrain(const Mat2d ¤t_local_almansi_strain, const Real &nu_); Mat3d getCorrectedAlmansiStrain(const Mat3d ¤t_local_almansi_strain, const Real &nu_); /** get the correction matrix. */ - Mat2d getCorrectionMatrix(const Mat2d &local_deformation_part_one); Mat3d getCorrectionMatrix(const Mat3d &local_deformation_part_one); Mat3d getCorrectionMatrix_beam(const Mat3d &local_deformation_part_one); } diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp index 6c0433af8f..130aa4e791 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp @@ -5,7 +5,20 @@ namespace SPH //=====================================================================================================// namespace thin_structure_dynamics { +//=================================================================================================// +Vec2d getVectorAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles) +{ + /**The rotation matrix. */ + Real sin_angle = sin(rotation_angles[0]); + Real cos_angle = cos(rotation_angles[0]); + Mat2d rotation_matrix{ + {cos_angle, sin_angle}, // First row + {-sin_angle, cos_angle}, // Second row + }; + + return rotation_matrix * initial_vector; +} //=================================================================================================// Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles) { @@ -29,7 +42,11 @@ Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec return rotation_matrix * initial_vector; } - +//=================================================================================================// +Vec2d getVectorChangeRateAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles, const Vec2d &angular_vel) +{ + return Vec2d(cos(rotation_angles[0]) * angular_vel[0], -sin(rotation_angles[0]) * angular_vel[0]); +} //=================================================================================================// Vec3d getVectorChangeRateAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles, const Vec3d &angular_vel) { @@ -46,6 +63,16 @@ Vec3d getVectorChangeRateAfterThinStructureRotation(const Vec3d &initial_vector, return Vec3d(dpseudo_n_dt_0, dpseudo_n_dt_1, dpseudo_n_dt_2); } //=================================================================================================// +Vec2d getRotationFromPseudoNormalForFiniteDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt) +{ + Real cos_rotation_0 = cos(rotation[0]); + Real sin_rotation_0 = sin(rotation[0]); + + Real angle_vel_dt_0 = cos_rotation_0 * (dpseudo_n_d2t[0] + sin_rotation_0 * angular_vel[0] * angular_vel[0]) - sin_rotation_0 * (dpseudo_n_d2t[1] + cos_rotation_0 * angular_vel[0] * angular_vel[0]); + + return Vec2d(angle_vel_dt_0, 0.0); +} +//=================================================================================================// Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt) { Real sin_rotation_0 = sin(rotation[0]); @@ -65,13 +92,20 @@ Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t return Vec3d(angle_vel_dt_0, angle_vel_dt_1, 0.0); } //=================================================================================================// - +Vec2d getRotationFromPseudoNormalForSmallDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt) +{ + return Vec2d(dpseudo_n_d2t[0], 0); +} //=================================================================================================// Vec3d getRotationFromPseudoNormalForSmallDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt) { return Vec3d(-dpseudo_n_d2t[1], dpseudo_n_d2t[0], 0.0); } - +//=================================================================================================// +Vec2d getNormalFromDeformationGradientTensor(const Mat2d &F) +{ + return Vec2d(-F.col(0)[1], F.col(0)[0]).normalized(); +} //=================================================================================================// Vec3d getNormalFromDeformationGradientTensor(const Mat3d &F) { @@ -155,7 +189,13 @@ Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_ma local_rotation_jump[1] = local_pseuodo_n_jump[1]; return transformation_matrix.transpose() * local_rotation_jump; } - +//=================================================================================================// +Mat2d getCorrectedAlmansiStrain(const Mat2d ¤t_local_almansi_strain, const Real &nu_) +{ + Mat2d corrected_almansi_strain = current_local_almansi_strain; + corrected_almansi_strain(1, 1) = -nu_ * current_local_almansi_strain(0, 0) / (1.0 - nu_); + return corrected_almansi_strain; +} //=================================================================================================// Mat3d getCorrectedAlmansiStrain(const Mat3d ¤t_local_almansi_strain, const Real &nu_) { @@ -163,7 +203,15 @@ Mat3d getCorrectedAlmansiStrain(const Mat3d ¤t_local_almansi_strain, const corrected_almansi_strain(2, 2) = -nu_ * (current_local_almansi_strain(0, 0) + current_local_almansi_strain(1, 1)) / (1.0 - nu_); return corrected_almansi_strain; } - +//=================================================================================================// +Mat2d getCorrectionMatrix(const Mat2d &local_deformation_part_one) +{ + Real one_over_local_deformation = 1.0 / local_deformation_part_one(0, 0); + return Mat2d{ + {one_over_local_deformation, 0}, + {0, 0}, + }; +} //=================================================================================================// Mat3d getCorrectionMatrix(const Mat3d &local_deformation_part_one) { diff --git a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h index 309b4f3233..81fe69b861 100644 --- a/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h +++ b/src/shared/particle_dynamics/solid_dynamics/thin_structure_math.h @@ -44,18 +44,23 @@ namespace thin_structure_dynamics * when the axis about which they occur points toward the observer, * and the coordinate system is right-handed. */ +Vec2d getVectorAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles); Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles); /** Vector change rate after rotation. */ +Vec2d getVectorChangeRateAfterThinStructureRotation(const Vec2d &initial_vector, const Vec2d &rotation_angles, const Vec2d &angular_vel); Vec3d getVectorChangeRateAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles, const Vec3d &angular_vel); /** get the rotation from pseudo-normal for finite deformation. */ +Vec2d getRotationFromPseudoNormalForFiniteDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt); Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt); /** get the rotation from pseudo-normal for small deformation. */ +Vec2d getRotationFromPseudoNormalForSmallDeformation(const Vec2d &dpseudo_n_d2t, const Vec2d &rotation, const Vec2d &angular_vel, Real dt); Vec3d getRotationFromPseudoNormalForSmallDeformation(const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt); /** get the current normal direction from deformation gradient tensor. */ +Vec2d getNormalFromDeformationGradientTensor(const Mat2d &F); Vec3d getNormalFromDeformationGradientTensor(const Mat3d &F); /** get variable jump form gradient tensor. */ @@ -72,12 +77,15 @@ Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_ const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value); /** get the artificial rotation from the pseudo-normal jump. */ +Vec2d getRotationJump(const Vec2d &pseudo_n_jump, const Mat2d &transformation_matrix); Vec3d getRotationJump(const Vec3d &pseudo_n_jump, const Mat3d &transformation_matrix); /** get the corrected Eulerian Almansi strain tensor according to plane stress problem. */ +Mat2d getCorrectedAlmansiStrain(const Mat2d ¤t_local_almansi_strain, const Real &nu_); Mat3d getCorrectedAlmansiStrain(const Mat3d ¤t_local_almansi_strain, const Real &nu_); /** get the correction matrix. */ +Mat2d getCorrectionMatrix(const Mat2d &local_deformation_part_one); Mat3d getCorrectionMatrix(const Mat3d &local_deformation_part_one); } // namespace thin_structure_dynamics } // namespace SPH