add getRotationFromPseudoNormalForSmallDeformation.

src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp

@@ -231,8 +231,9 @@ void BarStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)
 void BarStressRelaxationFirstHalf::update(size_t index_i, Real dt)
 {
     vel_[index_i] += (acc_prior_[index_i] + acc_[index_i]) * dt;
-    angular_vel_[index_i] += (dangular_vel_dt_[index_i]) * dt;
-    angular_b_vel_[index_i] += (dangular_b_vel_dt_[index_i]) * dt;
+    Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i], pseudo_b_n_[index_i]);
+    angular_vel_[index_i] += (transformation_matrix_[index_i] * (current_transformation_matrix.transpose() * dangular_vel_dt_[index_i])) * dt;
+    angular_b_vel_[index_i] += (transformation_matrix_[index_i] * (current_transformation_matrix.transpose() * dangular_b_vel_dt_[index_i])) * dt;
 }
 //=================================================================================================//
 void BarStressRelaxationSecondHalf::initialization(size_t index_i, Real dt)

src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h

@@ -232,9 +232,9 @@ class BarStressRelaxationFirstHalf : public BaseBarRelaxation
         acc_[index_i] = acceleration * inv_rho0_ / (thickness_[index_i] * width_[index_i]);
         dpseudo_n_d2t_[index_i] = pseudo_normal_acceleration * inv_rho0_ * 12.0 / pow(thickness_[index_i], 4);
         dpseudo_b_n_d2t_[index_i] = -pseudo_b_normal_acceleration * inv_rho0_ * 12.0 / pow(thickness_[index_i], 4);
-
-        Vecd local_dpseudo_n_d2t = transformation_matrix_[index_i] * dpseudo_n_d2t_[index_i];
-        Vecd local_dpseudo_b_n_d2t = transformation_matrix_[index_i] * dpseudo_b_n_d2t_[index_i];
+        Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i], pseudo_b_n_[index_i]);
+        Vecd local_dpseudo_n_d2t = current_transformation_matrix * dpseudo_n_d2t_[index_i];
+        Vecd local_dpseudo_b_n_d2t = current_transformation_matrix * dpseudo_b_n_d2t_[index_i];
         dangular_b_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation_b(
             Vec3d(local_dpseudo_b_n_d2t), Vec3d(local_dpseudo_n_d2t), Vec3d(rotation_b_[index_i]), Vec3d(angular_b_vel_[index_i]), dt);
         dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation(

src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp

@@ -77,14 +77,14 @@ Vec3d getRotationFromPseudoNormalForFiniteDeformation(const Vec3d &dpseudo_n_d2t
 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);
+    return Vec3d(-dpseudo_n_d2t[1] - dpseudo_b_n_d2t[2], dpseudo_n_d2t[0], dpseudo_b_n_d2t[0]);
 }
 //=================================================================================================//
 
 Vec3d getRotationFromPseudoNormalForSmallDeformation_b(
     const Vec3d &dpseudo_b_n_d2t, const Vec3d &dpseudo_n_d2t, const Vec3d &rotation, const Vec3d &angular_vel, Real dt)
 {
-    return Vec3d(0.0, 0.0, dpseudo_b_n_d2t[0]);
+    return Vec3d(-dpseudo_n_d2t[1] - dpseudo_b_n_d2t[2], dpseudo_n_d2t[0], dpseudo_b_n_d2t[0]);
 }
 //=================================================================================================//
 

src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp

@@ -177,16 +177,26 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)
 void ShellStressRelaxationFirstHalf::update(size_t index_i, Real dt)
 {
     vel_[index_i] += (acc_prior_[index_i] + acc_[index_i]) * dt;
-    angular_vel_[index_i] += dangular_vel_dt_[index_i] * dt;
+    Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]);
+    angular_vel_[index_i] += (transformation_matrix_[index_i] * (current_transformation_matrix.transpose() * dangular_vel_dt_[index_i])) * dt;
 }
 //=================================================================================================//
 void ShellStressRelaxationSecondHalf::initialization(size_t index_i, Real dt)
 {
     pos_[index_i] += vel_[index_i] * dt * 0.5;
     rotation_[index_i] += angular_vel_[index_i] * dt * 0.5;
-    dpseudo_n_dt_[index_i] = transformation_matrix_[index_i].transpose() *
-                             getVectorChangeRateAfterThinStructureRotation(local_pseudo_n_0, rotation_[index_i], angular_vel_[index_i]);
-    pseudo_n_[index_i] += dpseudo_n_dt_[index_i] * dt * 0.5;
+    Vecd pseudo_n_temp = pseudo_n_[index_i];
+
+    pseudo_n_[index_i] = transformation_matrix_[index_i].transpose() * getVectorAfterThinStructureRotation(local_pseudo_n_0, rotation_[index_i]);
+
+    if (dt < 1e-10)
+    {
+        dpseudo_n_dt_[index_i] = Vecd::Zero();
+    }
+    else
+    {
+        dpseudo_n_dt_[index_i] = (pseudo_n_[index_i] - pseudo_n_temp) / (0.5 * dt);
+    }
 }
 //=================================================================================================//
 void ShellStressRelaxationSecondHalf::update(size_t index_i, Real dt)

src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h

@@ -230,8 +230,9 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation
         dpseudo_n_d2t_[index_i] = pseudo_normal_acceleration * inv_rho0_ * 12.0 / pow(thickness_[index_i], 3);
 
         /** the relation between pseudo-normal and rotations */
-        Vecd local_dpseudo_n_d2t = transformation_matrix_[index_i] * dpseudo_n_d2t_[index_i];
-        dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForFiniteDeformation(local_dpseudo_n_d2t, rotation_[index_i], angular_vel_[index_i], dt);
+        Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]);
+        Vecd local_dpseudo_n_d2t = current_transformation_matrix * dpseudo_n_d2t_[index_i];
+        dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation(local_dpseudo_n_d2t, rotation_[index_i], angular_vel_[index_i], dt);
     };
 
     void update(size_t index_i, Real dt = 0.0);

src/shared/particle_dynamics/solid_dynamics/thin_structure_math.cpp

@@ -23,22 +23,17 @@ Vec2d getVectorAfterThinStructureRotation(const Vec2d &initial_vector, const Vec
 Vec3d getVectorAfterThinStructureRotation(const Vec3d &initial_vector, const Vec3d &rotation_angles)
 {
     /**The rotation matrix is the rotation around Y-axis \times rotation around X-axis. */
-    Real sin_angle_x = sin(rotation_angles[0]);
-    Real cos_angle_x = cos(rotation_angles[0]);
-
-    Real sin_angle_y = sin(rotation_angles[1]);
-    Real cos_angle_y = cos(rotation_angles[1]);
-
-    Mat3d rotation_matrix = Mat3d{
-                                {cos_angle_y, 0.0, sin_angle_y},
-                                {0.0, 1.0, 0.0},
-                                {-sin_angle_y, 0.0, cos_angle_y},
-                            } *
-                            Mat3d{
-                                {1.0, 0.0, 0.0},
-                                {0.0, cos_angle_x, -sin_angle_x},
-                                {0.0, sin_angle_x, cos_angle_x},
-                            };
+
+    Real theta = sqrt(rotation_angles[0] * rotation_angles[0] + rotation_angles[1] * rotation_angles[1] + rotation_angles[2] * rotation_angles[2]);
+    Mat3d R1 = Mat3d::Zero();
+    R1(0, 1) = -rotation_angles[2];
+    R1(0, 2) = rotation_angles[1];
+    R1(1, 0) = rotation_angles[2];
+    R1(1, 2) = -rotation_angles[0];
+    R1(2, 0) = -rotation_angles[1];
+    R1(2, 1) = rotation_angles[0];
+
+    Mat3d rotation_matrix = Mat3d::Identity() + sin(theta) / (theta + Eps) * R1 + (1 - cos(theta)) / (theta * theta + Eps) * R1 * R1;
 
     return rotation_matrix * initial_vector;
 }