diff --git a/src/for_2D_build/common/data_type.h b/src/for_2D_build/common/data_type.h index 3486dd8a21..74bea8c271 100644 --- a/src/for_2D_build/common/data_type.h +++ b/src/for_2D_build/common/data_type.h @@ -60,6 +60,7 @@ const Matd reduced_unit_matrix{ /** initial local normal, only works for thin structure dynamics. */ const Vecd local_pseudo_n_0 = Vecd(0.0, 1.0); +const Vecd local_pseudo_b_n_0 = Vecd(0.0, 1.0); const Vecd ZeroVecd = Vec2d::Zero(); } // namespace SPH diff --git a/src/for_3D_build/common/data_type.h b/src/for_3D_build/common/data_type.h index a9e9278784..f417adf4b6 100644 --- a/src/for_3D_build/common/data_type.h +++ b/src/for_3D_build/common/data_type.h @@ -59,6 +59,7 @@ const Matd reduced_unit_matrix{ }; /** initial local normal, only works for thin structure dynamics. */ const Vecd local_pseudo_n_0 = Vecd(0.0, 0.0, 1.0); +const Vecd local_pseudo_b_n_0 = Vecd(0.0, 1.0, 0.0); const Vecd ZeroVecd = Vec3d::Zero(); } // namespace SPH diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp new file mode 100644 index 0000000000..b3df7d3f07 --- /dev/null +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.cpp @@ -0,0 +1,358 @@ +#include "slender_structure_dynamics.h" + +namespace SPH +{ + //=====================================================================================================// + namespace slender_structure_dynamics + { + //=================================================================================================// + BarDynamicsInitialCondition::BarDynamicsInitialCondition(SPHBody &sph_body) + : LocalDynamics(sph_body), BarDataSimple(sph_body), + n0_(particles_->n0_), n_(particles_->n_), pseudo_n_(particles_->pseudo_n_), + pos0_(particles_->pos0_), transformation_matrix_(particles_->transformation_matrix_), + b_n0_(particles_->b_n0_), b_n_(particles_->b_n_), pseudo_b_n_(particles_->pseudo_b_n_) {} + //=================================================================================================// + BarAcousticTimeStepSize::BarAcousticTimeStepSize(SPHBody &sph_body, Real CFL) + : LocalDynamicsReduce(sph_body, Real(MaxRealNumber)), + BarDataSimple(sph_body), CFL_(CFL), + vel_(particles_->vel_), acc_(particles_->acc_), + angular_vel_(particles_->angular_vel_), dangular_vel_dt_(particles_->dangular_vel_dt_), + angular_b_vel_(particles_->angular_b_vel_), dangular_b_vel_dt_(particles_->dangular_b_vel_dt_), + acc_prior_(particles_->acc_prior_), + thickness_(particles_->thickness_), + width_(particles_->thickness_), + rho0_(particles_->elastic_solid_.ReferenceDensity()), + E0_(particles_->elastic_solid_.YoungsModulus()), + nu_(particles_->elastic_solid_.PoissonRatio()), + c0_(particles_->elastic_solid_.ReferenceSoundSpeed()), + smoothing_length_(sph_body.sph_adaptation_->ReferenceSmoothingLength()) {} + //=================================================================================================// + Real BarAcousticTimeStepSize::reduce(size_t index_i, Real dt) + { + // Since the particle does not change its configuration in pressure relaxation step, + // I chose a time-step size according to Eulerian method. + Real time_setp_0 = SMIN(sqrt(smoothing_length_ / ((acc_[index_i] + acc_prior_[index_i]).norm() + TinyReal)), + smoothing_length_ / (c0_ + vel_[index_i].norm())); + Real time_setp_1 = SMIN(sqrt(1.0 / (dangular_vel_dt_[index_i].norm() + TinyReal)), + 1.0 / (angular_vel_[index_i].norm() + TinyReal)); + Real time_setp_2 = smoothing_length_ * sqrt(rho0_ * (1.0 - nu_ * nu_) / E0_ / + (2.0 + (Pi * Pi / 12.0) * (1.0 - nu_) * + (1.0 + 1.5 * pow(smoothing_length_ / thickness_[index_i], 2)))); + return CFL_ * SMIN(time_setp_0, time_setp_1, time_setp_2); + } + //=================================================================================================// + BarCorrectConfiguration:: + BarCorrectConfiguration(BaseInnerRelation &inner_relation) + : LocalDynamics(inner_relation.getSPHBody()), BarDataInner(inner_relation), + B_(particles_->B_), + n0_(particles_->n0_), transformation_matrix_(particles_->transformation_matrix_), b_n0_(particles_->b_n0_) {} + //=================================================================================================// + BarDeformationGradientTensor:: + BarDeformationGradientTensor(BaseInnerRelation &inner_relation) + : LocalDynamics(inner_relation.getSPHBody()), BarDataInner(inner_relation), + pos_(particles_->pos_), pseudo_n_(particles_->pseudo_n_), n0_(particles_->n0_), + pseudo_b_n_(particles_->pseudo_n_), b_n0_(particles_->n0_), + B_(particles_->B_), F_(particles_->F_), F_bending_(particles_->F_bending_), + F_b_bending_(particles_->F_b_bending_), + transformation_matrix_(particles_->transformation_matrix_) {} + //=================================================================================================// + BaseBarRelaxation::BaseBarRelaxation(BaseInnerRelation &inner_relation) + : LocalDynamics(inner_relation.getSPHBody()), BarDataInner(inner_relation), + rho_(particles_->rho_), + thickness_(particles_->thickness_), + pos_(particles_->pos_), vel_(particles_->vel_), + acc_(particles_->acc_), + acc_prior_(particles_->acc_prior_), + n0_(particles_->n0_), pseudo_n_(particles_->pseudo_n_), + dpseudo_n_dt_(particles_->dpseudo_n_dt_), dpseudo_n_d2t_(particles_->dpseudo_n_d2t_), + rotation_(particles_->rotation_), angular_vel_(particles_->angular_vel_), + dangular_vel_dt_(particles_->dangular_vel_dt_), + B_(particles_->B_), F_(particles_->F_), dF_dt_(particles_->dF_dt_), + F_bending_(particles_->F_bending_), dF_bending_dt_(particles_->dF_bending_dt_), + transformation_matrix_(particles_->transformation_matrix_), + F_b_bending_(particles_->F_b_bending_), dF_b_bending_dt_(particles_->dF_b_bending_dt_), + b_n0_(particles_->b_n0_), pseudo_b_n_(particles_->pseudo_b_n_), + dpseudo_b_n_dt_(particles_->dpseudo_b_n_dt_), dpseudo_b_n_d2t_(particles_->dpseudo_b_n_d2t_), + rotation_b_(particles_->rotation_b_), angular_b_vel_(particles_->angular_b_vel_), + dangular_b_vel_dt_(particles_->dangular_b_vel_dt_), width_(particles_->width_) {} + //=================================================================================================// + BarStressRelaxationFirstHalf:: + BarStressRelaxationFirstHalf(BaseInnerRelation &inner_relation, + int number_of_gaussian_points, bool hourglass_control) + : BaseBarRelaxation(inner_relation), + elastic_solid_(particles_->elastic_solid_), + global_stress_(particles_->global_stress_), + global_moment_(particles_->global_moment_), + mid_surface_cauchy_stress_(particles_->mid_surface_cauchy_stress_), + numerical_damping_scaling_(particles_->numerical_damping_scaling_), + global_shear_stress_(particles_->global_shear_stress_), + n_(particles_->n_), + rho0_(elastic_solid_.ReferenceDensity()), + inv_rho0_(1.0 / rho0_), + smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()), + E0_(elastic_solid_.YoungsModulus()), + G0_(elastic_solid_.ShearModulus()), + nu_(elastic_solid_.PoissonRatio()), + hourglass_control_(hourglass_control), + number_of_gaussian_points_(number_of_gaussian_points), + global_b_stress_(particles_->global_b_stress_), + global_b_moment_(particles_->global_b_moment_), + global_b_shear_stress_(particles_->global_b_shear_stress_), + b_n_(particles_->b_n_) + { + + /** Note that, only three-point and five-point Gaussian quadrature rules are defined. */ + switch (number_of_gaussian_points) + { + + case 9: + gaussian_point_x = nine_gaussian_points_2d_vector_x; + gaussian_point_y = nine_gaussian_points_2d_vector_y; + gaussian_weight_ = nine_gaussian_weights_2d_; + break; + default: + gaussian_point_x = four_gaussian_points_2d_vector_x; + gaussian_point_y = four_gaussian_points_2d_vector_y; + gaussian_weight_ = four_gaussian_weights_2d_; + } + /** Define the factor of hourglass control algorithm. */ + hourglass_control_factor_ = 0.005; + } + //=================================================================================================// + void BarStressRelaxationFirstHalf::initialization(size_t index_i, Real dt) + { + // Note that F_[index_i], F_bending_[index_i], dF_dt_[index_i], dF_bending_dt_[index_i] + // and rotation_[index_i], angular_vel_[index_i], dangular_vel_dt_[index_i], B_[index_i] + // are defined in local coordinates, while others in global coordinates. + pos_[index_i] += vel_[index_i] * dt * 0.5; + rotation_[index_i] += angular_vel_[index_i] * dt * 0.5; + pseudo_n_[index_i] += dpseudo_n_dt_[index_i] * dt * 0.5; + + rotation_b_[index_i] += angular_b_vel_[index_i] * dt * 0.5; + pseudo_b_n_[index_i] += dpseudo_b_n_dt_[index_i] * dt * 0.5; + + F_[index_i] += dF_dt_[index_i] * dt * 0.5; + F_bending_[index_i] += dF_bending_dt_[index_i] * dt * 0.5; + F_b_bending_[index_i] += dF_b_bending_dt_[index_i] * dt * 0.5; + + Real J = F_[index_i].determinant(); + Matd inverse_F = F_[index_i].inverse(); + + rho_[index_i] = rho0_ / J; + + /** Calculate the current normal direction of mid-surface. */ + n_[index_i] = transformation_matrix_[index_i].transpose() * getNormalFromDeformationGradientTensor(F_[index_i]); + b_n_[index_i] = transformation_matrix_[index_i].transpose() * getBinormalFromDeformationGradientTensor(F_[index_i]); + + /** Get transformation matrix from global coordinates to current local coordinates. */ + Matd current_transformation_matrix = getTransformationMatrix(n_[index_i], b_n_[index_i]); + + Matd resultant_stress = Matd::Zero(); + Matd resultant_moment = Matd::Zero(); + Vecd resultant_shear_stress = Vecd::Zero(); + Matd resultant_b_stress = Matd::Zero(); + Matd resultant_b_moment = Matd::Zero(); + Vecd resultant_b_shear_stress = Vecd::Zero(); + + for (int i = 0; i != number_of_gaussian_points_; ++i) + { + Matd F_gaussian_point = F_[index_i] + gaussian_point_y[i] * F_bending_[index_i] * thickness_[index_i] * 0.5+gaussian_point_x[i]*F_b_bending_[index_i]*width_[index_i]*0.5; + Matd dF_gaussian_point_dt = dF_dt_[index_i] + gaussian_point_y[i] * dF_bending_dt_[index_i] * thickness_[index_i] * 0.5 + gaussian_point_x[i] * dF_b_bending_dt_[index_i] * width_[index_i] * 0.5; + Matd inverse_F_gaussian_point = F_gaussian_point.inverse(); + Matd current_local_almansi_strain = current_transformation_matrix * transformation_matrix_[index_i].transpose() * 0.5 * + (Matd::Identity() - inverse_F_gaussian_point.transpose() * inverse_F_gaussian_point) * + transformation_matrix_[index_i] * current_transformation_matrix.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. */ + Matd cauchy_stress = elastic_solid_.StressCauchy(current_local_almansi_strain, F_gaussian_point, index_i) + + current_transformation_matrix * transformation_matrix_[index_i].transpose() * F_gaussian_point * + elastic_solid_.NumericalDampingRightCauchy(F_gaussian_point, dF_gaussian_point_dt, numerical_damping_scaling_[index_i], index_i) + * F_gaussian_point.transpose() * transformation_matrix_[index_i] * current_transformation_matrix.transpose() + / F_gaussian_point.determinant(); + + /** Impose modeling assumptions. */ + cauchy_stress(0, 1) *= shear_correction_factor_; + cauchy_stress(1, 0) *= shear_correction_factor_; + cauchy_stress(0, 2) *= shear_correction_factor_; + cauchy_stress(2, 0) *= shear_correction_factor_; + cauchy_stress(1, 2) *= shear_correction_factor_; + cauchy_stress(2, 1) *= shear_correction_factor_; + + cauchy_stress(Dimensions - 1, Dimensions - 1) = 0.0; + cauchy_stress(Dimensions - 2, Dimensions - 2) = 0.0; + if (i == 0) + { + mid_surface_cauchy_stress_[index_i] = cauchy_stress; + } + + /** Integrate Cauchy stress along thickness. */ + resultant_stress += + 0.5 * width_[index_i] * 0.5 * thickness_[index_i] * gaussian_weight_[i] * cauchy_stress; + resultant_moment += + 0.5 * width_[index_i] * 0.5 * thickness_[index_i] * gaussian_weight_[i] * ((cauchy_stress)*gaussian_point_y[i] * thickness_[index_i] * 0.5); + resultant_b_moment += + 0.5 * width_[index_i] * 0.5 * thickness_[index_i] * gaussian_weight_[i] * ((cauchy_stress)*gaussian_point_x[i] * width_[index_i] * 0.5); + resultant_shear_stress -= + 0.5 * width_[index_i] * 0.5 * thickness_[index_i] * gaussian_weight_[i] * cauchy_stress.col(Dimensions - 1); + resultant_b_shear_stress -= + 0.5 * width_[index_i] * 0.5 * thickness_[index_i] * gaussian_weight_[i] * cauchy_stress.col(Dimensions - 2); + + + resultant_stress.row(Dimensions - 1) = Vecd::Zero().transpose(); + resultant_stress.col(Dimensions - 1) = Vecd::Zero(); + + resultant_stress.row(Dimensions - 2) = Vecd::Zero().transpose(); + resultant_stress.col(Dimensions - 2) = Vecd::Zero(); + + resultant_moment.row(Dimensions - 1) = Vecd::Zero().transpose(); + resultant_moment.col(Dimensions - 1) = Vecd::Zero(); + + resultant_b_moment.row(Dimensions - 2) = Vecd::Zero().transpose(); + resultant_b_moment.col(Dimensions - 2) = Vecd::Zero(); + + resultant_shear_stress[Dimensions - 2] = 0.0; + resultant_b_shear_stress[Dimensions - 1] = 0.0; + } + + /** stress and moment in global coordinates for pair interaction */ + global_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_stress * current_transformation_matrix * transformation_matrix_[index_i].transpose() * inverse_F.transpose() * transformation_matrix_[index_i]; + global_moment_[index_i] = J * current_transformation_matrix.transpose() * resultant_moment * current_transformation_matrix * transformation_matrix_[index_i].transpose() * inverse_F.transpose() * transformation_matrix_[index_i]; + global_shear_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_shear_stress; + + global_b_moment_[index_i] = J * current_transformation_matrix.transpose() * resultant_b_moment * current_transformation_matrix * transformation_matrix_[index_i].transpose() * inverse_F.transpose() * transformation_matrix_[index_i]; + global_b_shear_stress_[index_i] = J * current_transformation_matrix.transpose() * resultant_b_shear_stress; + } + //=================================================================================================// + 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; + } + //=================================================================================================// + void BarStressRelaxationSecondHalf::initialization(size_t index_i, Real dt) + { + pos_[index_i] += vel_[index_i] * dt * 0.5; + rotation_b_[index_i] += angular_b_vel_[index_i] * dt * 0.5; + rotation_[index_i] += angular_vel_[index_i] * dt * 0.5; + Vecd pseudo_n_temp = pseudo_n_[index_i]; + Vecd pseudo_b_n_temp = pseudo_b_n_[index_i]; + + pseudo_n_[index_i] = transformation_matrix_[index_i].transpose() * getVectorAfterThinStructureRotation(local_pseudo_n_0, rotation_[index_i] ); + + pseudo_b_n_[index_i] = transformation_matrix_[index_i].transpose() * getVectorAfterThinStructureRotation(local_pseudo_b_n_0, rotation_b_[index_i] ); + + if (dt < 1e-10) + { + dpseudo_n_dt_[index_i] = Vecd::Zero(); + dpseudo_b_n_dt_[index_i] = Vecd::Zero(); + } + else + { + dpseudo_n_dt_[index_i] = (pseudo_n_[index_i] - pseudo_n_temp) / (0.5 * dt); + dpseudo_b_n_dt_[index_i] = (pseudo_b_n_[index_i] - pseudo_b_n_temp) / (0.5 * dt); + } + } + //=================================================================================================// + void BarStressRelaxationSecondHalf::update(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; + F_b_bending_[index_i] += dF_b_bending_dt_[index_i] * dt * 0.5; + + } + //=================================================================================================// + ConstrainBarBodyRegion:: + ConstrainBarBodyRegion(BodyPartByParticle &body_part) + : BaseLocalDynamics(body_part), BarDataSimple(sph_body_), + vel_(particles_->vel_), angular_vel_(particles_->angular_vel_), angular_b_vel_(particles_->angular_b_vel_) {} + //=================================================================================================// + void ConstrainBarBodyRegion::update(size_t index_i, Real dt) + { + vel_[index_i] = Vecd::Zero(); + angular_vel_[index_i] = Vecd::Zero(); + angular_b_vel_[index_i] = Vecd::Zero(); + } + //=================================================================================================// + ConstrainBarBodyRegionAlongAxis::ConstrainBarBodyRegionAlongAxis(BodyPartByParticle &body_part, int axis) + : BaseLocalDynamics(body_part), BarDataSimple(sph_body_), + axis_(axis), pos_(particles_->pos_), pos0_(particles_->pos0_), vel_(particles_->vel_), + acc_(particles_->acc_), rotation_(particles_->rotation_), angular_vel_(particles_->angular_vel_), + dangular_vel_dt_(particles_->dangular_vel_dt_), rotation_b_(particles_->rotation_b_), angular_b_vel_(particles_->angular_b_vel_), dangular_b_vel_dt_(particles_->dangular_b_vel_dt_) {} + //=================================================================================================// + void ConstrainBarBodyRegionAlongAxis::update(size_t index_i, Real dt) + { + vel_[index_i] = Vecd::Zero(); + } + //=================================================================================================// + DistributingPointForcesToBar:: + DistributingPointForcesToBar(SPHBody &sph_body, std::vector point_forces, + std::vector reference_positions, Real time_to_full_external_force, + Real particle_spacing_ref, Real h_spacing_ratio) + : LocalDynamics(sph_body), BarDataSimple(sph_body), + point_forces_(point_forces), reference_positions_(reference_positions), + time_to_full_external_force_(time_to_full_external_force), + particle_spacing_ref_(particle_spacing_ref), h_spacing_ratio_(h_spacing_ratio), + pos0_(particles_->pos0_), acc_prior_(particles_->acc_prior_), + thickness_(particles_->thickness_) + { + for (size_t i = 0; i < point_forces_.size(); i++) + { + weight_.push_back(StdLargeVec(0.0)); + time_dependent_point_forces_.push_back(Vecd::Zero()); + sum_of_weight_.push_back(0.0); + particles_->registerVariable(weight_[i], "Weight_" + std::to_string(i)); + } + + getWeight(); // TODO: should be revised and parallelized, using SimpleDynamics + } + //=================================================================================================// + void DistributingPointForcesToBar::getWeight() + { + Kernel *kernel_ = sph_body_.sph_adaptation_->getKernel(); + Real reference_smoothing_length = sph_body_.sph_adaptation_->ReferenceSmoothingLength(); + Real smoothing_length = h_spacing_ratio_ * particle_spacing_ref_; + Real h_ratio = reference_smoothing_length / smoothing_length; + Real cutoff_radius_sqr = pow(2.0 * smoothing_length, 2); + for (size_t i = 0; i < point_forces_.size(); ++i) + { + sum_of_weight_[i] = 0.0; + for (size_t index = 0; index < particles_->total_real_particles_; ++index) + { + weight_[i][index] = 0.0; + Vecd displacement = reference_positions_[i] - pos0_[index]; + if (displacement.squaredNorm() <= cutoff_radius_sqr) + { + weight_[i][index] = kernel_->W(h_ratio, displacement.norm(), displacement); + sum_of_weight_[i] += weight_[i][index]; + } + } + } + } + //=================================================================================================// + void DistributingPointForcesToBar::setupDynamics(Real dt) + { + Real current_time = GlobalStaticVariables::physical_time_; + for (size_t i = 0; i < point_forces_.size(); ++i) + { + time_dependent_point_forces_[i] = current_time < time_to_full_external_force_ + ? current_time * point_forces_[i] / time_to_full_external_force_ + : point_forces_[i]; + } + } + //=================================================================================================// + void DistributingPointForcesToBar::update(size_t index_i, Real dt) + { + acc_prior_[index_i] = Vecd::Zero(); + for (size_t i = 0; i < point_forces_.size(); ++i) + { + Vecd force = weight_[i][index_i] / (sum_of_weight_[i] + TinyReal) * time_dependent_point_forces_[i]; + acc_prior_[index_i] += force / particles_->ParticleMass(index_i); + } + } + //=================================================================================================// + } +} \ No newline at end of file diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h new file mode 100644 index 0000000000..b16d63c015 --- /dev/null +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_dynamics.h @@ -0,0 +1,408 @@ +/* -------------------------------------------------------------------------* + * SPHinXsys * + * -------------------------------------------------------------------------* + * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle* + * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * + * physical accurate simulation and aims to model coupled industrial dynamic* + * systems including fluid, solid, multi-body dynamics and beyond with SPH * + * (smoothed particle hydrodynamics), a meshless computational method using * + * particle discretization. * + * * + * SPHinXsys is partially funded by German Research Foundation * + * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * + * HU1527/12-1 and HU1527/12-4 * + * * + * Portions copyright (c) 2017-2022 Technical University of Munich and * + * the authors' affiliations. * + * * + * Licensed under the Apache License, Version 2.0 (the "License"); you may * + * not use this file except in compliance with the License. You may obtain a* + * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * + * * + * ------------------------------------------------------------------------*/ +/** + * @file thin_structure_dynamics.h + * @brief Here, we define the algorithm classes for slender structure dynamics based on thin structure dynamics. + * @details We consider here a weakly compressible solids. + * @author Xipeng Lyu + */ + +#ifndef SLENDER_STRUCTURE_DYNAMICS_H +#define SLENDER_STRUCTURE_DYNAMICS_H + +#include "all_particle_dynamics.h" +#include "base_kernel.h" +#include "all_body_relations.h" +#include "solid_body.h" +#include "solid_particles.h" +#include "elastic_solid.h" +#include "slender_structure_math.h" + +namespace SPH +{ + namespace slender_structure_dynamics + { + typedef DataDelegateSimple BarDataSimple; + typedef DataDelegateInner BarDataInner; + + /** + * @class BarDynamicsInitialCondition + * @brief set initial condition for bar particles + * This is a abstract class to be override for case specific initial conditions. + */ + class BarDynamicsInitialCondition : public LocalDynamics, public BarDataSimple + { + public: + explicit BarDynamicsInitialCondition(SPHBody &sph_body); + virtual ~BarDynamicsInitialCondition(){}; + + protected: + StdLargeVec &b_n0_, &b_n_, &pseudo_b_n_; + StdLargeVec &n0_, &n_, &pseudo_n_, &pos0_; + StdLargeVec &transformation_matrix_; + }; + + /** + * @class BarAcousticTimeStepSize + * @brief Computing the acoustic time step size for bar + */ + class BarAcousticTimeStepSize : public LocalDynamicsReduce, + public BarDataSimple + { + protected: + Real CFL_; + StdLargeVec &vel_, &acc_, &angular_vel_, &dangular_vel_dt_, &acc_prior_; + StdLargeVec &thickness_; + Real rho0_, E0_, nu_, c0_; + Real smoothing_length_; + + StdLargeVec &angular_b_vel_, &dangular_b_vel_dt_; + StdLargeVec &width_; + + public: + explicit BarAcousticTimeStepSize(SPHBody &sph_body, Real CFL = 0.6); + virtual ~BarAcousticTimeStepSize(){}; + + Real reduce(size_t index_i, Real dt = 0.0); + }; + + /** + * @class BarCorrectConfiguration + * @brief obtain the corrected initial configuration in strong form + */ + class BarCorrectConfiguration : public LocalDynamics, public BarDataInner + { + public: + explicit BarCorrectConfiguration(BaseInnerRelation &inner_relation); + virtual ~BarCorrectConfiguration(){}; + + inline void interaction(size_t index_i, Real dt = 0.0) + { + /** A small number is added to diagonal to avoid dividing by zero. */ + Matd global_configuration = Eps * Matd::Identity(); + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + Vecd gradW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + Vecd r_ji = -inner_neighborhood.r_ij_[n] * inner_neighborhood.e_ij_[n]; + global_configuration += r_ji * gradW_ijV_j.transpose(); + } + Matd local_configuration = + transformation_matrix_[index_i] * global_configuration * transformation_matrix_[index_i].transpose(); + /** correction matrix is obtained from local configuration. */ + B_[index_i] = getCorrectionMatrix_beam(local_configuration); + }; + + protected: + StdLargeVec &B_; + StdLargeVec &n0_; + StdLargeVec &b_n0_; + StdLargeVec &transformation_matrix_; + }; + + /** + * @class BarDeformationGradientTensor + * @brief computing deformation gradient tensor for bar + * TODO: need a test case for this. + */ + class BarDeformationGradientTensor : public LocalDynamics, public BarDataInner + { + public: + explicit BarDeformationGradientTensor(BaseInnerRelation &inner_relation); + virtual ~BarDeformationGradientTensor(){}; + + inline void interaction(size_t index_i, Real dt = 0.0) + { + const Vecd &pseudo_n_i = pseudo_n_[index_i]; + const Vecd &pseudo_b_n_i = pseudo_b_n_[index_i]; + const Vecd &pos_n_i = pos_[index_i]; + const Matd &transformation_matrix_i = transformation_matrix_[index_i]; + + Matd deformation_part_one = Matd::Zero(); + Matd deformation_part_two = Matd::Zero(); + Matd deformation_part_three = Matd::Zero(); + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + Vecd gradW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + deformation_part_one -= (pos_n_i - pos_[index_j]) * gradW_ijV_j.transpose(); + deformation_part_two -= ((pseudo_n_i - n0_[index_i]) - (pseudo_n_[index_j] - n0_[index_j])) * gradW_ijV_j.transpose(); + deformation_part_three -= ((pseudo_b_n_i - b_n0_[index_i]) - (pseudo_b_n_[index_j] - b_n0_[index_j])) * gradW_ijV_j.transpose(); + } + F_[index_i] = transformation_matrix_i * deformation_part_one * transformation_matrix_i.transpose() * B_[index_i]; + F_[index_i].col(Dimensions - 1) = transformation_matrix_i * pseudo_n_[index_i]; + F_[index_i].col(Dimensions - 2) = transformation_matrix_i * pseudo_b_n_[index_i]; + F_bending_[index_i] = transformation_matrix_i * deformation_part_two * transformation_matrix_i.transpose() * B_[index_i]; + F_b_bending_[index_i] = transformation_matrix_i * deformation_part_three * transformation_matrix_i.transpose() * B_[index_i]; + }; + + protected: + StdLargeVec &pos_, &pseudo_n_, &n0_; + StdLargeVec &B_, &F_, &F_bending_; + StdLargeVec &transformation_matrix_; + + StdLargeVec &pseudo_b_n_, &b_n0_; + StdLargeVec &F_b_bending_; + }; + + /** + * @class BaseBarRelaxation + * @brief abstract class for preparing bar relaxation + */ + class BaseBarRelaxation : public LocalDynamics, public BarDataInner + { + public: + explicit BaseBarRelaxation(BaseInnerRelation &inner_relation); + virtual ~BaseBarRelaxation(){}; + + protected: + StdLargeVec &rho_, &thickness_; + StdLargeVec &pos_, &vel_, &acc_, &acc_prior_; + StdLargeVec &n0_, &pseudo_n_, &dpseudo_n_dt_, &dpseudo_n_d2t_, &rotation_, + &angular_vel_, &dangular_vel_dt_; + StdLargeVec &B_, &F_, &dF_dt_, &F_bending_, &dF_bending_dt_; + StdLargeVec &transformation_matrix_; + + StdLargeVec &width_; + StdLargeVec &b_n0_, &pseudo_b_n_, &dpseudo_b_n_dt_, &dpseudo_b_n_d2t_, &rotation_b_, + &angular_b_vel_, dangular_b_vel_dt_; + StdLargeVec &F_b_bending_, &dF_b_bending_dt_; + }; + + /** + * @class BarStressRelaxationFirstHalf + * @brief computing stress relaxation process by Verlet time stepping + * This is the first step + */ + class BarStressRelaxationFirstHalf : public BaseBarRelaxation + { + public: + explicit BarStressRelaxationFirstHalf(BaseInnerRelation &inner_relation, + int number_of_gaussian_points = 4, bool hourglass_control = false); + virtual ~BarStressRelaxationFirstHalf(){}; + void initialization(size_t index_i, Real dt = 0.0); + + inline void interaction(size_t index_i, Real dt = 0.0) + { + const Vecd &global_shear_stress_i = global_shear_stress_[index_i]; + const Matd &global_stress_i = global_stress_[index_i]; + const Matd &global_moment_i = global_moment_[index_i]; + const Matd &global_b_moment_i = global_b_moment_[index_i]; + const Vecd &global_b_shear_stress_i = global_b_shear_stress_[index_i]; + + Vecd acceleration = Vecd::Zero(); + Vecd pseudo_normal_acceleration = global_shear_stress_i; + Vecd pseudo_b_normal_acceleration = global_b_shear_stress_i; + + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + + if (hourglass_control_) + { + Vecd e_ij = inner_neighborhood.e_ij_[n]; + Real r_ij = inner_neighborhood.r_ij_[n]; + Real dim_inv_r_ij = Dimensions / r_ij; + Real weight = inner_neighborhood.W_ij_[n] * inv_W0_; + Vecd pos_jump = getLinearVariableJump(e_ij, r_ij, pos_[index_i], + transformation_matrix_[index_i].transpose() * F_[index_i] * transformation_matrix_[index_i], + pos_[index_j], + transformation_matrix_[index_i].transpose() * F_[index_j] * transformation_matrix_[index_i]); + acceleration += hourglass_control_factor_ * weight * E0_ * pos_jump * dim_inv_r_ij * + inner_neighborhood.dW_ijV_j_[n] * thickness_[index_i]; + + Vecd pseudo_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_n_[index_i] - n0_[index_i], + transformation_matrix_[index_i].transpose() * F_bending_[index_i] * transformation_matrix_[index_i], + pseudo_n_[index_j] - n0_[index_j], + transformation_matrix_[index_j].transpose() * F_bending_[index_j] * transformation_matrix_[index_j]); + Vecd rotation_jump = getRotationJump(pseudo_n_jump, transformation_matrix_[index_i]); + pseudo_normal_acceleration += hourglass_control_factor_ * weight * G0_ * rotation_jump * dim_inv_r_ij * + inner_neighborhood.dW_ijV_j_[n] * pow(thickness_[index_i], 4); + + Vecd pseudo_b_n_jump = getLinearVariableJump(e_ij, r_ij, pseudo_b_n_[index_i] - b_n0_[index_i], + transformation_matrix_[index_i].transpose() * F_b_bending_[index_i] * transformation_matrix_[index_i], + pseudo_b_n_[index_j] - b_n0_[index_j], + transformation_matrix_[index_j].transpose() * F_b_bending_[index_j] * transformation_matrix_[index_j]); + Vecd rotation_b_jump = getRotationJump(pseudo_b_n_jump, transformation_matrix_[index_i]); + pseudo_b_normal_acceleration += hourglass_control_factor_ * weight * G0_ * rotation_b_jump * dim_inv_r_ij * + inner_neighborhood.dW_ijV_j_[n] * pow(thickness_[index_i], 4); + } + acceleration += (global_stress_i + global_stress_[index_j]) * inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + pseudo_normal_acceleration += (global_moment_i + global_moment_[index_j]) * inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + pseudo_b_normal_acceleration += (global_b_moment_i + global_b_moment_[index_j]) * inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + } + + 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]; + dangular_b_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation_b(local_dpseudo_b_n_d2t, local_dpseudo_n_d2t, rotation_b_[index_i], angular_b_vel_[index_i], dt); + dangular_vel_dt_[index_i] = getRotationFromPseudoNormalForSmallDeformation(local_dpseudo_b_n_d2t, local_dpseudo_n_d2t, rotation_[index_i], angular_vel_[index_i], dt); + }; + + void update(size_t index_i, Real dt = 0.0); + + protected: + ElasticSolid &elastic_solid_; + StdLargeVec &global_stress_, &global_moment_, &mid_surface_cauchy_stress_, &numerical_damping_scaling_; + StdLargeVec &global_shear_stress_, &n_; + + StdLargeVec &global_b_moment_, &global_b_stress_; + StdLargeVec &b_n_, &global_b_shear_stress_; + Real rho0_, inv_rho0_; + Real smoothing_length_, E0_, G0_, nu_, hourglass_control_factor_; + bool hourglass_control_; + const Real inv_W0_ = 1.0 / sph_body_.sph_adaptation_->getKernel()->W0(ZeroVecd); + const Real shear_correction_factor_ = 5.0 / 6.0; + + + Real gpt = sqrt(3.0 / 5.0); + const StdVec nine_gaussian_points_2d_vector_x{-gpt, gpt, gpt, -gpt, 0, gpt, 0, -gpt, 0}; + const StdVec nine_gaussian_points_2d_vector_y{-gpt, -gpt, gpt, gpt, -gpt, 0, gpt, 0, 0}; + + const StdVec nine_gaussian_weights_2d_ = {25. / 81., 25. / 81., 25. / 81., 25. / 81., 40. / 81., 40. / 81., 40. / 81., 40. / 81., 64. / 81.}; + + Real gpt_4 = sqrt(1. / 3.); + const StdVec four_gaussian_points_2d_vector_x{-gpt_4, gpt_4, gpt_4, -gpt_4}; + const StdVec four_gaussian_points_2d_vector_y{-gpt_4, -gpt_4, gpt_4, gpt_4}; + + const StdVec four_gaussian_weights_2d_ = {1, 1, 1, 1}; + + + int number_of_gaussian_points_; + StdVec gaussian_point_x; + StdVec gaussian_point_y; + StdVec gaussian_weight_; + }; + + /** + * @class BarStressRelaxationSecondHalf + * @brief computing stress relaxation process by verlet time stepping + * This is the second step + */ + class BarStressRelaxationSecondHalf : public BaseBarRelaxation + { + public: + explicit BarStressRelaxationSecondHalf(BaseInnerRelation &inner_relation) + : BaseBarRelaxation(inner_relation){}; + virtual ~BarStressRelaxationSecondHalf(){}; + void initialization(size_t index_i, Real dt = 0.0); + + inline void interaction(size_t index_i, Real dt = 0.0) + { + const Vecd &vel_n_i = vel_[index_i]; + const Vecd &dpseudo_n_dt_i = dpseudo_n_dt_[index_i]; + const Vecd &dpseudo_b_n_dt_i = dpseudo_b_n_dt_[index_i]; + const Matd &transformation_matrix_i = transformation_matrix_[index_i]; + + Matd deformation_gradient_change_rate_part_one = Matd::Zero(); + Matd deformation_gradient_change_rate_part_three = Matd::Zero(); + Matd deformation_gradient_change_rate_part_two = Matd::Zero(); + const Neighborhood &inner_neighborhood = inner_configuration_[index_i]; + for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) + { + size_t index_j = inner_neighborhood.j_[n]; + + Vecd gradW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n]; + deformation_gradient_change_rate_part_one -= (vel_n_i - vel_[index_j]) * gradW_ijV_j.transpose(); + deformation_gradient_change_rate_part_two -= (dpseudo_n_dt_i - dpseudo_n_dt_[index_j]) * gradW_ijV_j.transpose(); + deformation_gradient_change_rate_part_three -= (dpseudo_b_n_dt_i - dpseudo_b_n_dt_[index_j]) * gradW_ijV_j.transpose(); + + } + dF_dt_[index_i] = transformation_matrix_i * deformation_gradient_change_rate_part_one * transformation_matrix_i.transpose() * B_[index_i]; + dF_dt_[index_i].col(Dimensions - 1) = transformation_matrix_i * dpseudo_n_dt_[index_i]; + dF_dt_[index_i].col(Dimensions - 2) = transformation_matrix_i * dpseudo_b_n_dt_[index_i]; + dF_bending_dt_[index_i] = transformation_matrix_i * deformation_gradient_change_rate_part_two * transformation_matrix_i.transpose() * B_[index_i]; + dF_b_bending_dt_[index_i] = transformation_matrix_i * deformation_gradient_change_rate_part_three * transformation_matrix_i.transpose() * B_[index_i]; + }; + + void update(size_t index_i, Real dt = 0.0); + }; + + /**@class ConstrainBarBodyRegion + * @brief Fix the position and angle of a bar body part. + */ + class ConstrainBarBodyRegion : public BaseLocalDynamics, public BarDataSimple + { + public: + ConstrainBarBodyRegion(BodyPartByParticle &body_part); + virtual ~ConstrainBarBodyRegion(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + StdLargeVec &vel_, &angular_vel_, &angular_b_vel_; + }; + + /**@class ConstrainBarBodyRegionAlongAxis + * @brief The boundary conditions are denoted by SS1 according to the references. + * The axis must be 0 or 1. + * Note that the average values for FSI are prescribed also. + */ + class ConstrainBarBodyRegionAlongAxis : public BaseLocalDynamics, public BarDataSimple + { + public: + ConstrainBarBodyRegionAlongAxis(BodyPartByParticle &body_part, int axis); + virtual ~ConstrainBarBodyRegionAlongAxis(){}; + void update(size_t index_i, Real dt = 0.0); + + protected: + const int axis_; /**< the axis direction for bounding*/ + StdLargeVec &pos_, &pos0_; + StdLargeVec &vel_, &acc_; + StdLargeVec &rotation_b_, &angular_b_vel_, &dangular_b_vel_dt_; + StdLargeVec &rotation_, &angular_vel_, &dangular_vel_dt_; + }; + + /** + * @class DistributingPointForcesToBar + * @brief Distribute a series of point forces to its contact Bar bodies. + */ + class DistributingPointForcesToBar : public LocalDynamics, public BarDataSimple + { + protected: + std::vector point_forces_, reference_positions_, time_dependent_point_forces_; + Real time_to_full_external_force_; + Real particle_spacing_ref_, h_spacing_ratio_; + StdLargeVec &pos0_, &acc_prior_; + StdLargeVec &thickness_; + std::vector> weight_; + std::vector sum_of_weight_; + + void getWeight(); + + public: + DistributingPointForcesToBar(SPHBody &sph_body, std::vector point_forces, + std::vector reference_positions, Real time_to_full_external_force, + Real particle_spacing_ref, Real h_spacing_ratio = 1.6); + virtual ~DistributingPointForcesToBar(){}; + + virtual void setupDynamics(Real dt = 0.0) override; + void update(size_t index_i, Real dt = 0.0); + }; + } +} +#endif // THIN_STRUCTURE_DYNAMICS_H diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp new file mode 100644 index 0000000000..c1aef91895 --- /dev/null +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.cpp @@ -0,0 +1,264 @@ +#include "slender_structure_math.h" + +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]); + Real cos_angle_x = cos(rotation_angles[0]); + + Real sin_angle_y = sin(rotation_angles[1]); + Real cos_angle_y = cos(rotation_angles[1]); + + Real sin_angle_z = sin(rotation_angles[2]); + Real cos_angle_z = cos(rotation_angles[2]); + + 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; + } + //=================================================================================================// + 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) + { + Real sin_rotation_0 = sin(rotation_angles[0]); + Real cos_rotation_0 = cos(rotation_angles[0]); + + Real sin_rotation_1 = sin(rotation_angles[1]); + Real cos_rotation_1 = cos(rotation_angles[1]); + + Real dpseudo_n_dt_0 = -sin_rotation_0 * sin_rotation_1 * angular_vel[0] + cos_rotation_0 * cos_rotation_1 * angular_vel[1]; + Real dpseudo_n_dt_1 = -cos_rotation_0 * angular_vel[0]; + Real dpseudo_n_dt_2 = -sin_rotation_0 * cos_rotation_1 * angular_vel[0] - cos_rotation_0 * sin_rotation_1 * angular_vel[1]; + + 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]); + Real cos_rotation_0 = cos(rotation[0]); + Real sin_rotation_1 = sin(rotation[1]); + Real cos_rotation_1 = cos(rotation[1]); + + Real rotation_0_a = -(dpseudo_n_d2t[2] * cos_rotation_1 + dpseudo_n_d2t[0] * sin_rotation_1 + + angular_vel[1] * angular_vel[1] * cos_rotation_0 + + angular_vel[0] * angular_vel[0] * cos_rotation_0); + Real rotation_0_b = sin_rotation_0 * angular_vel[0] * angular_vel[0] - dpseudo_n_d2t[1]; + Real angle_vel_dt_0 = sin_rotation_0 * rotation_0_a + cos_rotation_0 * rotation_0_b; + + Real rotation_1_a = dpseudo_n_d2t[0] * cos_rotation_1 - dpseudo_n_d2t[2] * sin_rotation_1 + + 2.0 * angular_vel[1] * angular_vel[0] * sin_rotation_0; + Real rotation_1_b1 = dpseudo_n_d2t[0] * cos_rotation_0 + + angular_vel[1] * angular_vel[1] * cos_rotation_0 * cos_rotation_0 * sin_rotation_1 + + angular_vel[0] * angular_vel[0] * sin_rotation_1 + - dpseudo_n_d2t[1] * sin_rotation_1 * sin_rotation_0 + + 2.0 * angular_vel[1] * angular_vel[0] * cos_rotation_1 * cos_rotation_0 * sin_rotation_0; + Real rotation_1_b2 = -(dpseudo_n_d2t[2] * cos_rotation_0 + + angular_vel[1] * angular_vel[1] * cos_rotation_1 * cos_rotation_0 * cos_rotation_0 + + angular_vel[0] * angular_vel[0] * cos_rotation_1 + - dpseudo_n_d2t[1] * cos_rotation_1 * sin_rotation_0 + - 2.0 * angular_vel[1] * angular_vel[0] * cos_rotation_0 * sin_rotation_1 * sin_rotation_0); + Real angle_vel_dt_1 = rotation_1_a * rotation_1_a * (rotation_1_b1 * cos_rotation_1 + rotation_1_b2 * sin_rotation_1) + / (rotation_1_b1 * rotation_1_b1 + rotation_1_b2 * rotation_1_b2 + Eps); + + 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); + } + //=================================================================================================// + + 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]); + } + //=================================================================================================// + + 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(); + } + Vec3d getBinormalFromDeformationGradientTensor(const Mat3d &F) + { + return (F.col(2).cross(F.col(0))).normalized(); + } + //=================================================================================================// + Vecd getLinearVariableJump(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value) + { + return particle_i_value - particle_j_value + - 0.5 * r_ij * (gradient_particle_i_value + gradient_particle_j_value) * e_ij; + } + //=================================================================================================// + Vecd getWENOVariableJump(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value) + { + return getWENOLeftState(e_ij, r_ij, particle_i_value, + gradient_particle_i_value, particle_j_value, gradient_particle_j_value) + - getWENORightState(e_ij, r_ij, particle_i_value, + gradient_particle_i_value, particle_j_value, gradient_particle_j_value); + } + //=================================================================================================// + Vecd getWENOStateWithStencilPoints(const Vecd &v1, const Vecd &v2, const Vecd &v3, const Vecd &v4) + { + Vecd f1 = 0.5 * v2 + 0.5 * v3; + Vecd f2 = -0.5 * v1 + 1.5 * v2; + Vecd f3 = v2 / 3.0 + 5.0 * v3 / 6.0 - v4 / 6.0; + + Real epsilon = 1.0e-6; + Real s1 = (v2 - v3).dot(v2 - v3) + epsilon; + Real s2 = (v2 - v1).dot(v2 - v1) + epsilon; + Real s3 = (3.0 * v2 - 4.0 * v3 + v4).dot(3.0 * v2 - 4.0 * v3 + v4) / 4.0 + + 13.0 * (v2 - 2.0 * v3 + v4).dot(v2 - 2.0 * v3 + v4) / 12.0 + epsilon; + Real s12 = 13.0 * (v1 - 2.0 * v2 + v3).dot(v1 - 2.0 * v2 + v3) / 12.0 + + (v1 - v3).dot(v1 - v3) / 4.0 + epsilon; + Real s4 = (v1.dot(6649.0 * v1 - 30414.0 * v2 + 23094.0 * v3 - 5978.0 * v4) + + 3.0 * v2.dot(13667.0 * v2 - 23534.0 * v3 + 6338.0 * v4) + + 3.0 * v3.dot(11147.0 * v3 - 6458.0 * v4) + + 3169.0 * v4.dot(v4)) / 2880.0; + Real tau_4 = s4 - 0.5 * (s1 + s2); + + Real alpha_1 = (1.0 + (tau_4 / s1) * (tau_4 / s12)) / 3.0; + Real alpha_2 = (1.0 + (tau_4 / s2) * (tau_4 / s12)) / 6.0; + Real alpha_3 = (1.0 + tau_4 / s3) / 2.0; + Real w_1 = alpha_1 / (alpha_1 + alpha_2 + alpha_3); + Real w_2 = alpha_2 / (alpha_1 + alpha_2 + alpha_3); + Real w_3 = alpha_3 / (alpha_1 + alpha_2 + alpha_3); + + return w_1 * f1 + w_2 * f2 + w_3 * f3; + } + //=================================================================================================// + Vecd getWENOLeftState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value) + { + Vecd v1 = particle_i_value + gradient_particle_i_value * e_ij * r_ij; + Vecd v2 = particle_i_value; + Vecd v3 = particle_j_value; + Vecd v4 = particle_j_value - gradient_particle_j_value * e_ij * r_ij; + + return getWENOStateWithStencilPoints(v1, v2, v3, v4); + } + //=================================================================================================// + Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value) + { + Vecd v1 = particle_j_value - gradient_particle_j_value * e_ij * r_ij; + Vecd v2 = particle_j_value; + Vecd v3 = particle_i_value; + Vecd v4 = particle_i_value + gradient_particle_i_value * e_ij * r_ij; + + 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(); + Vec3d local_pseuodo_n_jump = transformation_matrix * pseudo_n_jump; + local_rotation_jump[0] = local_pseuodo_n_jump[0]; + 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_) + { + Mat3d corrected_almansi_strain = current_local_almansi_strain; + 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) + { + Mat3d correction_matrix = Mat3d::Zero(); + correction_matrix.block<2,2>(0,0) = local_deformation_part_one.block<2,2>(0,0).inverse(); + return correction_matrix; + } + + Mat3d getCorrectionMatrix_beam(const Mat3d &local_deformation_part_one) + { + Mat3d correction_matrix = Mat3d::Zero(); + correction_matrix.block<1, 1>(0, 0) = local_deformation_part_one.block<1, 1>(0, 0).inverse(); + return correction_matrix; + } + //=================================================================================================// + } +} diff --git a/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.h b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.h new file mode 100644 index 0000000000..9e65aaeda6 --- /dev/null +++ b/src/for_3D_build/particle_dynamics/solid_dynamics/slender_structure_math.h @@ -0,0 +1,95 @@ +/* -------------------------------------------------------------------------* + * SPHinXsys * + * -------------------------------------------------------------------------* + * SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle* + * Hydrodynamics for industrial compleX systems. It provides C++ APIs for * + * physical accurate simulation and aims to model coupled industrial dynamic* + * systems including fluid, solid, multi-body dynamics and beyond with SPH * + * (smoothed particle hydrodynamics), a meshless computational method using * + * particle discretization. * + * * + * SPHinXsys is partially funded by German Research Foundation * + * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * + * HU1527/12-1 and HU1527/12-4 * + * * + * Portions copyright (c) 2017-2022 Technical University of Munich and * + * the authors' affiliations. * + * * + * Licensed under the Apache License, Version 2.0 (the "License"); you may * + * not use this file except in compliance with the License. You may obtain a* + * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. * + * * + * ------------------------------------------------------------------------*/ +/** + * @file thin_structure_dynamics.h + * @brief Here, we define the math operation for thin structure dynamics. + * @details We consider here a weakly compressible solids. + * @author Dong Wu, Chi Zhang and Xiangyu Hu + */ + +#ifndef SLENDER_STRUCTURE_MATH_H +#define SLENDER_STRUCTURE_MATH_H + +#include "all_particle_dynamics.h" +#include "elastic_solid.h" +#include "weakly_compressible_fluid.h" +#include "base_kernel.h" + +namespace SPH +{ + namespace slender_structure_dynamics + { + /** + * Each of these basic vector rotations appears counterclockwise + * 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); + + /** get variable jump form gradient tensor. */ + Vecd getLinearVariableJump(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value); + Vecd getWENOVariableJump(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value); + + /** predict mid-point value by applying WENO reconstruction. */ + Vecd getWENOStateWithStencilPoints(const Vecd &v1, const Vecd &v2, const Vecd &v3, const Vecd &v4); + Vecd getWENOLeftState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + const Matd &gradient_particle_i_value, const Vecd &particle_j_value, const Matd &gradient_particle_j_value); + Vecd getWENORightState(const Vecd &e_ij, const Real &r_ij, const Vecd &particle_i_value, + 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); + } +} +#endif //THIN_STRUCTURE_MATH_H \ No newline at end of file diff --git a/src/shared/common/vector_functions.cpp b/src/shared/common/vector_functions.cpp index d382adf062..8bd1410e37 100644 --- a/src/shared/common/vector_functions.cpp +++ b/src/shared/common/vector_functions.cpp @@ -197,6 +197,31 @@ Mat3d getTransformationMatrix(const Vec3d &direction_of_z) return transformation_matrix; } //=================================================================================================// +//=================================================================================================// +Mat3d getTransformationMatrix(const Vec3d &direction_of_z, const Vec3d &direction_of_y) +{ + Mat3d transformation_matrix = Mat3d::Zero(); + Vec3d direction_of_x = direction_of_y.cross(direction_of_z); + transformation_matrix.row(0) = direction_of_x.transpose(); + transformation_matrix.row(1) = direction_of_y.transpose(); + transformation_matrix.row(2) = direction_of_z.transpose(); + + return transformation_matrix; +} + + +Mat2d getTransformationMatrix(const Vec2d &direction_of_z, const Vec2d &direction_of_y) +{ + Mat2d transformation_matrix = Mat2d::Zero(); + transformation_matrix.row(0) = direction_of_y.transpose(); + transformation_matrix.row(1) = direction_of_z.transpose(); + + return transformation_matrix; +} + + +//=================================================================================================// + Mat2d getDiagonal(const Mat2d &A) { Mat2d diag = Mat2d::Identity(); diff --git a/src/shared/common/vector_functions.h b/src/shared/common/vector_functions.h index 5e13068b9d..d7ef1c02d8 100644 --- a/src/shared/common/vector_functions.h +++ b/src/shared/common/vector_functions.h @@ -62,7 +62,8 @@ Real CalculateBiDotProduct(Mat3d Matrix1, Mat3d Matrix2); // calculate Real dot /** get transformation matrix. */ Mat2d getTransformationMatrix(const Vec2d &direction_of_y); Mat3d getTransformationMatrix(const Vec3d &direction_of_z); - +Mat3d getTransformationMatrix(const Vec3d &direction_of_z, const Vec3d &direction_of_y); +Mat2d getTransformationMatrix(const Vec2d &direction_of_z, const Vec2d &direction_of_y); /** get angle between two vectors. */ Real getCosineOfAngleBetweenTwoVectors(const Vec2d &vector_1, const Vec2d &vector_2); Real getCosineOfAngleBetweenTwoVectors(const Vec3d &vector_1, const Vec3d &vector_2); diff --git a/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h index 472b745979..76a3b00195 100644 --- a/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h +++ b/src/shared/particle_dynamics/solid_dynamics/all_solid_dynamics.h @@ -38,3 +38,4 @@ #include "loading_dynamics.h" #include "thin_structure_dynamics.h" #include "thin_structure_math.h" + diff --git a/src/shared/particle_generator/base_particle_generator.cpp b/src/shared/particle_generator/base_particle_generator.cpp index 3386fe3f8e..89edc1cd92 100644 --- a/src/shared/particle_generator/base_particle_generator.cpp +++ b/src/shared/particle_generator/base_particle_generator.cpp @@ -47,11 +47,27 @@ void SurfaceParticleGenerator::initializeSurfaceProperties(const Vecd &surface_n thickness_.push_back(thickness); } //=================================================================================================// + +LineParticleGenerator::LineParticleGenerator(SPHBody &sph_body) + : ParticleGenerator(sph_body), + n_(*base_particles_.getVariableByName("NormalDirection")), + thickness_(*base_particles_.getVariableByName("Thickness")), + b_n_(*base_particles_.getVariableByName("BinormalDirection")), + width_(*base_particles_.getVariableByName("Width")) {} +//=================================================================================================// +void LineParticleGenerator::initializeLineProperties(const Vecd &line_normal, const Vecd &line_binormal, Real thickness, Real width) +{ + n_.push_back(line_normal); + thickness_.push_back(thickness); + b_n_.push_back(line_binormal); + width_.push_back(width); +} +//=================================================================================================// void ObserverParticleGenerator::initializeGeometricVariables() { for (size_t i = 0; i < positions_.size(); ++i) { - initializePositionAndVolumetricMeasure(positions_[i], 0.0); + initializePositionAndVolumetricMeasure(positions_[i], 0.0); } } //=================================================================================================// diff --git a/src/shared/particle_generator/base_particle_generator.h b/src/shared/particle_generator/base_particle_generator.h index 37b79fd56d..2c8c49404d 100644 --- a/src/shared/particle_generator/base_particle_generator.h +++ b/src/shared/particle_generator/base_particle_generator.h @@ -97,6 +97,26 @@ class SurfaceParticleGenerator : public ParticleGenerator virtual void initializeSurfaceProperties(const Vecd &surface_normal, Real thickness); }; +/** + * @class LineParticleGenerator + * @brief Generate volumetric particles by initialize extra Line variables. + */ +class LineParticleGenerator : public ParticleGenerator +{ + public: + explicit LineParticleGenerator(SPHBody &sph_body); + virtual ~LineParticleGenerator(){}; + + protected: + StdLargeVec &n_; /**< line normal */ + StdLargeVec &thickness_; /**< line thickness */ + StdLargeVec &b_n_; /**< line binormal */ + StdLargeVec &width_; /**< line width */ + + /** Initialize line particle. */ + virtual void initializeLineProperties(const Vecd &line_normal, const Vecd &line_binormal, Real thickness, Real width); +}; + /** * @class ObserverParticleGenerator * @brief Generate particle directly from position-and-volume data. diff --git a/src/shared/particles/solid_particles.cpp b/src/shared/particles/solid_particles.cpp index 03ded438b9..058c6cbb86 100644 --- a/src/shared/particles/solid_particles.cpp +++ b/src/shared/particles/solid_particles.cpp @@ -10,14 +10,14 @@ namespace SPH { //=============================================================================================// SolidParticles::SolidParticles(SPHBody &sph_body, Solid *solid) - : BaseParticles(sph_body, solid), solid_(*solid) {} + : BaseParticles(sph_body, solid), solid_(*solid) {} //=================================================================================================// void SolidParticles::initializeOtherVariables() { - BaseParticles::initializeOtherVariables(); + BaseParticles::initializeOtherVariables(); registerVariable(pos0_, "InitialPosition", [&](size_t i) -> Vecd { return pos_[i]; }); - registerVariable(n_, "NormalDirection"); + registerVariable(n_, "NormalDirection"); registerVariable(n0_, "InitialNormalDirection", [&](size_t i) -> Vecd { return n_[i]; }); registerVariable(B_, "CorrectionMatrix", [&](size_t i) -> Matd @@ -25,244 +25,244 @@ void SolidParticles::initializeOtherVariables() } //=============================================================================================// ElasticSolidParticles:: - ElasticSolidParticles(SPHBody &sph_body, ElasticSolid *elastic_solid) - : SolidParticles(sph_body, elastic_solid), - elastic_solid_(*elastic_solid) {} + ElasticSolidParticles(SPHBody &sph_body, ElasticSolid *elastic_solid) + : SolidParticles(sph_body, elastic_solid), + elastic_solid_(*elastic_solid) {} //=================================================================================================// void ElasticSolidParticles::initializeOtherVariables() { - SolidParticles::initializeOtherVariables(); - /** - * register particle data - */ + SolidParticles::initializeOtherVariables(); + /** + * register particle data + */ registerVariable(F_, "DeformationGradient", [&](size_t i) -> Matd { return Matd::Identity(); }); - registerVariable(dF_dt_, "DeformationRate"); - /** - * register FSI data - */ - registerVariable(vel_ave_, "AverageVelocity"); - registerVariable(acc_ave_, "AverageAcceleration"); - /** - * add restart output particle data + registerVariable(dF_dt_, "DeformationRate"); + /** + * register FSI data + */ + registerVariable(vel_ave_, "AverageVelocity"); + registerVariable(acc_ave_, "AverageAcceleration"); + /** + * add restart output particle data + */ + addVariableToRestart("DeformationGradient"); + /** + * add restart output particle data */ - addVariableToRestart("DeformationGradient"); - /** - * add restart output particle data - */ - addVariableToWrite("NormalDirection"); - addDerivedVariableToWrite(); - addDerivedVariableToWrite(); - addDerivedVariableToWrite(); - addVariableToRestart("DeformationGradient"); - // get which stress measure is relevant for the material - stress_measure_ = elastic_solid_.getRelevantStressMeasureName(); + addVariableToWrite("NormalDirection"); + addDerivedVariableToWrite(); + addDerivedVariableToWrite(); + addDerivedVariableToWrite(); + addVariableToRestart("DeformationGradient"); + // get which stress measure is relevant for the material + stress_measure_ = elastic_solid_.getRelevantStressMeasureName(); } //=================================================================================================// Matd ElasticSolidParticles::getGreenLagrangeStrain(size_t particle_i) { - Matd F = F_[particle_i]; - return 0.5 * (F.transpose() * F - Matd::Identity()); + Matd F = F_[particle_i]; + return 0.5 * (F.transpose() * F - Matd::Identity()); } //=================================================================================================// Vecd ElasticSolidParticles::getPrincipalStrains(size_t particle_i) { Matd epsilon = getGreenLagrangeStrain(particle_i); - return getPrincipalValuesFromMatrix(epsilon); + return getPrincipalValuesFromMatrix(epsilon); } //=================================================================================================// Matd ElasticSolidParticles::getStressCauchy(size_t particle_i) { - Matd F = F_[particle_i]; - Matd stress_PK2 = elastic_solid_.StressPK2(F, particle_i); - return (1.0 / F.determinant()) * F * stress_PK2 * F.transpose(); + Matd F = F_[particle_i]; + Matd stress_PK2 = elastic_solid_.StressPK2(F, particle_i); + return (1.0 / F.determinant()) * F * stress_PK2 * F.transpose(); } //=================================================================================================// Matd ElasticSolidParticles::getStressPK2(size_t particle_i) { - return elastic_solid_.StressPK2(F_[particle_i], particle_i); + return elastic_solid_.StressPK2(F_[particle_i], particle_i); } //=================================================================================================// Vecd ElasticSolidParticles::getPrincipalStresses(size_t particle_i) { - Matd sigma; - if (stress_measure_ == "Cauchy") - { - sigma = getStressCauchy(particle_i); // Cauchy stress - } - else if (stress_measure_ == "PK2") - { - sigma = getStressPK2(particle_i); // Second Piola-Kirchhoff stress - } - else - { - throw std::runtime_error("get_Principal_stresses: wrong input"); - } + Matd sigma; + if (stress_measure_ == "Cauchy") + { + sigma = getStressCauchy(particle_i); // Cauchy stress + } + else if (stress_measure_ == "PK2") + { + sigma = getStressPK2(particle_i); // Second Piola-Kirchhoff stress + } + else + { + throw std::runtime_error("get_Principal_stresses: wrong input"); + } - return getPrincipalValuesFromMatrix(sigma); + return getPrincipalValuesFromMatrix(sigma); } //=================================================================================================// Real ElasticSolidParticles::getVonMisesStress(size_t particle_i) { - Matd sigma; - if (stress_measure_ == "Cauchy") - { - sigma = getStressCauchy(particle_i); // Cauchy stress - } - else if (stress_measure_ == "PK2") - { - sigma = getStressPK2(particle_i); // Second Piola-Kirchhoff stress - } - else - { - throw std::runtime_error("get_von_Mises_stress: wrong input"); - } + Matd sigma; + if (stress_measure_ == "Cauchy") + { + sigma = getStressCauchy(particle_i); // Cauchy stress + } + else if (stress_measure_ == "PK2") + { + sigma = getStressPK2(particle_i); // Second Piola-Kirchhoff stress + } + else + { + throw std::runtime_error("get_von_Mises_stress: wrong input"); + } - return getVonMisesStressFromMatrix(sigma); + return getVonMisesStressFromMatrix(sigma); } //=================================================================================================// StdLargeVec ElasticSolidParticles::getVonMisesStrainVector(std::string strain_measure) { - StdLargeVec strain_vector = {}; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - Real strain = 0.0; - if (strain_measure == "static") - { - strain = getVonMisesStrain(index_i); - } - else if (strain_measure == "dynamic") - { - strain = getVonMisesStrainDynamic(index_i, elastic_solid_.PoissonRatio()); - } - else - { - throw std::runtime_error("getVonMisesStrainVector: wrong input"); - } + StdLargeVec strain_vector = {}; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + Real strain = 0.0; + if (strain_measure == "static") + { + strain = getVonMisesStrain(index_i); + } + else if (strain_measure == "dynamic") + { + strain = getVonMisesStrainDynamic(index_i, elastic_solid_.PoissonRatio()); + } + else + { + throw std::runtime_error("getVonMisesStrainVector: wrong input"); + } - strain_vector.push_back(strain); - } + strain_vector.push_back(strain); + } - return strain_vector; + return strain_vector; } //=================================================================================================// Real ElasticSolidParticles::getVonMisesStrainMax(std::string strain_measure) { - Real strain_max = 0; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - Real strain = 0.0; - if (strain_measure == "static") - { - strain = getVonMisesStrain(index_i); - } - else if (strain_measure == "dynamic") - { - strain = getVonMisesStrainDynamic(index_i, elastic_solid_.PoissonRatio()); - } - else - { - throw std::runtime_error("getVonMisesStrainMax: wrong input"); - } - if (strain_max < strain) - strain_max = strain; - } - return strain_max; + Real strain_max = 0; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + Real strain = 0.0; + if (strain_measure == "static") + { + strain = getVonMisesStrain(index_i); + } + else if (strain_measure == "dynamic") + { + strain = getVonMisesStrainDynamic(index_i, elastic_solid_.PoissonRatio()); + } + else + { + throw std::runtime_error("getVonMisesStrainMax: wrong input"); + } + if (strain_max < strain) + strain_max = strain; + } + return strain_max; } //=================================================================================================// Real ElasticSolidParticles::getPrincipalStressMax() { - Real stress_max = 0.0; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - /** take the max. component, which is the first one, this represents the max. tension. */ + Real stress_max = 0.0; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + /** take the max. component, which is the first one, this represents the max. tension. */ Real stress = getPrincipalStresses(index_i)[0]; - if (stress_max < stress) - stress_max = stress; - } - return stress_max; + if (stress_max < stress) + stress_max = stress; + } + return stress_max; }; //=================================================================================================// StdLargeVec ElasticSolidParticles::getVonMisesStressVector() { - StdLargeVec stress_vector = {}; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - Real stress = getVonMisesStress(index_i); - stress_vector.push_back(stress); - } - return stress_vector; + StdLargeVec stress_vector = {}; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + Real stress = getVonMisesStress(index_i); + stress_vector.push_back(stress); + } + return stress_vector; } //=================================================================================================// Real ElasticSolidParticles::getVonMisesStressMax() { - Real stress_max = 0.0; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - Real stress = getVonMisesStress(index_i); - if (stress_max < stress) - stress_max = stress; - } - return stress_max; + Real stress_max = 0.0; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + Real stress = getVonMisesStress(index_i); + if (stress_max < stress) + stress_max = stress; + } + return stress_max; } //=================================================================================================// Vecd ElasticSolidParticles::displacement(size_t particle_i) { - return pos_[particle_i] - pos0_[particle_i]; + return pos_[particle_i] - pos0_[particle_i]; } //=================================================================================================// Vecd ElasticSolidParticles::normal(size_t particle_i) { - return n_[particle_i]; + return n_[particle_i]; } //=================================================================================================// StdLargeVec ElasticSolidParticles::getDisplacement() { - StdLargeVec displacement_vector = {}; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - displacement_vector.push_back(displacement(index_i)); - } - return displacement_vector; + StdLargeVec displacement_vector = {}; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + displacement_vector.push_back(displacement(index_i)); + } + return displacement_vector; } //=================================================================================================// Real ElasticSolidParticles::getMaxDisplacement() { - Real displ_max = 0.0; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - Real displ = displacement(index_i).norm(); - if (displ_max < displ) - displ_max = displ; - } - return displ_max; + Real displ_max = 0.0; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + Real displ = displacement(index_i).norm(); + if (displ_max < displ) + displ_max = displ; + } + return displ_max; }; //=================================================================================================// StdLargeVec ElasticSolidParticles::getNormal() { - StdLargeVec normal_vector = {}; - for (size_t index_i = 0; index_i < pos0_.size(); index_i++) - { - normal_vector.push_back(normal(index_i)); - } - return normal_vector; + StdLargeVec normal_vector = {}; + for (size_t index_i = 0; index_i < pos0_.size(); index_i++) + { + normal_vector.push_back(normal(index_i)); + } + return normal_vector; } //=============================================================================================// ShellParticles::ShellParticles(SPHBody &sph_body, ElasticSolid *elastic_solid) - : ElasticSolidParticles(sph_body, elastic_solid), thickness_ref_(1.0) + : ElasticSolidParticles(sph_body, elastic_solid), thickness_ref_(1.0) { - //---------------------------------------------------------------------- - // modify kernel function for surface particles - //---------------------------------------------------------------------- - sph_body.sph_adaptation_->getKernel()->reduceOnce(); - //---------------------------------------------------------------------- - // register geometric data only - //---------------------------------------------------------------------- - registerVariable(n_, "NormalDirection"); - registerVariable(thickness_, "Thickness"); - /** - * add particle reload data - */ + //---------------------------------------------------------------------- + // modify kernel function for surface particles + //---------------------------------------------------------------------- + //sph_body.sph_adaptation_->getKernel()->reduceOnce(); + //---------------------------------------------------------------------- + // register geometric data only + //---------------------------------------------------------------------- + registerVariable(n_, "NormalDirection"); + registerVariable(thickness_, "Thickness"); + /** + * add particle reload data + */ addVariableToList(variables_to_reload_, "NormalDirection"); addVariableToList(variables_to_reload_, "Thickness"); } @@ -329,10 +329,75 @@ void ShellParticles::initializeOtherVariables() */ for (size_t i = 0; i != real_particles_bound_; ++i) { - transformation_matrix_[i] = getTransformationMatrix(n_[i]); - numerical_damping_scaling_[i](Dimensions - 1, Dimensions - 1) = - thickness_[i] < sph_body_.sph_adaptation_->ReferenceSmoothingLength() ? thickness_[i] : sph_body_.sph_adaptation_->ReferenceSmoothingLength(); + transformation_matrix_[i] = getTransformationMatrix(n_[i]); + numerical_damping_scaling_[i](Dimensions - 1, Dimensions - 1) = + thickness_[i] < sph_body_.sph_adaptation_->ReferenceSmoothingLength() ? thickness_[i] : sph_body_.sph_adaptation_->ReferenceSmoothingLength(); } } -//=================================================================================================// -} // namespace SPH + BarParticles::BarParticles(SPHBody &sph_body, ElasticSolid *elastic_solid) + : ShellParticles(sph_body, elastic_solid), width_ref_(1.0) + { + //---------------------------------------------------------------------- + // modify kernel function for surface particles + //---------------------------------------------------------------------- + // sph_body.sph_adaptation_->getKernel()->reduceTwice(); + //---------------------------------------------------------------------- + // register geometric data only + //---------------------------------------------------------------------- + registerVariable(b_n_, "BinormalDirection"); + registerVariable(width_, "Width"); + /** + * add particle reload data + */ + //addVariableNameToList(variables_to_reload_, "BinormalDirection"); + //addVariableNameToList(variables_to_reload_, "Width"); + } + //=================================================================================================// + void BarParticles::initializeOtherVariables() + { + ShellParticles::initializeOtherVariables(); + /** + * register particle data + */ + registerVariable(b_n0_, "InitialBinormalDirection", + [&](size_t i) -> Vecd + { return b_n_[i]; }); + registerVariable(pseudo_b_n_, "PseudoBinormal", + [&](size_t i) -> Vecd + { return pseudo_b_n_[i]; }); + registerVariable(dpseudo_b_n_dt_, "PseudoBinormalChangeRate"); + registerVariable(dpseudo_b_n_d2t_, "PseudoBinormal2ndOrderTimeDerivative"); + registerVariable(rotation_b_, "Rotation_b"); + registerVariable(angular_b_vel_, "AngularVelocity_b"); + registerVariable(dangular_b_vel_dt_, "AngularAccelerationofBinormal"); + registerVariable(F_b_bending_, "b_BendingDeformationGradient"); + registerVariable(dF_b_bending_dt_, "b_BendingDeformationGradientChangeRate"); + + registerVariable(global_b_shear_stress_, "Global_b_ShearStress"); + registerVariable(global_b_stress_, "Global_b_Stress"); + registerVariable(global_b_moment_, "Global_b_Moment"); + + /** + * for rotation. + */ + + addVariableToRestart("PseudoBinormal"); + addVariableToRestart("Rotation_b"); + addVariableToRestart("AngularVelocity_b"); + /** + * add basic output particle data + */ + addVariableToWrite("BinormalDirection"); + addVariableToWrite("Rotation_b"); + /** + * initialize transformation matrix + */ + for (size_t i = 0; i != real_particles_bound_; ++i) + { + + transformation_matrix_[i] = getTransformationMatrix(n_[i], b_n_[i]); + + } + } + } + \ No newline at end of file diff --git a/src/shared/particles/solid_particles.h b/src/shared/particles/solid_particles.h index 192957532e..ba7df7041d 100644 --- a/src/shared/particles/solid_particles.h +++ b/src/shared/particles/solid_particles.h @@ -10,9 +10,9 @@ * * * SPHinXsys is partially funded by German Research Foundation * * (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, * - * HU1527/12-1 and HU1527/12-4. * + * HU1527/12-1 and HU1527/12-4 * * * - * Portions copyright (c) 2017-2023 Technical University of Munich and * + * Portions copyright (c) 2017-2022 Technical University of Munich and * * the authors' affiliations. * * * * Licensed under the Apache License, Version 2.0 (the "License"); you may * @@ -195,5 +195,43 @@ class ShellParticles : public ElasticSolidParticles /** Return this pointer. */ virtual ShellParticles *ThisObjectPtr() override { return this; }; }; -} // namespace SPH + + class BarParticles : public ShellParticles + { + public: + BarParticles(SPHBody &sph_body, ElasticSolid *elastic_solid); + virtual ~BarParticles(){}; + + Real width_ref_; + StdLargeVec b_n_; /**< binormal direction */ + StdLargeVec width_; + StdLargeVec b_n0_; /**< initial binormal direction */ + + StdLargeVec pseudo_b_n_; /**< current pseudo-binormal vector */ + StdLargeVec dpseudo_b_n_dt_; /**< pseudo-binormal vector change rate */ + StdLargeVec dpseudo_b_n_d2t_; /**< pseudo-binormal vector second order time derivation */ + StdLargeVec global_b_shear_stress_; /**< global b shear stress */ + StdLargeVec global_b_stress_; /**< global b stress for pair interaction */ + StdLargeVec global_b_moment_; /**< global b bending moment for pair interaction */ + //---------------------------------------------------------------------- + // extra generalized coordinate and velocity in local coordinate + //---------------------------------------------------------------------- + StdLargeVec rotation_b_; /**< rotation angle of the initial binormal respective to each axis */ + StdLargeVec angular_b_vel_; /**< angular velocity respective to each axis */ + StdLargeVec dangular_b_vel_dt_; /**< angular acceleration of respective to each axis*/ + // extra deformation and deformation rate in local coordinate + //---------------------------------------------------------------------- + StdLargeVec F_b_bending_; /**< bending deformation gradient */ + StdLargeVec dF_b_bending_dt_; /**< bending deformation gradient change rate */ + //---------------------------------------------------------------------- + /** get particle volume. */ + virtual Real ParticleVolume(size_t index_i) override { return Vol_[index_i] * thickness_[index_i]*width_[index_i]; } + /** get particle mass. */ + virtual Real ParticleMass(size_t index_i) override { return mass_[index_i] * thickness_[index_i] * width_[index_i]; } + /** Initialize variable for shell particles. */ + virtual void initializeOtherVariables() override; + virtual BarParticles *ThisObjectPtr() override{return this;}; + }; + +} #endif // SOLID_PARTICLES_H diff --git a/tests/3d_examples/test_3d_slender_beam/CMakeLists.txt b/tests/3d_examples/test_3d_slender_beam/CMakeLists.txt new file mode 100644 index 0000000000..7ad4e1cf8a --- /dev/null +++ b/tests/3d_examples/test_3d_slender_beam/CMakeLists.txt @@ -0,0 +1,17 @@ +STRING(REGEX REPLACE ".*/(.*)" "\\1" CURRENT_FOLDER ${CMAKE_CURRENT_SOURCE_DIR}) +PROJECT("${CURRENT_FOLDER}") + +SET(LIBRARY_OUTPUT_PATH ${PROJECT_BINARY_DIR}/lib) +SET(EXECUTABLE_OUTPUT_PATH "${PROJECT_BINARY_DIR}/bin/") +SET(BUILD_INPUT_PATH "${EXECUTABLE_OUTPUT_PATH}/input") +SET(BUILD_RELOAD_PATH "${EXECUTABLE_OUTPUT_PATH}/reload") + +add_executable(${PROJECT_NAME}) +aux_source_directory(. DIR_SRCS) +target_sources(${PROJECT_NAME} PRIVATE ${DIR_SRCS}) +target_link_libraries(${PROJECT_NAME} sphinxsys_3d) +set_target_properties(${PROJECT_NAME} PROPERTIES VS_DEBUGGER_WORKING_DIRECTORY "${EXECUTABLE_OUTPUT_PATH}") + +add_test(NAME ${PROJECT_NAME} + COMMAND ${PROJECT_NAME} + WORKING_DIRECTORY ${EXECUTABLE_OUTPUT_PATH}) diff --git a/tests/3d_examples/test_3d_slender_beam/test_3d_slender_beam.cpp b/tests/3d_examples/test_3d_slender_beam/test_3d_slender_beam.cpp new file mode 100644 index 0000000000..481dc42610 --- /dev/null +++ b/tests/3d_examples/test_3d_slender_beam/test_3d_slender_beam.cpp @@ -0,0 +1,262 @@ +/** + * @file test_3d_slender_beam.cpp + * @brief This is the benchmark test of the bar. + * @details We consider the body force applied on a quasi-static stright beam. + * @author Xipeng Lyu + * @ref doi.org/10.1016/j.ijnonlinmec.2014.04.009, doi.org/10.1201/9780849384165 + */ +#include "sphinxsys.h" +#include +#include "slender_structure_dynamics.h" +#include "slender_structure_math.h" +using namespace SPH; + + +/** + * @brief Basic geometry parameters and numerical setup. + */ +Real PL = 10.0; /** Length of the square plate. */ +Real PH = 10.0; /** Width of the square plate. */ +Real PT = 1.0; /** Thickness of the square plate. */ +Real PW = 1.0; /** Thickness of the square plate. */ +Vec3d n_0 = Vec3d(0.0, 0.0, 1.0); /** Pseudo-normal. */ +Vec3d b_n_0 = Vec3d(0.0, 1.0, 0.0); /** Pseudo-normal. */ +int particle_number = 40; /** Particle number in the direction of the length */ +Real resolution_ref = PL / (Real)particle_number; /** Initial reference particle spacing. */ +int BWD = 1; /** Width of the boundary layer measured by number of particles. */ +Real BW = resolution_ref * (Real)BWD; /** Boundary width, determined by specific layer of boundary particles. */ +/** Domain bounds of the system. */ +BoundingBox system_domain_bounds(Vec3d(-BW, -BW, -0.5 * resolution_ref), + Vec3d(PL + BW, PH + BW, 0.5 * resolution_ref)); +// Observer location +StdVec observation_location = {Vecd(0.5 * PL, 0.0, 0.0)}; + +/** For material properties of the solid. */ +Real rho0_s = 1.0; /** Normalized density. */ +Real Youngs_modulus = 1.3024653e6; /** Normalized Youngs Modulus. */ +Real poisson = 0.3; /** Poisson ratio. */ +Real physical_viscosity = 200.0; /** physical damping, here we choose the same value as numerical viscosity. */ + +Real q = 100.0 * Youngs_modulus * 1.0e-4; /** Total distributed load. */ +Real time_to_full_external_force = 0.1; + +Real gravitational_acceleration = 0.009646; + +Real observed_quantity_0(0.0); +Real observed_quantity_n(0.0); +Real displ_max_reference = 1.8687; +TEST(Beam, MaxDisplacement) +{ + Real displ_max = observed_quantity_n - observed_quantity_0; + EXPECT_NEAR(displ_max, displ_max_reference, displ_max_reference * 0.1); + std::cout << "displ_max: " << displ_max << std::endl; +} + +/** Define application dependent particle generator for thin structure. */ +class BarParticleGenerator : public LineParticleGenerator +{ +public: + explicit BarParticleGenerator(SPHBody &sph_body) : LineParticleGenerator(sph_body){ + sph_body.sph_adaptation_->getKernel()->reduceTwice(); + }; + virtual void initializeGeometricVariables() override + { + // the beam and boundary + for (int i = 0; i < (particle_number + 2 * BWD); i++) + { + + Real x = resolution_ref * i - BW + resolution_ref * 0.5; + Real y = 0.0; + Real z = 0.0; + initializePositionAndVolumetricMeasure(Vecd(x, 0.0, 0.0), resolution_ref ); + initializeLineProperties(n_0,b_n_0, PT,PW); + + } + } +}; +/** Define the boundary geometry. */ +class BoundaryGeometryParallelToXAxis : public BodyPartByParticle +{ +public: + BoundaryGeometryParallelToXAxis(SPHBody &body, const std::string &body_part_name) + : BodyPartByParticle(body, body_part_name) + { + TaggingParticleMethod tagging_particle_method = std::bind(&BoundaryGeometryParallelToXAxis::tagManually, this, _1); + tagParticles(tagging_particle_method); + }; + virtual ~BoundaryGeometryParallelToXAxis() {}; + +private: + void tagManually(size_t index_i) + { + if (base_particles_.pos_[index_i][1] < 0.0 || base_particles_.pos_[index_i][1] > PH) + { + body_part_particles_.push_back(index_i); + } + }; +}; +class BoundaryGeometryParallelToYAxis : public BodyPartByParticle +{ +public: + BoundaryGeometryParallelToYAxis(SPHBody &body, const std::string &body_part_name) + : BodyPartByParticle(body, body_part_name) + { + TaggingParticleMethod tagging_particle_method = std::bind(&BoundaryGeometryParallelToYAxis::tagManually, this, _1); + tagParticles(tagging_particle_method); + }; + virtual ~BoundaryGeometryParallelToYAxis() {}; + +private: + void tagManually(size_t index_i) + { + if (base_particles_.pos_[index_i][0] < 0.0 || base_particles_.pos_[index_i][0] > PL) + { + body_part_particles_.push_back(index_i); + } + }; +}; +/** + * define time dependent external force + */ +class TimeDependentExternalForce : public Gravity +{ +public: + explicit TimeDependentExternalForce(Vecd external_force) + : Gravity(external_force) {} + virtual Vecd InducedAcceleration(Vecd &position) override + { + Real current_time = GlobalStaticVariables::physical_time_; + return current_time < time_to_full_external_force + ? current_time * global_acceleration_ / time_to_full_external_force + : global_acceleration_; + } +}; +/** + * The main program + */ +int main(int ac, char *av[]) +{ + /** Setup the system. */ + SPHSystem system(system_domain_bounds, resolution_ref); + + /** create a bar body. */ + SolidBody bar_body(system, makeShared("BarBody")); + bar_body.defineParticlesAndMaterial(rho0_s, Youngs_modulus, poisson); + bar_body.generateParticles(); + + /** Define Observer. */ + ObserverBody bar_observer(system, "BarObserver"); + bar_observer.defineParticlesAndMaterial(); + bar_observer.generateParticles(observation_location); + + /** Set body contact map + * The contact map gives the data connections between the bodies + * basically the the range of bodies to build neighbor particle lists + */ + InnerRelation bar_body_inner(bar_body); + ContactRelation bar_observer_contact(bar_observer, {&bar_body}); + + /** Common particle dynamics. */ + SimpleDynamics initialize_external_force(bar_body, + makeShared(Vec3d(0.0, 0.0, q / (PT * rho0_s) - gravitational_acceleration))); + + /** + * This section define all numerical methods will be used in this case. + */ + /** Corrected configuration. */ + InteractionDynamics + corrected_configuration(bar_body_inner); + /** Time step size calculation. */ + ReduceDynamics computing_time_step_size(bar_body); + /** active-passive stress relaxation. */ + Dynamics1Level + stress_relaxation_first_half(bar_body_inner); + Dynamics1Level + stress_relaxation_second_half(bar_body_inner); + /** Constrain the Boundary. */ + BoundaryGeometryParallelToXAxis boundary_geometry_x(bar_body, "BoundaryGeometryParallelToXAxis"); + SimpleDynamics + constrain_holder_x(boundary_geometry_x, 0); + BoundaryGeometryParallelToYAxis boundary_geometry_y(bar_body, "BoundaryGeometryParallelToYAxis"); + SimpleDynamics + constrain_holder_y(boundary_geometry_y, 1); + DampingWithRandomChoice>> + bar_position_damping(0.5, bar_body_inner, "Velocity", physical_viscosity); + DampingWithRandomChoice>> + bar_rotation_damping(0.5, bar_body_inner, "AngularVelocity", physical_viscosity); + DampingWithRandomChoice>> + bar_rotation_b_damping(0.5, bar_body_inner, "AngularVelocity_b", physical_viscosity); + /** Output */ + IOEnvironment io_environment(system); + BodyStatesRecordingToVtp write_states(io_environment, system.real_bodies_); + ObservedQuantityRecording write_beam_max_displacement("Position", io_environment, bar_observer_contact); + + /** Apply initial condition. */ + system.initializeSystemCellLinkedLists(); + system.initializeSystemConfigurations(); + corrected_configuration.exec(); + + /** + * From here the time stepping begins. + * Set the starting time. + */ + GlobalStaticVariables::physical_time_ = 0.0; + write_states.writeToFile(0); + write_beam_max_displacement.writeToFile(0); + observed_quantity_0 = (*write_beam_max_displacement.getObservedQuantity())[0][2]; + + /** Setup physical parameters. */ + int ite = 0; + Real end_time = 0.8; + Real output_period = end_time / 100.0; + Real dt = 0.0; + /** Statistics for computing time. */ + TickCount t1 = TickCount::now(); + TimeInterval interval; + /** + * Main loop + */ + while (GlobalStaticVariables::physical_time_ < end_time) + { + Real integral_time = 0.0; + while (integral_time < output_period) + { + if (ite % 100 == 0) + { + std::cout << "N=" << ite << " Time: " + << GlobalStaticVariables::physical_time_ << " dt: " + << dt << "\n"; + } + initialize_external_force.exec(dt); + stress_relaxation_first_half.exec(dt); + constrain_holder_x.exec(dt); + constrain_holder_y.exec(dt); + bar_position_damping.exec(dt); + bar_rotation_damping.exec(dt); + bar_rotation_b_damping.exec(dt); + constrain_holder_x.exec(dt); + constrain_holder_y.exec(dt); + stress_relaxation_second_half.exec(dt); + + ite++; + dt = computing_time_step_size.exec(); + integral_time += dt; + GlobalStaticVariables::physical_time_ += dt; + } + write_beam_max_displacement.writeToFile(ite); + TickCount t2 = TickCount::now(); + write_states.writeToFile(); + TickCount t3 = TickCount::now(); + interval += t3 - t2; + } + TickCount t4 = TickCount::now(); + + TimeInterval tt; + tt = t4 - t1 - interval; + std::cout << "Total wall time for computation: " << tt.seconds() << " seconds." << std::endl; + + observed_quantity_n = (*write_beam_max_displacement.getObservedQuantity())[0][2]; + + testing::InitGoogleTest(&ac, av); + return RUN_ALL_TESTS(); +}