Skip to content

Commit

Permalink
Merge pull request #604 from DongWuTUM/shell_code_optimization
Browse files Browse the repository at this point in the history 
Shell code simplification
  • Loading branch information
@Xiangyu-Hu
Xiangyu-Hu authored Jun 27, 2024
2 parents 64956ed + 1494786 commit 7900ebf
Show file tree
Hide file tree
Showing 2 changed files with 18 additions and 19 deletions.
25 changes: 14 additions & 11 deletions src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,8 @@ ShellStressRelaxationFirstHalf::
global_moment_(*particles_->registerSharedVariable<Matd>("GlobalMoment")),
mid_surface_cauchy_stress_(*particles_->registerSharedVariable<Matd>("MidSurfaceCauchyStress")),
global_shear_stress_(*particles_->registerSharedVariable<Vecd>("GlobalShearStress")),
global_F_(*particles_->registerSharedVariable<Matd>("GlobalDeformationGradient")),
global_F_bending_(*particles_->registerSharedVariable<Matd>("GlobalBendingDeformationGradient")),
E0_(elastic_solid_.YoungsModulus()),
G0_(elastic_solid_.ShearModulus()),
nu_(elastic_solid_.PoissonRatio()),
Expand Down Expand Up @@ -144,13 +146,18 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)
F_[index_i] += dF_dt_[index_i] * dt * 0.5;
F_bending_[index_i] += dF_bending_dt_[index_i] * dt * 0.5;

global_F_[index_i] = transformation_matrix0_[index_i].transpose() * F_[index_i] * transformation_matrix0_[index_i];
global_F_bending_[index_i] = transformation_matrix0_[index_i].transpose() * F_bending_[index_i] * transformation_matrix0_[index_i];

Real J = F_[index_i].determinant();
Matd inverse_F = F_[index_i].inverse();
Matd inverse_transpose_global_F = global_F_[index_i].inverse().transpose();

rho_[index_i] = rho0_ / J;

/** Get transformation matrix from global coordinates to current local coordinates. */
Matd current_transformation_matrix = getTransformationMatrix(pseudo_n_[index_i]);
/** Get transformation matrix from initial local coordinates to current local coordinates. */
Matd transformation_matrix_0_to_current = current_transformation_matrix * transformation_matrix0_[index_i].transpose();

Matd resultant_stress = Matd::Zero();
Matd resultant_moment = Matd::Zero();
Expand All @@ -161,19 +168,19 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)
Matd F_gaussian_point = F_[index_i] + gaussian_point_[i] * F_bending_[index_i] * thickness_[index_i] * 0.5;
Matd dF_gaussian_point_dt = dF_dt_[index_i] + gaussian_point_[i] * dF_bending_dt_[index_i] * thickness_[index_i] * 0.5;
Matd inverse_F_gaussian_point = F_gaussian_point.inverse();
Matd current_local_almansi_strain = current_transformation_matrix * transformation_matrix0_[index_i].transpose() * 0.5 *
Matd current_local_almansi_strain = transformation_matrix_0_to_current * 0.5 *
(Matd::Identity() - inverse_F_gaussian_point.transpose() * inverse_F_gaussian_point) *
transformation_matrix0_[index_i] * current_transformation_matrix.transpose();
transformation_matrix_0_to_current.transpose();

/** correct Almansi strain tensor according to plane stress problem. */
current_local_almansi_strain = getCorrectedAlmansiStrain(current_local_almansi_strain, nu_);

/** correct out-plane numerical damping. */
numerical_damping_scaling_matrix_(Dimensions - 1, Dimensions - 1) = thickness_[i] < smoothing_length_ ? thickness_[i] : smoothing_length_;
Matd cauchy_stress = elastic_solid_.StressCauchy(current_local_almansi_strain, index_i) +
current_transformation_matrix * transformation_matrix0_[index_i].transpose() * F_gaussian_point *
transformation_matrix_0_to_current * F_gaussian_point *
elastic_solid_.NumericalDampingRightCauchy(F_gaussian_point, dF_gaussian_point_dt, numerical_damping_scaling_matrix_, index_i) *
F_gaussian_point.transpose() * transformation_matrix0_[index_i] * current_transformation_matrix.transpose() / F_gaussian_point.determinant();
F_gaussian_point.transpose() * transformation_matrix_0_to_current.transpose() / F_gaussian_point.determinant();

/** Impose modeling assumptions. */
cauchy_stress.col(Dimensions - 1) *= shear_correction_factor_;
Expand All @@ -199,13 +206,9 @@ void ShellStressRelaxationFirstHalf::initialization(size_t index_i, Real dt)

/** stress and moment in global coordinates for pair interaction */
global_stress_[index_i] = J * current_transformation_matrix.transpose() *
resultant_stress * current_transformation_matrix *
transformation_matrix0_[index_i].transpose() *
inverse_F.transpose() * transformation_matrix0_[index_i];
resultant_stress * current_transformation_matrix * inverse_transpose_global_F;
global_moment_[index_i] = J * current_transformation_matrix.transpose() *
resultant_moment * current_transformation_matrix *
transformation_matrix0_[index_i].transpose() * inverse_F.transpose() *
transformation_matrix0_[index_i];
resultant_moment * current_transformation_matrix * inverse_transpose_global_F;
global_shear_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_shear_stress;
}
//=================================================================================================//
Expand Down
12 changes: 4 additions & 8 deletions src/shared/particle_dynamics/solid_dynamics/thin_structure_dynamics.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -205,20 +205,15 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation
Vecd e_ij = inner_neighborhood.e_ij_[n];
Real r_ij = inner_neighborhood.r_ij_[n];
Real weight = inner_neighborhood.W_ij_[n] * inv_W0_;
Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i],
transformation_matrix0_[index_i].transpose() * F_[index_i] * transformation_matrix0_[index_i],
pos_[index_j],
transformation_matrix0_[index_i].transpose() * F_[index_j] * transformation_matrix0_[index_i]);
Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i], global_F_[index_i], pos_[index_j], global_F_[index_j]);
Real limiter_pos = SMIN(2.0 * pos_jump.norm() / r_ij, 1.0);
force += mass_[index_i] * hourglass_control_factor_ * weight * G0_ * pos_jump * Dimensions *
inner_neighborhood.dW_ij_[n] * Vol_[index_j] * limiter_pos;

Vecd pseudo_n_variation_i = pseudo_n_[index_i] - n0_[index_i];
Vecd pseudo_n_variation_j = pseudo_n_[index_j] - n0_[index_j];
Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_variation_i,
transformation_matrix0_[index_i].transpose() * F_bending_[index_i] * transformation_matrix0_[index_i],
pseudo_n_variation_j,
transformation_matrix0_[index_j].transpose() * F_bending_[index_j] * transformation_matrix0_[index_j]);
Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_variation_i, global_F_bending_[index_i],
pseudo_n_variation_j, global_F_bending_[index_j]);
Real limiter_pseudo_n = SMIN(2.0 * pseudo_n_jump.norm() / ((pseudo_n_variation_i - pseudo_n_variation_j).norm() + Eps), 1.0);
pseudo_normal_acceleration += hourglass_control_factor_ * weight * G0_ * pseudo_n_jump * Dimensions *
inner_neighborhood.dW_ij_[n] * Vol_[index_j] * pow(thickness_[index_i], 2) * limiter_pseudo_n;
Expand Down Expand Up @@ -246,6 +241,7 @@ class ShellStressRelaxationFirstHalf : public BaseShellRelaxation
StdLargeVec<Real> &rho_, &mass_;
StdLargeVec<Matd> &global_stress_, &global_moment_, &mid_surface_cauchy_stress_;
StdLargeVec<Vecd> &global_shear_stress_;
StdLargeVec<Matd> &global_F_, &global_F_bending_;
Real E0_, G0_, nu_, hourglass_control_factor_;
bool hourglass_control_;
const Real inv_W0_ = 1.0 / sph_body_.sph_adaptation_->getKernel()->W0(ZeroVecd);
Expand Down

0 comments on commit 7900ebf

Please sign in to comment.