From 98707e10e7850001ac3a7ec3d3d5b7f60740210f Mon Sep 17 00:00:00 2001 From: Xiangyu Hu Date: Mon, 8 Apr 2024 16:18:19 +0200 Subject: [PATCH 1/2] all compiled --- .../structural_simulation_class.cpp | 4 +- .../structural_simulation_class.h | 4 +- .../continuum_integration.h | 6 +- .../general_dynamics/general_constraint.h | 71 ++ .../solid_dynamics/constraint_dynamics.cpp | 9 +- .../solid_dynamics/constraint_dynamics.h | 108 +-- src/shared/particles/base_particles.h | 5 +- src/shared/particles/base_particles.hpp | 7 +- .../test_2d_anisotropic_beam.cpp | 2 +- .../test_2d_elastic_gate/elastic_gate.cpp | 2 +- tests/2d_examples/test_2d_fsi2/fsi2.cpp | 2 +- .../hydrostatic_fsi.cpp | 2 +- .../nonlinear_wave_fsi.h | 888 +++++++++--------- .../oscillating_beam.cpp | 2 +- .../oscillating_beam_UL.cpp | 2 +- .../cauchy_oscillating_beam.cpp | 2 +- tests/2d_examples/test_2d_owsc/owsc.h | 12 +- .../test_2d_self_contact/self_contact.cpp | 2 +- .../shell_beam_collision.cpp | 2 +- .../test_2d_stretching/stretching.cpp | 12 +- .../viscous_cream_drop.cpp | 2 +- .../beam_pulling_pressure_load.cpp | 2 +- .../test_3d_dynamic_plate.cpp | 2 +- .../excitation-contraction.cpp | 2 +- .../excitation_contraction.cpp | 2 +- .../src/muscle_activation.cpp | 2 +- .../muscle_soft_body_contact.cpp | 2 +- .../muscle_solid_contact.cpp | 2 +- .../nonlinear_wave_fsi.h | 12 +- .../oscillating_plate_UL.cpp | 5 +- .../passive_cantilever.cpp | 2 +- .../passive_cantilever_neohookean.cpp | 2 +- .../pkj_lv_electrocontraction.cpp | 2 +- tests/3d_examples/test_3d_roof/3d_roof.cpp | 2 +- .../test_3d_roof_parametric_cvt.cpp | 2 +- .../twisting_column.cpp | 2 +- .../lid_driven_cavity.cpp | 4 +- .../passive_cantilever_LG.cpp | 2 +- .../2d_velocity_gradient.cpp | 4 +- 39 files changed, 610 insertions(+), 589 deletions(-) diff --git a/modules/structural_simulation/structural_simulation_class.cpp b/modules/structural_simulation/structural_simulation_class.cpp index 01bd6ec116..85e4a57198 100755 --- a/modules/structural_simulation/structural_simulation_class.cpp +++ b/modules/structural_simulation/structural_simulation_class.cpp @@ -503,7 +503,7 @@ void StructuralSimulation::initializeConstrainSolidBody() for (size_t i = 0; i < body_indices_fixed_constraint_.size(); i++) { int body_index = body_indices_fixed_constraint_[i]; - fixed_constraint_body_.emplace_back(makeShared>(*solid_body_list_[body_index]->getSolidBodyFromMesh())); + fixed_constraint_body_.emplace_back(makeShared>(*solid_body_list_[body_index]->getSolidBodyFromMesh())); } } @@ -527,7 +527,7 @@ void StructuralSimulation::initializeConstrainSolidBodyRegion() // create the triangle mesh of the box BodyPartFromMesh *bp = body_part_tri_mesh_ptr_keeper_.createPtr( *solid_body_list_[body_index]->getSolidBodyFromMesh(), makeShared(halfsize_bbox, resolution, center, imported_stl_list_[body_index])); - fixed_constraint_region_.emplace_back(makeShared>(*bp)); + fixed_constraint_region_.emplace_back(makeShared>(*bp)); } } diff --git a/modules/structural_simulation/structural_simulation_class.h b/modules/structural_simulation/structural_simulation_class.h index a58b4232a4..c679eeb4ee 100755 --- a/modules/structural_simulation/structural_simulation_class.h +++ b/modules/structural_simulation/structural_simulation_class.h @@ -228,10 +228,10 @@ class StructuralSimulation StdVec>> surface_spring_; StdVec surface_spring_tuple_; // for ConstrainSolidBody - StdVec>> fixed_constraint_body_; + StdVec>> fixed_constraint_body_; StdVec body_indices_fixed_constraint_; // for ConstrainSolidBodyRegion - StdVec>> fixed_constraint_region_; + StdVec>> fixed_constraint_region_; StdVec body_indices_fixed_constraint_region_; // for PositionSolidBody StdVec>> position_solid_body_; diff --git a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h index 9147140539..3e1a2fc5c4 100644 --- a/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h +++ b/src/shared/particle_dynamics/continuum_dynamics/continuum_integration.h @@ -23,7 +23,7 @@ /** * @file continuum_integration.h * @brief Here, we define the algorithm classes for continuum dynamics within the body. - * @details We consider here weakly compressible assumption to model elastic and + * @details We consider here weakly compressible assumption to model elastic and * plastic materials with the updated Lagrangian framework. * @author Shuaihao Zhang and Xiangyu Hu */ @@ -99,10 +99,6 @@ class ShearStressRelaxation : public fluid_dynamics::BaseIntegration &B_; }; -using FixBodyPartConstraint = solid_dynamics::FixConstraint; -using FixedInAxisDirection = solid_dynamics::BaseFixedInAxisDirection; -using ConstrainSolidBodyMassCenter = solid_dynamics::BaseConstrainSolidBodyMassCenter; - template class BasePlasticIntegration : public fluid_dynamics::BaseIntegration { diff --git a/src/shared/particle_dynamics/general_dynamics/general_constraint.h b/src/shared/particle_dynamics/general_dynamics/general_constraint.h index 1bf547163f..588f3c51ea 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_constraint.h +++ b/src/shared/particle_dynamics/general_dynamics/general_constraint.h @@ -55,5 +55,76 @@ class ShapeSurfaceBounding : public BaseLocalDynamics, LevelSetShape *level_set_shape_; Real constrained_distance_; }; + +/** + * @class MotionConstraint + * @brief Base class for constraining with prescribed motion. + * Exact motion function will be defined in derive class. + * Note that, we do not impose acceleration, so that this constraint + * must be imposed after updating particle velocity by forces + * and before updating particle position. + */ +template +class MotionConstraint : public BaseLocalDynamics, public GeneralDataDelegateSimple +{ + public: + explicit MotionConstraint(DynamicsIdentifier &identifier) + : BaseLocalDynamics(identifier), + GeneralDataDelegateSimple(identifier.getSPHBody()), + pos_(this->particles_->pos_), + pos0_(*this->particles_->registerSharedVariable( + "InitialPosition", [&](size_t index_i) + { return pos_[index_i]; })), + vel_(this->particles_->vel_){}; + + virtual ~MotionConstraint(){}; + + protected: + StdLargeVec &pos_, &pos0_, &vel_; +}; +using BodyPartMotionConstraint = MotionConstraint; + +/**@class FixConstraint + * @brief Constraint with zero velocity. + */ +template +class FixConstraint : public MotionConstraint +{ + public: + explicit FixConstraint(DynamicsIdentifier &identifier) + : MotionConstraint(identifier){}; + virtual ~FixConstraint(){}; + + void update(size_t index_i, Real dt = 0.0) + { + this->pos_[index_i] = this->pos0_[index_i]; + this->vel_[index_i] = Vecd::Zero(); + }; +}; +using FixBodyConstraint = FixConstraint; +using FixBodyPartConstraint = FixConstraint; + +/** + * @class FixedInAxisDirection + * @brief Constrain the velocity of a solid body part. + */ +class FixedInAxisDirection : public MotionConstraint +{ + public: + explicit FixedInAxisDirection(BodyPartByParticle &body_part, Vecd constrained_axises = Vecd::Zero()) + : MotionConstraint(body_part), constrain_matrix_(Matd::Identity()) + { + for (int k = 0; k != Dimensions; ++k) + constrain_matrix_(k, k) = constrained_axises[k]; + }; + virtual ~FixedInAxisDirection(){}; + void update(size_t index_i, Real dt = 0.0) + { + this->vel_[index_i] = constrain_matrix_ * this->vel_[index_i]; + }; + + protected: + Matd constrain_matrix_; +}; } // namespace SPH #endif // GENERAL_CONSTRAINT_H diff --git a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.cpp b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.cpp index 88b0a97aae..04f584e9cd 100644 --- a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.cpp +++ b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.cpp @@ -7,8 +7,9 @@ namespace solid_dynamics { //=================================================================================================// SpringConstrain::SpringConstrain(BodyPartByParticle &body_part, Real stiffness) - : BaseMotionConstraint(body_part), - stiffness_(stiffness * Vecd::Ones()) {} + : MotionConstraint(body_part), + stiffness_(stiffness * Vecd::Ones()), + mass_(*particles_->getVariableByName("Mass")) {} //=================================================================================================// Vecd SpringConstrain::getAcceleration(Vecd &disp, Real mass) { @@ -28,7 +29,7 @@ void SpringConstrain::update(size_t index_i, Real dt) //=================================================================================================// PositionSolidBody:: PositionSolidBody(SPHBody &sph_body, Real start_time, Real end_time, Vecd pos_end_center) - : BaseMotionConstraint(sph_body), + : MotionConstraint(sph_body), start_time_(start_time), end_time_(end_time), pos_end_center_(pos_end_center) { BoundingBox bounds = sph_body.getSPHBodyBounds(); @@ -56,7 +57,7 @@ void PositionSolidBody::update(size_t index_i, Real dt) //=================================================================================================// PositionScaleSolidBody:: PositionScaleSolidBody(SPHBody &sph_body, Real start_time, Real end_time, Real end_scale) - : BaseMotionConstraint(sph_body), + : MotionConstraint(sph_body), start_time_(start_time), end_time_(end_time), end_scale_(end_scale) { BoundingBox bounds = sph_body.getSPHBodyBounds(); diff --git a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h index 55a69e54f9..1558f594c3 100644 --- a/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h +++ b/src/shared/particle_dynamics/solid_dynamics/constraint_dynamics.h @@ -30,11 +30,11 @@ #ifndef CONSTRAINT_DYNAMICS_H #define CONSTRAINT_DYNAMICS_H -#include "all_body_relations.h" #include "all_particle_dynamics.h" #include "all_simbody.h" #include "base_kernel.h" #include "elastic_solid.h" +#include "general_constraint.h" #include "general_reduce.h" #include "solid_body.h" #include "solid_particles.h" @@ -49,57 +49,11 @@ namespace solid_dynamics typedef DataDelegateSimple SolidDataSimple; typedef DataDelegateInner SolidDataInner; -/** - * @class BaseMotionConstraint - * @brief Base class for constraining with prescribed motion. - * Exact motion function will be defined in derive class. - * Note that, we do not impose acceleration, so that this constraint - * must be imposed after updating particle velocity by forces - * and before updating particle position. - * TODO: to clarify the treatment of particle position, - * how to achieve consistency between velocity and position constraints. - */ -template -class BaseMotionConstraint : public BaseLocalDynamics, public DataDelegationType -{ - public: - explicit BaseMotionConstraint(DynamicsIdentifier &identifier) - : BaseLocalDynamics(identifier), DataDelegationType(identifier.getSPHBody()), - pos_(this->particles_->pos_), pos0_(this->particles_->pos0_), - n_(this->particles_->n_), n0_(this->particles_->n0_), - vel_(this->particles_->vel_), force_(this->particles_->force_), mass_(this->particles_->mass_){}; - - virtual ~BaseMotionConstraint(){}; - - protected: - StdLargeVec &pos_, &pos0_; - StdLargeVec &n_, &n0_; - StdLargeVec &vel_, &force_; - StdLargeVec &mass_; -}; -using MotionConstraint = BaseMotionConstraint; - -/**@class FixConstraint - * @brief Constraint with zero velocity. - */ -template -class FixConstraint : public BaseMotionConstraint -{ - public: - explicit FixConstraint(DynamicsIdentifier &identifier) - : BaseMotionConstraint(identifier){}; - virtual ~FixConstraint(){}; - - void update(size_t index_i, Real dt = 0.0) { this->vel_[index_i] = Vecd::Zero(); }; -}; -using FixBodyConstraint = FixConstraint; -using FixBodyPartConstraint = FixConstraint; - /**@class SpringConstrain * @brief Constrain with a spring for each constrained particles to its original position. * //TODO: a test case is required for this class. */ -class SpringConstrain : public BaseMotionConstraint +class SpringConstrain : public MotionConstraint { public: SpringConstrain(BodyPartByParticle &body_part, Real stiffness); @@ -109,6 +63,7 @@ class SpringConstrain : public BaseMotionConstraint &mass_; virtual Vecd getAcceleration(Vecd &disp, Real mass); }; @@ -118,7 +73,7 @@ class SpringConstrain : public BaseMotionConstraint +class PositionSolidBody : public MotionConstraint { public: PositionSolidBody(SPHBody &sph_body, Real start_time, Real end_time, Vecd pos_end_center); @@ -139,7 +94,7 @@ class PositionSolidBody : public BaseMotionConstraint * can be considered as a quasi-static position driven boundary condition. * Note that, this constraint is not for a elastic solid body. */ -class PositionScaleSolidBody : public BaseMotionConstraint +class PositionScaleSolidBody : public MotionConstraint { public: PositionScaleSolidBody(SPHBody &sph_body, Real start_time, Real end_time, Real end_scale); @@ -161,11 +116,11 @@ class PositionScaleSolidBody : public BaseMotionConstraint -class PositionTranslate : public BaseMotionConstraint +class PositionTranslate : public MotionConstraint { public: PositionTranslate(DynamicsIdentifier &identifier, Real start_time, Real end_time, Vecd translation) - : BaseMotionConstraint(identifier), + : MotionConstraint(identifier), start_time_(start_time), end_time_(end_time), translation_(translation){}; virtual ~PositionTranslate(){}; void update(size_t index_i, Real dt = 0.0) @@ -192,37 +147,11 @@ class PositionTranslate : public BaseMotionConstraint; using TranslateSolidBodyPart = PositionTranslate; -/** - * @class FixedInAxisDirection - * @brief Constrain the velocity of a solid body part. - */ -template -class BaseFixedInAxisDirection : public BaseMotionConstraint -{ - public: - explicit BaseFixedInAxisDirection(BodyPartByParticle &body_part, Vecd constrained_axises = Vecd::Zero()) - : BaseMotionConstraint(body_part), constrain_matrix_(Matd::Identity()) - { - for (int k = 0; k != Dimensions; ++k) - constrain_matrix_(k, k) = constrained_axises[k]; - }; - virtual ~BaseFixedInAxisDirection(){}; - void update(size_t index_i, Real dt = 0.0) - { - this->vel_[index_i] = constrain_matrix_ * this->vel_[index_i]; - }; - - protected: - Matd constrain_matrix_; -}; -using FixedInAxisDirection = BaseFixedInAxisDirection; - /** * @class ConstrainSolidBodyMassCenter * @brief Constrain the mass center of a solid body. */ -template -class BaseConstrainSolidBodyMassCenter : public LocalDynamics, public DataDelegationType +class ConstrainSolidBodyMassCenter : public MotionConstraint { private: Real total_mass_; @@ -239,8 +168,8 @@ class BaseConstrainSolidBodyMassCenter : public LocalDynamics, public DataDelega } public: - explicit BaseConstrainSolidBodyMassCenter(SPHBody &sph_body, Vecd constrain_direction = Vecd::Ones()) - : LocalDynamics(sph_body), DataDelegationType(sph_body), + explicit ConstrainSolidBodyMassCenter(SPHBody &sph_body, Vecd constrain_direction = Vecd::Ones()) + : MotionConstraint(sph_body), correction_matrix_(Matd::Identity()), vel_(this->particles_->vel_), compute_total_momentum_(sph_body, "Velocity") { @@ -249,28 +178,30 @@ class BaseConstrainSolidBodyMassCenter : public LocalDynamics, public DataDelega ReduceDynamics> compute_total_mass_(sph_body, "Mass"); total_mass_ = compute_total_mass_.exec(); } - virtual ~BaseConstrainSolidBodyMassCenter(){}; + virtual ~ConstrainSolidBodyMassCenter(){}; void update(size_t index_i, Real dt = 0.0) { this->vel_[index_i] -= velocity_correction_; } }; -using ConstrainSolidBodyMassCenter = BaseConstrainSolidBodyMassCenter; + /** * @class ConstraintBySimBody * @brief Constrain by the motion computed from Simbody. */ template -class ConstraintBySimBody : public BaseMotionConstraint +class ConstraintBySimBody : public MotionConstraint { public: ConstraintBySimBody(DynamicsIdentifier &identifier, SimTK::MultibodySystem &MBsystem, SimTK::MobilizedBody &mobod, SimTK::RungeKuttaMersonIntegrator &integ) - : BaseMotionConstraint(identifier), - MBsystem_(MBsystem), mobod_(mobod), integ_(integ) + : MotionConstraint(identifier), + MBsystem_(MBsystem), mobod_(mobod), integ_(integ), + n_(*this->particles_->template getVariableByName("NormalDirection")), + n0_(*this->particles_->template getVariableByName("InitialNormalDirection")) { simbody_state_ = &integ_.getState(); MBsystem_.realize(*simbody_state_, SimTK::Stage::Acceleration); @@ -298,14 +229,15 @@ class ConstraintBySimBody : public BaseMotionConstraintpos_[index_i] = degradeToVecd(SimTKToEigen(pos)); this->vel_[index_i] = degradeToVecd(SimTKToEigen(vel)); - SimTKVec3 n = (mobod_.getBodyRotation(*simbody_state_) * EigenToSimTK(upgradeToVec3d(this->n0_[index_i]))); - this->n_[index_i] = degradeToVecd(SimTKToEigen(n)); + SimTKVec3 n = (mobod_.getBodyRotation(*simbody_state_) * EigenToSimTK(upgradeToVec3d(n0_[index_i]))); + n_[index_i] = degradeToVecd(SimTKToEigen(n)); }; protected: SimTK::MultibodySystem &MBsystem_; SimTK::MobilizedBody &mobod_; SimTK::RungeKuttaMersonIntegrator &integ_; + StdLargeVec &n_, &n0_; const SimTK::State *simbody_state_; SimTKVec3 initial_mobod_origin_location_; }; diff --git a/src/shared/particles/base_particles.h b/src/shared/particles/base_particles.h index 1b2c4d4259..945b6a083c 100644 --- a/src/shared/particles/base_particles.h +++ b/src/shared/particles/base_particles.h @@ -130,9 +130,8 @@ class BaseParticles template void registerVariable(StdLargeVec &variable_addrs, const std::string &variable_name, const InitializationFunction &initialization); - template - StdLargeVec *registerSharedVariable( - const std::string &variable_name, const DataType &default_value = ZeroData::value); + template + StdLargeVec *registerSharedVariable(const std::string &variable_name, Args &&...args); template StdLargeVec *getVariableByName(const std::string &variable_name); ParticleVariables &AllDiscreteVariables() { return all_discrete_variables_; }; diff --git a/src/shared/particles/base_particles.hpp b/src/shared/particles/base_particles.hpp index f79fb8d6d0..00bdc00a5d 100644 --- a/src/shared/particles/base_particles.hpp +++ b/src/shared/particles/base_particles.hpp @@ -79,9 +79,8 @@ DataType *BaseParticles::getSingleVariableByName(const std::string &variable_nam return nullptr; } //=================================================================================================// -template -StdLargeVec *BaseParticles:: - registerSharedVariable(const std::string &variable_name, const DataType &default_value) +template +StdLargeVec *BaseParticles::registerSharedVariable(const std::string &variable_name, Args &&...args) { DiscreteVariable *variable = findVariableByName(all_discrete_variables_, variable_name); @@ -91,7 +90,7 @@ StdLargeVec *BaseParticles:: { UniquePtrsKeeper> &container = std::get(shared_particle_data_ptrs_); StdLargeVec *contained_data = container.template createPtr>(); - registerVariable(*contained_data, variable_name, default_value); + registerVariable(*contained_data, variable_name, std::forward(args)...); return contained_data; } else diff --git a/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp b/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp index 618a91cef4..454e548d88 100644 --- a/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp +++ b/tests/2d_examples/test_2d_anisotropic_beam/test_2d_anisotropic_beam.cpp @@ -218,7 +218,7 @@ int main(int ac, char *av[]) Dynamics1Level stress_relaxation_second_half(beam_body_inner); // clamping a solid body part. BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_beam_base(beam_base); + SimpleDynamics constraint_beam_base(beam_base); //----------------------------------------------------------------------------- // outputs //----------------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_elastic_gate/elastic_gate.cpp b/tests/2d_examples/test_2d_elastic_gate/elastic_gate.cpp index 25610c0139..b3f3defcae 100644 --- a/tests/2d_examples/test_2d_elastic_gate/elastic_gate.cpp +++ b/tests/2d_examples/test_2d_elastic_gate/elastic_gate.cpp @@ -210,7 +210,7 @@ int main(int ac, char *av[]) ReduceDynamics gate_computing_time_step_size(gate); BodyRegionByParticle gate_constraint_part(gate, makeShared(createGateConstrainShape())); - SimpleDynamics gate_constraint(gate_constraint_part); + SimpleDynamics gate_constraint(gate_constraint_part); SimpleDynamics gate_update_normal(gate); //---------------------------------------------------------------------- // Define the methods for I/O operations and observations of the simulation. diff --git a/tests/2d_examples/test_2d_fsi2/fsi2.cpp b/tests/2d_examples/test_2d_fsi2/fsi2.cpp index 9bbc36536f..7b11be1b63 100644 --- a/tests/2d_examples/test_2d_fsi2/fsi2.cpp +++ b/tests/2d_examples/test_2d_fsi2/fsi2.cpp @@ -156,7 +156,7 @@ int main(int ac, char *av[]) Dynamics1Level insert_body_stress_relaxation_second_half(insert_body_inner); /** Constrain region of the inserted body. */ BodyRegionByParticle beam_base(insert_body, makeShared(createBeamBaseShape())); - SimpleDynamics constraint_beam_base(beam_base); + SimpleDynamics constraint_beam_base(beam_base); /** Update surface norm direction.*/ SimpleDynamics insert_body_update_normal(insert_body); //---------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_hydrostatic_fsi/hydrostatic_fsi.cpp b/tests/2d_examples/test_2d_hydrostatic_fsi/hydrostatic_fsi.cpp index b09863bd3f..d98858daa2 100644 --- a/tests/2d_examples/test_2d_hydrostatic_fsi/hydrostatic_fsi.cpp +++ b/tests/2d_examples/test_2d_hydrostatic_fsi/hydrostatic_fsi.cpp @@ -261,7 +261,7 @@ int main(int ac, char *av[]) Dynamics1Level gate_stress_relaxation_second_half(gate_inner); /**Constrain a solid body part. */ BodyRegionByParticle gate_constraint_part(gate, makeShared(createGateConstrainShape())); - SimpleDynamics gate_constraint(gate_constraint_part); + SimpleDynamics gate_constraint(gate_constraint_part); /** Update the norm of elastic gate. */ SimpleDynamics gate_update_normal(gate); /** Compute the average velocity of gate. */ diff --git a/tests/2d_examples/test_2d_nonlinear_wave_fsi/nonlinear_wave_fsi.h b/tests/2d_examples/test_2d_nonlinear_wave_fsi/nonlinear_wave_fsi.h index 2dcae5c502..40ccdf3c44 100644 --- a/tests/2d_examples/test_2d_nonlinear_wave_fsi/nonlinear_wave_fsi.h +++ b/tests/2d_examples/test_2d_nonlinear_wave_fsi/nonlinear_wave_fsi.h @@ -10,41 +10,41 @@ using namespace SPH; // Basic geometry parameters and numerical setup. //---------------------------------------------------------------------- Real total_physical_time = 20.0; /**< TOTAL SIMULATION TIME*/ -Real DL = 35.0; /**< Tank length. */ -Real DH = 2.0; /**< Tank height. */ -Real WL = 20.0; /**< Water block width. */ -Real WH = 0.8; /**< Water block height. */ -Real TB = 15.0; /**< Beach start. */ -Real DB = 25.0; /**< Beach end. */ -Real BEH = 2.0; /**< Beach end height. */ -Real Wmk_p = 0.0; /**< Wavemaker initial position. */ -Real EXS =2.0; /**< etra space behind the wavemaker*/ -Real StructureBasePlateH=0.12; +Real DL = 35.0; /**< Tank length. */ +Real DH = 2.0; /**< Tank height. */ +Real WL = 20.0; /**< Water block width. */ +Real WH = 0.8; /**< Water block height. */ +Real TB = 15.0; /**< Beach start. */ +Real DB = 25.0; /**< Beach end. */ +Real BEH = 2.0; /**< Beach end height. */ +Real Wmk_p = 0.0; /**< Wavemaker initial position. */ +Real EXS = 2.0; /**< etra space behind the wavemaker*/ +Real StructureBasePlateH = 0.12; /** * Initial reference particle spacing * It is a multiple of the structure baseplate height. * */ -Real particle_spacing_ref = StructureBasePlateH/4; -Real BW = particle_spacing_ref * 4.0; /**< Extending width for BCs. */ -Real Maker_width = particle_spacing_ref * 4.0; /**< Width of the wavemaker. */ +Real particle_spacing_ref = StructureBasePlateH / 4; +Real BW = particle_spacing_ref * 4.0; /**< Extending width for BCs. */ +Real Maker_width = particle_spacing_ref * 4.0; /**< Width of the wavemaker. */ -BoundingBox system_domain_bounds(Vec2d(-EXS-BW, -BW), Vec2d(DL + BW, DH + BW)); +BoundingBox system_domain_bounds(Vec2d(-EXS - BW, -BW), Vec2d(DL + BW, DH + BW)); Vec2d offset = Vec2d::Zero(); // water block parameters -Vec2d watS_lb(0.0, 0.0); /**< Left bottom. */ -Vec2d watS_lt(0.0, WH); /**< Left top. */ +Vec2d watS_lb(0.0, 0.0); /**< Left bottom. */ +Vec2d watS_lt(0.0, WH); /**< Left top. */ /** Define the corner points of the gate geometry. */ -Vec2d Wmak_lb(Wmk_p - Maker_width, 0.0); /**< Left bottom. */ -Vec2d Wmak_lt(Wmk_p - Maker_width, 1.5); /**< Left top. */ -Vec2d Wmak_rt(Wmk_p , 1.5); /**< Right top. */ -Vec2d Wmak_rb(Wmk_p , 0.0); /**< Right bottom. */ +Vec2d Wmak_lb(Wmk_p - Maker_width, 0.0); /**< Left bottom. */ +Vec2d Wmak_lt(Wmk_p - Maker_width, 1.5); /**< Left top. */ +Vec2d Wmak_rt(Wmk_p, 1.5); /**< Right top. */ +Vec2d Wmak_rb(Wmk_p, 0.0); /**< Right bottom. */ /* BEACH */ -Vec2d b1a(TB,0); /**< beach start. */ -Vec2d bwh((WH)*10+TB,WH); /**< water height at beach. */ -Vec2d b1b(25,1); /**< beach end. */ - +Vec2d b1a(TB, 0); /**< beach start. */ +Vec2d bwh((WH) * 10 + TB, WH); /**< water height at beach. */ +Vec2d b1b(25, 1); /**< beach end. */ + /* * Buoyant Structure Shape */ @@ -55,66 +55,66 @@ Real bp_y = 0.573; Real bp_l = 1.3; Real bp_h = 0.12; -Vec2d bp_lb(bp_x, bp_y); /**< Left bottom. */ -Vec2d bp_lt(bp_x, bp_y+bp_h); /**< Left top. */ -Vec2d bp_rt(bp_x+bp_l, bp_y+bp_h); /**< Right top. */ -Vec2d bp_rb(bp_x+bp_l, bp_y); /**< Right bottom. */ +Vec2d bp_lb(bp_x, bp_y); /**< Left bottom. */ +Vec2d bp_lt(bp_x, bp_y + bp_h); /**< Left top. */ +Vec2d bp_rt(bp_x + bp_l, bp_y + bp_h); /**< Right top. */ +Vec2d bp_rb(bp_x + bp_l, bp_y); /**< Right bottom. */ /** Seaside pillar . */ -Real ssp_x = bp_x+0.25; -Real ssp_y = bp_y+bp_h; +Real ssp_x = bp_x + 0.25; +Real ssp_y = bp_y + bp_h; Real ssp_l = 0.2; Real ssp_h = 0.24; -Vec2d ssp_lb(ssp_x, ssp_y); /**< Left bottom. */ -Vec2d ssp_lt(ssp_x, ssp_y+ssp_h); /**< Left top. */ -Vec2d ssp_rt(ssp_x+ssp_l, ssp_y+ssp_h); /**< Right top. */ -Vec2d ssp_rb(ssp_x+ssp_l, ssp_y); /**< Right bottom. */ +Vec2d ssp_lb(ssp_x, ssp_y); /**< Left bottom. */ +Vec2d ssp_lt(ssp_x, ssp_y + ssp_h); /**< Left top. */ +Vec2d ssp_rt(ssp_x + ssp_l, ssp_y + ssp_h); /**< Right top. */ +Vec2d ssp_rb(ssp_x + ssp_l, ssp_y); /**< Right bottom. */ /** Portside pillar . */ -Real psp_x = bp_x+bp_l-0.45; -Real psp_y = bp_y+bp_h; +Real psp_x = bp_x + bp_l - 0.45; +Real psp_y = bp_y + bp_h; Real psp_l = 0.2; Real psp_h = 0.24; -Vec2d psp_lb(psp_x, psp_y); /**< Left bottom. */ -Vec2d psp_lt(psp_x, psp_y+psp_h); /**< Left top. */ -Vec2d psp_rt(psp_x+psp_l, psp_y+psp_h); /**< Right top. */ -Vec2d psp_rb(psp_x+psp_l, psp_y); /**< Right bottom. */ +Vec2d psp_lb(psp_x, psp_y); /**< Left bottom. */ +Vec2d psp_lt(psp_x, psp_y + psp_h); /**< Left top. */ +Vec2d psp_rt(psp_x + psp_l, psp_y + psp_h); /**< Right top. */ +Vec2d psp_rb(psp_x + psp_l, psp_y); /**< Right bottom. */ /** Topplate. */ -Real tp_x = bp_x+0.18; -Real tp_y = bp_y+0.36; +Real tp_x = bp_x + 0.18; +Real tp_y = bp_y + 0.36; Real tp_l = 0.94; Real tp_h = 0.11; -Vec2d tp_lb(tp_x, tp_y); /**< Left bottom. */ -Vec2d tp_lt(tp_x, tp_y+tp_h); /**< Left top. */ -Vec2d tp_rt(tp_x+tp_l, tp_y+tp_h); /**< Right top. */ -Vec2d tp_rb(tp_x+tp_l, tp_y); /**< Right bottom. */ +Vec2d tp_lb(tp_x, tp_y); /**< Left bottom. */ +Vec2d tp_lt(tp_x, tp_y + tp_h); /**< Left top. */ +Vec2d tp_rt(tp_x + tp_l, tp_y + tp_h); /**< Right top. */ +Vec2d tp_rb(tp_x + tp_l, tp_y); /**< Right bottom. */ /** * Topology of the tethers. * */ -Real cxA = bp_x+0.35; /**< Center of buoy in x direction. */ -Real cyA = bp_y; /**< Center of buoy in y direction. */ -Vec2d tethering_pointA(cxA, 0.0); /**< The ground tethering point. */ -Vec2d cable_endA(cxA,cyA); /**< The structure tethering point. */ +Real cxA = bp_x + 0.35; /**< Center of buoy in x direction. */ +Real cyA = bp_y; /**< Center of buoy in y direction. */ +Vec2d tethering_pointA(cxA, 0.0); /**< The ground tethering point. */ +Vec2d cable_endA(cxA, cyA); /**< The structure tethering point. */ -Real cxB = bp_x+bp_l-0.35; /**< Center of buoy in x direction. */ -Real cyB = bp_y; /**< Center of buoy in y direction. */ -Vec2d tethering_pointB(cxB, 0.0); /**< The ground tethering point. */ -Vec2d cable_endB(cxB,cyB); /**< The structure tethering point. */ +Real cxB = bp_x + bp_l - 0.35; /**< Center of buoy in x direction. */ +Real cyB = bp_y; /**< Center of buoy in y direction. */ +Vec2d tethering_pointB(cxB, 0.0); /**< The ground tethering point. */ +Vec2d cable_endB(cxB, cyB); /**< The structure tethering point. */ -Real cablength=bp_y; +Real cablength = bp_y; //---------------------------------------------------------------------- // Material properties of the fluid. //---------------------------------------------------------------------- -Real rho0_f = 1000.0; /**< Reference density of fluid. */ -Real gravity_g = 9.81; /**< Value of gravity. */ -Real U_f = 2.0 * sqrt(0.79*gravity_g); /**< Characteristic velocity. */ -Real c_f = 10.0 * U_f; /**< Reference sound speed. */ +Real rho0_f = 1000.0; /**< Reference density of fluid. */ +Real gravity_g = 9.81; /**< Value of gravity. */ +Real U_f = 2.0 * sqrt(0.79 * gravity_g); /**< Characteristic velocity. */ +Real c_f = 10.0 * U_f; /**< Reference sound speed. */ Real mu_f = 1.0e-3; //---------------------------------------------------------------------- @@ -123,79 +123,79 @@ Real mu_f = 1.0e-3; /* Weight of the solid structure*/ Real StructureMass = 62.036; -Real bp_area =bp_l*bp_h; -Real ssp_area =ssp_l*ssp_h; -Real psp_area =psp_l*psp_h; -Real tp_area =tp_l*tp_h; +Real bp_area = bp_l * bp_h; +Real ssp_area = ssp_l * ssp_h; +Real psp_area = psp_l * psp_h; +Real tp_area = tp_l * tp_h; /* Surface of the solid structure*/ -Real Area=bp_area+ssp_area+psp_area+tp_area; +Real Area = bp_area + ssp_area + psp_area + tp_area; /**< Density of the solid structure*/ -Real rho_s=StructureMass/Area; +Real rho_s = StructureMass / Area; -Real rho_bp=rho_s; -Vec2d bp_cm(bp_x+bp_l/2,bp_y+bp_h/2); +Real rho_bp = rho_s; +Vec2d bp_cm(bp_x + bp_l / 2, bp_y + bp_h / 2); Vec3d Ibp( - (rho_bp*bp_area)/12*(bp_l*bp_h*bp_h*bp_h), - (rho_bp*bp_area)/12*(bp_h*bp_l*bp_l*bp_l), - (rho_bp*bp_area)/12*(bp_l*bp_l+bp_h*bp_h)); + (rho_bp * bp_area) / 12 * (bp_l * bp_h * bp_h * bp_h), + (rho_bp * bp_area) / 12 * (bp_h * bp_l * bp_l * bp_l), + (rho_bp * bp_area) / 12 * (bp_l * bp_l + bp_h * bp_h)); -Real rho_ssp=rho_s; -Vec2d ssp_cm(ssp_x+ssp_l/2,ssp_y+ssp_h/2); +Real rho_ssp = rho_s; +Vec2d ssp_cm(ssp_x + ssp_l / 2, ssp_y + ssp_h / 2); Vec3d Issp( - (rho_s*ssp_area)/12*(ssp_l*ssp_h*ssp_h*ssp_h), - (rho_s*ssp_area)/12*(ssp_h*ssp_l*ssp_l*ssp_l), - (rho_s*ssp_area)/12*(ssp_l*ssp_l+ssp_h*ssp_h)); + (rho_s * ssp_area) / 12 * (ssp_l * ssp_h * ssp_h * ssp_h), + (rho_s * ssp_area) / 12 * (ssp_h * ssp_l * ssp_l * ssp_l), + (rho_s * ssp_area) / 12 * (ssp_l * ssp_l + ssp_h * ssp_h)); -Real rho_psp=rho_s; -Vec2d psp_cm(psp_x+psp_l/2,psp_y+psp_h/2); +Real rho_psp = rho_s; +Vec2d psp_cm(psp_x + psp_l / 2, psp_y + psp_h / 2); Vec3d Ipsp( - (rho_s*psp_area)/12*(psp_l*psp_h*psp_h*psp_h), - (rho_s*psp_area)/12*(psp_h*psp_l*psp_l*psp_l), - (rho_s*psp_area)/12*(psp_l*psp_l+psp_h*psp_h)); + (rho_s * psp_area) / 12 * (psp_l * psp_h * psp_h * psp_h), + (rho_s * psp_area) / 12 * (psp_h * psp_l * psp_l * psp_l), + (rho_s * psp_area) / 12 * (psp_l * psp_l + psp_h * psp_h)); -Real rho_tp=rho_s; -Vec2d tp_cm(tp_x+tp_l/2,tp_y+tp_h/2); +Real rho_tp = rho_s; +Vec2d tp_cm(tp_x + tp_l / 2, tp_y + tp_h / 2); Vec3d Itp( - (rho_s*tp_area)/12*(tp_l*tp_h*tp_h*tp_h), - (rho_s*tp_area)/12*(tp_h*tp_l*tp_l*tp_l), - (rho_s*tp_area)/12*(tp_l*tp_l+tp_h*tp_h)); - -Real bcmx=(rho_bp*bp_area*bp_cm[0]+ - rho_ssp*ssp_area*ssp_cm[0]+ - rho_psp*psp_area*psp_cm[0]+ - rho_tp*tp_area*tp_cm[0])/ - (rho_bp*bp_area+ - rho_ssp*ssp_area+ - rho_psp*psp_area+ - rho_tp*tp_area); -Real bcmy=(rho_bp*bp_area*bp_cm[1]+ - rho_ssp*ssp_area*ssp_cm[1]+ - rho_psp*psp_area*psp_cm[1]+ - rho_tp*tp_area*tp_cm[1])/ - (rho_bp*bp_area+ - rho_ssp*ssp_area+ - rho_psp*psp_area+ - rho_tp*tp_area); - -Vec2d G(bcmx,bcmy); -Real d_bp =sqrt((G[0]-bp_cm[0])*(G[0]-bp_cm[0])+(G[1]-bp_cm[1])*(G[1]-bp_cm[1])); -Real d_ssp=sqrt((G[0]-ssp_cm[0])*(G[0]-ssp_cm[0])+(G[1]-ssp_cm[1])*(G[1]-ssp_cm[1])); -Real d_psp=sqrt((G[0]-psp_cm[0])*(G[0]-psp_cm[0])+(G[1]-psp_cm[1])*(G[1]-psp_cm[1])); -Real d_tp =sqrt((G[0]-tp_cm[0])*(G[0]-tp_cm[0])+(G[1]-tp_cm[1])*(G[1]-tp_cm[1])); - -Real Ix = (Ibp[0]+rho_bp*bp_area*(d_bp*d_bp)+ - Issp[0]+rho_ssp*ssp_area*(d_ssp*d_ssp)+ - Ipsp[0]+rho_psp*psp_area*(d_psp*d_psp)+ - Itp[0]+rho_tp*tp_area*(d_tp*d_tp)); -Real Iy = (Ibp[1]+rho_bp*bp_area*(d_bp*d_bp)+ - Issp[1]+rho_ssp*ssp_area*(d_ssp*d_ssp)+ - Ipsp[1]+rho_psp*psp_area*(d_psp*d_psp)+ - Itp[1]+rho_tp*tp_area*(d_tp*d_tp)); -Real Iz = (Ibp[2]+rho_bp*bp_area*(d_bp*d_bp)+ - Issp[2]+rho_ssp*ssp_area*(d_ssp*d_ssp)+ - Ipsp[2]+rho_psp*psp_area*(d_psp*d_psp)+ - Itp[2]+rho_tp*tp_area*(d_tp*d_tp)); + (rho_s * tp_area) / 12 * (tp_l * tp_h * tp_h * tp_h), + (rho_s * tp_area) / 12 * (tp_h * tp_l * tp_l * tp_l), + (rho_s * tp_area) / 12 * (tp_l * tp_l + tp_h * tp_h)); + +Real bcmx = (rho_bp * bp_area * bp_cm[0] + + rho_ssp * ssp_area * ssp_cm[0] + + rho_psp * psp_area * psp_cm[0] + + rho_tp * tp_area * tp_cm[0]) / + (rho_bp * bp_area + + rho_ssp * ssp_area + + rho_psp * psp_area + + rho_tp * tp_area); +Real bcmy = (rho_bp * bp_area * bp_cm[1] + + rho_ssp * ssp_area * ssp_cm[1] + + rho_psp * psp_area * psp_cm[1] + + rho_tp * tp_area * tp_cm[1]) / + (rho_bp * bp_area + + rho_ssp * ssp_area + + rho_psp * psp_area + + rho_tp * tp_area); + +Vec2d G(bcmx, bcmy); +Real d_bp = sqrt((G[0] - bp_cm[0]) * (G[0] - bp_cm[0]) + (G[1] - bp_cm[1]) * (G[1] - bp_cm[1])); +Real d_ssp = sqrt((G[0] - ssp_cm[0]) * (G[0] - ssp_cm[0]) + (G[1] - ssp_cm[1]) * (G[1] - ssp_cm[1])); +Real d_psp = sqrt((G[0] - psp_cm[0]) * (G[0] - psp_cm[0]) + (G[1] - psp_cm[1]) * (G[1] - psp_cm[1])); +Real d_tp = sqrt((G[0] - tp_cm[0]) * (G[0] - tp_cm[0]) + (G[1] - tp_cm[1]) * (G[1] - tp_cm[1])); + +Real Ix = (Ibp[0] + rho_bp * bp_area * (d_bp * d_bp) + + Issp[0] + rho_ssp * ssp_area * (d_ssp * d_ssp) + + Ipsp[0] + rho_psp * psp_area * (d_psp * d_psp) + + Itp[0] + rho_tp * tp_area * (d_tp * d_tp)); +Real Iy = (Ibp[1] + rho_bp * bp_area * (d_bp * d_bp) + + Issp[1] + rho_ssp * ssp_area * (d_ssp * d_ssp) + + Ipsp[1] + rho_psp * psp_area * (d_psp * d_psp) + + Itp[1] + rho_tp * tp_area * (d_tp * d_tp)); +Real Iz = (Ibp[2] + rho_bp * bp_area * (d_bp * d_bp) + + Issp[2] + rho_ssp * ssp_area * (d_ssp * d_ssp) + + Ipsp[2] + rho_psp * psp_area * (d_psp * d_psp) + + Itp[2] + rho_tp * tp_area * (d_tp * d_tp)); /** * Structure observer position @@ -207,98 +207,97 @@ Vec2d obs = bp_cm; // geometric shape elements used in the case //------------------------------------------------------------------------------ MultiPolygon createStructureShape() -{ - /** Geometry definition. */ - std::vector sructure_bp; - sructure_bp.push_back(bp_lb); - sructure_bp.push_back(bp_lt); - sructure_bp.push_back(bp_rt); - sructure_bp.push_back(bp_rb); - sructure_bp.push_back(bp_lb); - - std::vector sructure_ssp; - sructure_ssp.push_back(ssp_lb); - sructure_ssp.push_back(ssp_lt); - sructure_ssp.push_back(ssp_rt); - sructure_ssp.push_back(ssp_rb); - sructure_ssp.push_back(ssp_lb); - - std::vector sructure_psp; - sructure_psp.push_back(psp_lb); - sructure_psp.push_back(psp_lt); - sructure_psp.push_back(psp_rt); - sructure_psp.push_back(psp_rb); - sructure_psp.push_back(psp_lb); - - std::vector sructure_tp; - sructure_tp.push_back(tp_lb); - sructure_tp.push_back(tp_lt); - sructure_tp.push_back(tp_rt); - sructure_tp.push_back(tp_rb); - sructure_tp.push_back(tp_lb); - - MultiPolygon multi_polygon_; - - multi_polygon_.addAPolygon(sructure_bp, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_ssp, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_psp, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_tp, ShapeBooleanOps::add); - - return multi_polygon_; +{ + /** Geometry definition. */ + std::vector sructure_bp; + sructure_bp.push_back(bp_lb); + sructure_bp.push_back(bp_lt); + sructure_bp.push_back(bp_rt); + sructure_bp.push_back(bp_rb); + sructure_bp.push_back(bp_lb); + + std::vector sructure_ssp; + sructure_ssp.push_back(ssp_lb); + sructure_ssp.push_back(ssp_lt); + sructure_ssp.push_back(ssp_rt); + sructure_ssp.push_back(ssp_rb); + sructure_ssp.push_back(ssp_lb); + + std::vector sructure_psp; + sructure_psp.push_back(psp_lb); + sructure_psp.push_back(psp_lt); + sructure_psp.push_back(psp_rt); + sructure_psp.push_back(psp_rb); + sructure_psp.push_back(psp_lb); + + std::vector sructure_tp; + sructure_tp.push_back(tp_lb); + sructure_tp.push_back(tp_lt); + sructure_tp.push_back(tp_rt); + sructure_tp.push_back(tp_rb); + sructure_tp.push_back(tp_lb); + + MultiPolygon multi_polygon_; + + multi_polygon_.addAPolygon(sructure_bp, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_ssp, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_psp, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_tp, ShapeBooleanOps::add); + + return multi_polygon_; } class FloatingStructure : public MultiPolygonShape { -public: - explicit FloatingStructure(const std::string &shape_name) : MultiPolygonShape(shape_name) - { - /** Geometry definition. */ - std::vector sructure_bp; - sructure_bp.push_back(bp_lb); - sructure_bp.push_back(bp_lt); - sructure_bp.push_back(bp_rt); - sructure_bp.push_back(bp_rb); - sructure_bp.push_back(bp_lb); - - std::vector sructure_ssp; - sructure_ssp.push_back(ssp_lb); - sructure_ssp.push_back(ssp_lt); - sructure_ssp.push_back(ssp_rt); - sructure_ssp.push_back(ssp_rb); - sructure_ssp.push_back(ssp_lb); - - std::vector sructure_psp; - sructure_psp.push_back(psp_lb); - sructure_psp.push_back(psp_lt); - sructure_psp.push_back(psp_rt); - sructure_psp.push_back(psp_rb); - sructure_psp.push_back(psp_lb); - - std::vector sructure_tp; - sructure_tp.push_back(tp_lb); - sructure_tp.push_back(tp_lt); - sructure_tp.push_back(tp_rt); - sructure_tp.push_back(tp_rb); - sructure_tp.push_back(tp_lb); - - multi_polygon_.addAPolygon(sructure_bp, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_ssp, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_psp, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_tp, ShapeBooleanOps::add); - - } + public: + explicit FloatingStructure(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + /** Geometry definition. */ + std::vector sructure_bp; + sructure_bp.push_back(bp_lb); + sructure_bp.push_back(bp_lt); + sructure_bp.push_back(bp_rt); + sructure_bp.push_back(bp_rb); + sructure_bp.push_back(bp_lb); + + std::vector sructure_ssp; + sructure_ssp.push_back(ssp_lb); + sructure_ssp.push_back(ssp_lt); + sructure_ssp.push_back(ssp_rt); + sructure_ssp.push_back(ssp_rb); + sructure_ssp.push_back(ssp_lb); + + std::vector sructure_psp; + sructure_psp.push_back(psp_lb); + sructure_psp.push_back(psp_lt); + sructure_psp.push_back(psp_rt); + sructure_psp.push_back(psp_rb); + sructure_psp.push_back(psp_lb); + + std::vector sructure_tp; + sructure_tp.push_back(tp_lb); + sructure_tp.push_back(tp_lt); + sructure_tp.push_back(tp_rt); + sructure_tp.push_back(tp_rb); + sructure_tp.push_back(tp_lb); + + multi_polygon_.addAPolygon(sructure_bp, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_ssp, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_psp, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_tp, ShapeBooleanOps::add); + } }; class StructureSystemForSimbody : public SolidBodyPartForSimbody { -public: - StructureSystemForSimbody(SPHBody &sph_body, SharedPtr shape_ptr) - : SolidBodyPartForSimbody(sph_body, shape_ptr) - { - body_part_mass_properties_ = - mass_properties_ptr_keeper_ - .createPtr(StructureMass, SimTK::Vec3(0.0), SimTK::UnitInertia(0, 0, Iz)); - } + public: + StructureSystemForSimbody(SPHBody &sph_body, SharedPtr shape_ptr) + : SolidBodyPartForSimbody(sph_body, shape_ptr) + { + body_part_mass_properties_ = + mass_properties_ptr_keeper_ + .createPtr(StructureMass, SimTK::Vec3(0.0), SimTK::UnitInertia(0, 0, Iz)); + } }; //---------------------------------------------------------------------- @@ -306,62 +305,61 @@ class StructureSystemForSimbody : public SolidBodyPartForSimbody //---------------------------------------------------------------------- class WaterBlock : public MultiPolygonShape { -public: - explicit WaterBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) - { - /** Geometry definition. */ - std::vector water_block_shape; - water_block_shape.push_back(watS_lb); - water_block_shape.push_back(watS_lt); - water_block_shape.push_back(bwh); - water_block_shape.push_back(b1a); - water_block_shape.push_back(watS_lb); - - /*Structure substract*/ - - std::vector sructure_bp; - sructure_bp.push_back(bp_lb); - sructure_bp.push_back(bp_lt); - sructure_bp.push_back(bp_rt); - sructure_bp.push_back(bp_rb); - sructure_bp.push_back(bp_lb); - - std::vector sructure_ssp; - sructure_ssp.push_back(ssp_lb); - sructure_ssp.push_back(ssp_lt); - sructure_ssp.push_back(ssp_rt); - sructure_ssp.push_back(ssp_rb); - sructure_ssp.push_back(ssp_lb); - - std::vector sructure_psp; - sructure_psp.push_back(psp_lb); - sructure_psp.push_back(psp_lt); - sructure_psp.push_back(psp_rt); - sructure_psp.push_back(psp_rb); - sructure_psp.push_back(psp_lb); - - std::vector sructure_tp; - sructure_tp.push_back(tp_lb); - sructure_tp.push_back(tp_lt); - sructure_tp.push_back(tp_rt); - sructure_tp.push_back(tp_rb); - sructure_tp.push_back(tp_lb); - - std::vector sructure_mdp; - sructure_mdp.push_back(ssp_rb); - sructure_mdp.push_back(ssp_rt); - sructure_mdp.push_back(psp_lt); - sructure_mdp.push_back(psp_lb); - sructure_mdp.push_back(ssp_rb); - - multi_polygon_.addAPolygon(water_block_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(sructure_bp, ShapeBooleanOps::sub); - multi_polygon_.addAPolygon(sructure_ssp, ShapeBooleanOps::sub); - multi_polygon_.addAPolygon(sructure_psp, ShapeBooleanOps::sub); - multi_polygon_.addAPolygon(sructure_tp, ShapeBooleanOps::sub); - multi_polygon_.addAPolygon(sructure_mdp, ShapeBooleanOps::sub); - - } + public: + explicit WaterBlock(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + /** Geometry definition. */ + std::vector water_block_shape; + water_block_shape.push_back(watS_lb); + water_block_shape.push_back(watS_lt); + water_block_shape.push_back(bwh); + water_block_shape.push_back(b1a); + water_block_shape.push_back(watS_lb); + + /*Structure substract*/ + + std::vector sructure_bp; + sructure_bp.push_back(bp_lb); + sructure_bp.push_back(bp_lt); + sructure_bp.push_back(bp_rt); + sructure_bp.push_back(bp_rb); + sructure_bp.push_back(bp_lb); + + std::vector sructure_ssp; + sructure_ssp.push_back(ssp_lb); + sructure_ssp.push_back(ssp_lt); + sructure_ssp.push_back(ssp_rt); + sructure_ssp.push_back(ssp_rb); + sructure_ssp.push_back(ssp_lb); + + std::vector sructure_psp; + sructure_psp.push_back(psp_lb); + sructure_psp.push_back(psp_lt); + sructure_psp.push_back(psp_rt); + sructure_psp.push_back(psp_rb); + sructure_psp.push_back(psp_lb); + + std::vector sructure_tp; + sructure_tp.push_back(tp_lb); + sructure_tp.push_back(tp_lt); + sructure_tp.push_back(tp_rt); + sructure_tp.push_back(tp_rb); + sructure_tp.push_back(tp_lb); + + std::vector sructure_mdp; + sructure_mdp.push_back(ssp_rb); + sructure_mdp.push_back(ssp_rt); + sructure_mdp.push_back(psp_lt); + sructure_mdp.push_back(psp_lb); + sructure_mdp.push_back(ssp_rb); + + multi_polygon_.addAPolygon(water_block_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(sructure_bp, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(sructure_ssp, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(sructure_psp, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(sructure_tp, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(sructure_mdp, ShapeBooleanOps::sub); + } }; //---------------------------------------------------------------------- @@ -369,40 +367,39 @@ class WaterBlock : public MultiPolygonShape //---------------------------------------------------------------------- class WallBoundary : public MultiPolygonShape { -public: - explicit WallBoundary(const std::string &shape_name) : MultiPolygonShape(shape_name) - { - /** Geometry definition. */ - std::vector outer_wall_shape; - outer_wall_shape.push_back(Vecd(-EXS, 0.0)+Vec2d(-BW,-BW)); - outer_wall_shape.push_back(Vec2d(-EXS,DH)+Vec2d(-BW,0.0)); - outer_wall_shape.push_back(Vec2d(DL,DH)+Vec2d(+BW,0.0)); - outer_wall_shape.push_back(Vec2d(DL,BEH)+Vec2d(+BW,-BW)); - outer_wall_shape.push_back(b1b+Vec2d(0.0,-BW)); - outer_wall_shape.push_back(b1a+Vec2d(-BW,-BW)); - outer_wall_shape.push_back(Vec2d(-EXS,0.0)+Vec2d(-BW,-BW)); - - std::vector inner_wall_shape; - inner_wall_shape.push_back(Vecd(-EXS, 0.0)); - inner_wall_shape.push_back(Vec2d(-EXS,DH)); - inner_wall_shape.push_back(Vec2d(DL,DH)); - inner_wall_shape.push_back(Vec2d(DL,BEH)); - inner_wall_shape.push_back(b1b); - inner_wall_shape.push_back(b1a); - inner_wall_shape.push_back(Vec2d(-EXS,0.0)); - - std::vector Wmak_shape; - Wmak_shape.push_back(Wmak_lb); - Wmak_shape.push_back(Wmak_lt); - Wmak_shape.push_back(Wmak_rt); - Wmak_shape.push_back(Wmak_rb); - Wmak_shape.push_back(Wmak_lb); - - multi_polygon_.addAPolygon(outer_wall_shape, ShapeBooleanOps::add); - multi_polygon_.addAPolygon(inner_wall_shape, ShapeBooleanOps::sub); - multi_polygon_.addAPolygon(Wmak_shape, ShapeBooleanOps::add); - - } + public: + explicit WallBoundary(const std::string &shape_name) : MultiPolygonShape(shape_name) + { + /** Geometry definition. */ + std::vector outer_wall_shape; + outer_wall_shape.push_back(Vecd(-EXS, 0.0) + Vec2d(-BW, -BW)); + outer_wall_shape.push_back(Vec2d(-EXS, DH) + Vec2d(-BW, 0.0)); + outer_wall_shape.push_back(Vec2d(DL, DH) + Vec2d(+BW, 0.0)); + outer_wall_shape.push_back(Vec2d(DL, BEH) + Vec2d(+BW, -BW)); + outer_wall_shape.push_back(b1b + Vec2d(0.0, -BW)); + outer_wall_shape.push_back(b1a + Vec2d(-BW, -BW)); + outer_wall_shape.push_back(Vec2d(-EXS, 0.0) + Vec2d(-BW, -BW)); + + std::vector inner_wall_shape; + inner_wall_shape.push_back(Vecd(-EXS, 0.0)); + inner_wall_shape.push_back(Vec2d(-EXS, DH)); + inner_wall_shape.push_back(Vec2d(DL, DH)); + inner_wall_shape.push_back(Vec2d(DL, BEH)); + inner_wall_shape.push_back(b1b); + inner_wall_shape.push_back(b1a); + inner_wall_shape.push_back(Vec2d(-EXS, 0.0)); + + std::vector Wmak_shape; + Wmak_shape.push_back(Wmak_lb); + Wmak_shape.push_back(Wmak_lt); + Wmak_shape.push_back(Wmak_rt); + Wmak_shape.push_back(Wmak_rb); + Wmak_shape.push_back(Wmak_lb); + + multi_polygon_.addAPolygon(outer_wall_shape, ShapeBooleanOps::add); + multi_polygon_.addAPolygon(inner_wall_shape, ShapeBooleanOps::sub); + multi_polygon_.addAPolygon(Wmak_shape, ShapeBooleanOps::add); + } }; //---------------------------------------------------------------------- @@ -410,159 +407,174 @@ class WallBoundary : public MultiPolygonShape //---------------------------------------------------------------------- MultiPolygon createWaveMakerShape() { - std::vector Wmak_shape; - Wmak_shape.push_back(Wmak_lb); - Wmak_shape.push_back(Wmak_lt); - Wmak_shape.push_back(Wmak_rt); - Wmak_shape.push_back(Wmak_rb); - Wmak_shape.push_back(Wmak_lb); - - MultiPolygon multi_polygon; - multi_polygon.addAPolygon(Wmak_shape, ShapeBooleanOps::add); - return multi_polygon; + std::vector Wmak_shape; + Wmak_shape.push_back(Wmak_lb); + Wmak_shape.push_back(Wmak_lt); + Wmak_shape.push_back(Wmak_rt); + Wmak_shape.push_back(Wmak_rb); + Wmak_shape.push_back(Wmak_lb); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(Wmak_shape, ShapeBooleanOps::add); + return multi_polygon; } //---------------------------------------------------------------------- // Boundary condition for wavemaker //---------------------------------------------------------------------- -class WaveMaking : public solid_dynamics::MotionConstraint -{ - Real h; - Real tf; - Real xf; - Real fmn; - Real fmx; - Real a; - int N; - std::vector om; - std::vector S; - std::vector k; - Real g; - - Vecd getDisplacement(const Real &time) - { - Real dp=0; - for (int j=0;j<(N);j++){ - dp = dp + 0.5*S[j]*cos(-k[j]*xf-om[j]*(time-tf)); - }; - - Vecd displacement{Vecd::Zero()}; - displacement[0] = dp; - return displacement; - } - - Vec2d getVelocity(const Real &time) - { - Real vl=0; - for (int j=0;j<(N);j++){ - vl = vl + 0.5*om[j]*S[j]*sin(-k[j]*xf-om[j]*(time-tf)); - }; - Vec2d velocity{Vecd::Zero()}; - velocity[0] = vl; - return velocity; - } - - Vec2d getAcceleration(const Real &time) - { - Real ax=0; - for (int j=0;j<(N);j++){ - ax = ax -0.5*om[j]*om[j]*S[j]*cos(-k[j]*xf-om[j]*(time-tf)); - }; - Vec2d acceleration{Vecd::Zero()}; - acceleration[0]=ax; - return acceleration; - } - - Real OBJ(Real wnmb,Real omsq){ - return omsq-g*wnmb*tanh(wnmb*h); - }; - - void ComputeWaveChar(){ - std::vector f(N); - - - for (int i=0;i<(N);i++){ - f[i]=fmn+i*(fmx-fmn)/N; - }; - - om=f; - for (int i=0;i<(N);i++){ - om[i]=2*PI*f[i]; - }; - - k=om; - std::vector OBJ_(N); - OBJ_=om; - for (int i=0;i<(N);i++){ - //Solve dispersion equation - Real omsq=om[i]*om[i]; - Real kmin=0; - Real kmax=20; - - Real tol=1E-15; - Real wnmb=kmin; - while ((kmax-kmin)>=tol){ - - wnmb=(kmin+kmax)/2; - OBJ_[i] = OBJ(wnmb,omsq); - - Real OBJmin = OBJ(kmin,omsq); - Real OBJmax = OBJ(wnmb,omsq); - - if (abs(OBJ_[i])<=tol){ - break; - } - else if(OBJmin*OBJmax<0){ - kmax=wnmb; - } - else{ - kmin=wnmb; - }; - }; - - k[i]=wnmb; - - }; - - S=f; - for (int i=0;i<(N);i++){ - S[i]=a*(sinh(k[i]*h)*cosh(k[i]*h)+k[i]*h)/(sinh(k[i]*h)*sinh(k[i]*h)); - }; - }; - - - -public: - WaveMaking(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part), - h(WH),tf(20.480),xf(12.0),fmn(0.32),fmx(0.96),a(0.0068),N(32),g(gravity_g) - { - ComputeWaveChar(); - } - - void update(size_t index_i, Real dt = 0.0) - { - Real time = GlobalStaticVariables::physical_time_; - pos_[index_i] = pos0_[index_i] + getDisplacement(time); - vel_[index_i] = getVelocity(time); - force_[index_i] = mass_[index_i] * getAcceleration(time); - }; +class WaveMaking : public BodyPartMotionConstraint +{ + Real h; + Real tf; + Real xf; + Real fmn; + Real fmx; + Real a; + int N; + std::vector om; + std::vector S; + std::vector k; + Real g; + + Vecd getDisplacement(const Real &time) + { + Real dp = 0; + for (int j = 0; j < (N); j++) + { + dp = dp + 0.5 * S[j] * cos(-k[j] * xf - om[j] * (time - tf)); + }; + + Vecd displacement{Vecd::Zero()}; + displacement[0] = dp; + return displacement; + } + + Vec2d getVelocity(const Real &time) + { + Real vl = 0; + for (int j = 0; j < (N); j++) + { + vl = vl + 0.5 * om[j] * S[j] * sin(-k[j] * xf - om[j] * (time - tf)); + }; + Vec2d velocity{Vecd::Zero()}; + velocity[0] = vl; + return velocity; + } + + Vec2d getAcceleration(const Real &time) + { + Real ax = 0; + for (int j = 0; j < (N); j++) + { + ax = ax - 0.5 * om[j] * om[j] * S[j] * cos(-k[j] * xf - om[j] * (time - tf)); + }; + Vec2d acceleration{Vecd::Zero()}; + acceleration[0] = ax; + return acceleration; + } + + Real OBJ(Real wnmb, Real omsq) + { + return omsq - g * wnmb * tanh(wnmb * h); + }; + + void ComputeWaveChar() + { + std::vector f(N); + + for (int i = 0; i < (N); i++) + { + f[i] = fmn + i * (fmx - fmn) / N; + }; + + om = f; + for (int i = 0; i < (N); i++) + { + om[i] = 2 * PI * f[i]; + }; + + k = om; + std::vector OBJ_(N); + OBJ_ = om; + for (int i = 0; i < (N); i++) + { + // Solve dispersion equation + Real omsq = om[i] * om[i]; + Real kmin = 0; + Real kmax = 20; + + Real tol = 1E-15; + Real wnmb = kmin; + while ((kmax - kmin) >= tol) + { + + wnmb = (kmin + kmax) / 2; + OBJ_[i] = OBJ(wnmb, omsq); + + Real OBJmin = OBJ(kmin, omsq); + Real OBJmax = OBJ(wnmb, omsq); + + if (abs(OBJ_[i]) <= tol) + { + break; + } + else if (OBJmin * OBJmax < 0) + { + kmax = wnmb; + } + else + { + kmin = wnmb; + }; + }; + + k[i] = wnmb; + }; + + S = f; + for (int i = 0; i < (N); i++) + { + S[i] = a * (sinh(k[i] * h) * cosh(k[i] * h) + k[i] * h) / (sinh(k[i] * h) * sinh(k[i] * h)); + }; + }; + + public: + WaveMaking(BodyPartByParticle &body_part) + : BodyPartMotionConstraint(body_part), + h(WH), tf(20.480), xf(12.0), fmn(0.32), fmx(0.96), a(0.0068), N(32), g(gravity_g), + mass_(*particles_->getVariableByName("Mass")), + force_(*particles_->getVariableByName("Force")) + { + ComputeWaveChar(); + } + + void update(size_t index_i, Real dt = 0.0) + { + Real time = GlobalStaticVariables::physical_time_; + pos_[index_i] = pos0_[index_i] + getDisplacement(time); + vel_[index_i] = getVelocity(time); + force_[index_i] = mass_[index_i] * getAcceleration(time); + }; + + protected: + StdLargeVec &mass_; + StdLargeVec &force_; }; //---------------------------------------------------------------------- -// create mesuring probes +// create measuring probes //---------------------------------------------------------------------- Real h = 1.3 * particle_spacing_ref; MultiPolygon createWaveGauge() { - std::vector pnts; - pnts.push_back(Vecd(10.848 - h, 0.0)); - pnts.push_back(Vecd(10.848 - h, 1.4)); - pnts.push_back(Vecd(10.848 + h, 1.4)); - pnts.push_back(Vecd(10.848 + h, 0.0)); - pnts.push_back(Vecd(10.848 - h, 0.0)); - - MultiPolygon multi_polygon; - multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); - return multi_polygon; + std::vector pnts; + pnts.push_back(Vecd(10.848 - h, 0.0)); + pnts.push_back(Vecd(10.848 - h, 1.4)); + pnts.push_back(Vecd(10.848 + h, 1.4)); + pnts.push_back(Vecd(10.848 + h, 0.0)); + pnts.push_back(Vecd(10.848 - h, 0.0)); + + MultiPolygon multi_polygon; + multi_polygon.addAPolygon(pnts, ShapeBooleanOps::add); + return multi_polygon; } \ No newline at end of file diff --git a/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp b/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp index b7a45de58c..a4d27b825c 100644 --- a/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam/oscillating_beam.cpp @@ -135,7 +135,7 @@ int main(int ac, char *av[]) Dynamics1Level stress_relaxation_second_half(beam_body_inner); // clamping a solid body part. This is softer than a direct constraint BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_beam_base(beam_base); + SimpleDynamics constraint_beam_base(beam_base); //----------------------------------------------------------------------------- // outputs //----------------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp index 25e0e30211..1ddd65b40b 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_UL/oscillating_beam_UL.cpp @@ -145,7 +145,7 @@ int main(int ac, char *av[]) ReduceDynamics fluid_acoustic_time_step(beam_body, 0.4); // clamping a solid body part. BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_beam_base(beam_base); + SimpleDynamics constraint_beam_base(beam_base); //----------------------------------------------------------------------------- // outputs //----------------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_oscillating_beam_cauchy/cauchy_oscillating_beam.cpp b/tests/2d_examples/test_2d_oscillating_beam_cauchy/cauchy_oscillating_beam.cpp index 0fb2ab6601..7d6d2c7fb9 100644 --- a/tests/2d_examples/test_2d_oscillating_beam_cauchy/cauchy_oscillating_beam.cpp +++ b/tests/2d_examples/test_2d_oscillating_beam_cauchy/cauchy_oscillating_beam.cpp @@ -136,7 +136,7 @@ int main(int ac, char *av[]) Dynamics1Level stress_relaxation_second_half(beam_body_inner); // clamping a solid body part. This is softer than a direct constraint BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_beam_base(beam_base); + SimpleDynamics constraint_beam_base(beam_base); //----------------------------------------------------------------------------- // outputs //----------------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_owsc/owsc.h b/tests/2d_examples/test_2d_owsc/owsc.h index a35e19a1fe..29dcefcd85 100644 --- a/tests/2d_examples/test_2d_owsc/owsc.h +++ b/tests/2d_examples/test_2d_owsc/owsc.h @@ -258,7 +258,7 @@ class FlapSystemForSimbody : public SolidBodyPartForSimbody } }; -class WaveMaking : public solid_dynamics::MotionConstraint +class WaveMaking : public BodyPartMotionConstraint { Real model_scale_; Real gravity_; @@ -325,9 +325,11 @@ class WaveMaking : public solid_dynamics::MotionConstraint public: WaveMaking(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part), + : BodyPartMotionConstraint(body_part), model_scale_(25.0), gravity_(gravity_g), water_depth_(Water_H), wave_height_(5.0), - wave_period_(10.0) + wave_period_(10.0), + mass_(*particles_->getVariableByName("Mass")), + force_(*particles_->getVariableByName("Force")) { computeWaveStrokeAndFrequency(); } @@ -339,6 +341,10 @@ class WaveMaking : public solid_dynamics::MotionConstraint vel_[index_i] = getVelocity(time); force_[index_i] = mass_[index_i] * getAcceleration(time); }; + + protected: + StdLargeVec &mass_; + StdLargeVec &force_; }; Real h = 1.3 * particle_spacing_ref; diff --git a/tests/2d_examples/test_2d_self_contact/self_contact.cpp b/tests/2d_examples/test_2d_self_contact/self_contact.cpp index 313df6da5a..23284311b3 100644 --- a/tests/2d_examples/test_2d_self_contact/self_contact.cpp +++ b/tests/2d_examples/test_2d_self_contact/self_contact.cpp @@ -146,7 +146,7 @@ int main(int ac, char *av[]) InteractionDynamics beam_self_contact_density(beam_self_contact); InteractionWithUpdate beam_self_contact_forces(beam_self_contact); BodyRegionByParticle beam_base(beam_body, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_beam_base(beam_base); + SimpleDynamics constraint_beam_base(beam_base); //----------------------------------------------------------------------------- // outputs //----------------------------------------------------------------------------- diff --git a/tests/2d_examples/test_2d_shell_beam_collision/shell_beam_collision.cpp b/tests/2d_examples/test_2d_shell_beam_collision/shell_beam_collision.cpp index 050e02a200..1c7692ca1f 100644 --- a/tests/2d_examples/test_2d_shell_beam_collision/shell_beam_collision.cpp +++ b/tests/2d_examples/test_2d_shell_beam_collision/shell_beam_collision.cpp @@ -193,7 +193,7 @@ int main(int ac, char *av[]) InteractionWithUpdate beam_compute_solid_contact_forces(beam_contact); InteractionWithUpdate shell_compute_solid_contact_forces(shell_contact); BodyRegionByParticle holder(beam, makeShared(createBeamConstrainShape())); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); /** Damping with the solid body*/ DampingWithRandomChoice>> beam_damping(0.5, beam_inner, "Velocity", physical_viscosity); diff --git a/tests/2d_examples/test_2d_stretching/stretching.cpp b/tests/2d_examples/test_2d_stretching/stretching.cpp index ed92832e05..fb9dace35c 100644 --- a/tests/2d_examples/test_2d_stretching/stretching.cpp +++ b/tests/2d_examples/test_2d_stretching/stretching.cpp @@ -82,12 +82,12 @@ class Beam : public MultiPolygonShape } }; -class LeftStretchSolidBodyRegion : public solid_dynamics::MotionConstraint +class LeftStretchSolidBodyRegion : public BodyPartMotionConstraint { public: // TODO: use only body part as argment since body can be referred from it already LeftStretchSolidBodyRegion(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part), + : BodyPartMotionConstraint(body_part), vel_(particles_->vel_), pos_(particles_->pos_){}; virtual ~LeftStretchSolidBodyRegion(){}; @@ -101,12 +101,12 @@ class LeftStretchSolidBodyRegion : public solid_dynamics::MotionConstraint }; }; -class RightStretchSolidBodyRegion : public solid_dynamics::MotionConstraint +class RightStretchSolidBodyRegion : public BodyPartMotionConstraint { public: // TODO: use only body part as argment since body can be referred from it already RightStretchSolidBodyRegion(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part), + : BodyPartMotionConstraint(body_part), vel_(particles_->vel_), pos_(particles_->pos_){}; virtual ~RightStretchSolidBodyRegion(){}; @@ -143,12 +143,12 @@ MultiPolygon createConstrainBeamShape() return multi_polygon; }; -class ConstrainXVelocity : public solid_dynamics::MotionConstraint +class ConstrainXVelocity : public BodyPartMotionConstraint { public: // TODO: use only body part as argment since body can be referred from it already ConstrainXVelocity(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part), + : BodyPartMotionConstraint(body_part), vel_(particles_->vel_), pos_(particles_->pos_){}; virtual ~ConstrainXVelocity(){}; diff --git a/tests/2d_examples/test_2d_viscous_cream_drop/viscous_cream_drop.cpp b/tests/2d_examples/test_2d_viscous_cream_drop/viscous_cream_drop.cpp index 640f7d11eb..68b5c38dc0 100644 --- a/tests/2d_examples/test_2d_viscous_cream_drop/viscous_cream_drop.cpp +++ b/tests/2d_examples/test_2d_viscous_cream_drop/viscous_cream_drop.cpp @@ -165,7 +165,7 @@ int main(int ac, char *av[]) Dynamics1Level cream_stress_relaxation_second_half(cream_inner); /** constraint for the cream. */ BodyRegionByParticle platform(cream, makeShared(createPlatformConstraint())); - SimpleDynamics platform_constraint(platform); + SimpleDynamics platform_constraint(platform); //---------------------------------------------------------------------- // Define the methods for I/O operations and observations of the simulation. //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp b/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp index 32d920aadd..83635188f9 100644 --- a/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp +++ b/tests/3d_examples/test_3d_beam_pulling_pressure_load/beam_pulling_pressure_load.cpp @@ -160,7 +160,7 @@ int main(int ac, char *av[]) /* create a brick to tag the region */ Vecd half_size_1(0.03, 0.03, 0.02); BodyRegionByParticle holder(beam_body, makeShared(half_size_1, 1, Vecd(0.0, 0.0, -0.02))); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); /** Damping with the solid body*/ DampingWithRandomChoice>> diff --git a/tests/3d_examples/test_3d_dynamic_plate/test_3d_dynamic_plate.cpp b/tests/3d_examples/test_3d_dynamic_plate/test_3d_dynamic_plate.cpp index a83d3e5eb7..dd55b598d1 100644 --- a/tests/3d_examples/test_3d_dynamic_plate/test_3d_dynamic_plate.cpp +++ b/tests/3d_examples/test_3d_dynamic_plate/test_3d_dynamic_plate.cpp @@ -139,7 +139,7 @@ int main(int ac, char *av[]) stress_relaxation_second_half(plate_body_inner); /** Constrain the Boundary. */ BoundaryGeometry boundary_geometry(plate_body, "BoundaryGeometry"); - SimpleDynamics constrain_holder(boundary_geometry); + SimpleDynamics constrain_holder(boundary_geometry); /** Output */ IOEnvironment io_environment(sph_system); BodyStatesRecordingToVtp write_states(sph_system.real_bodies_); diff --git a/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp b/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp index 29a493a1ac..1b9e740444 100644 --- a/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp +++ b/tests/3d_examples/test_3d_heart_electromechanics/excitation-contraction.cpp @@ -418,7 +418,7 @@ int main(int ac, char *av[]) /** Constrain region of the inserted body. */ MuscleBaseShapeParameters muscle_base_parameters; BodyRegionByParticle muscle_base(mechanics_heart, makeShared(muscle_base_parameters, "Holder")); - SimpleDynamics constraint_holder(muscle_base); + SimpleDynamics constraint_holder(muscle_base); //---------------------------------------------------------------------- // SPH Output section //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.cpp b/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.cpp index 6593cc4758..1a41ec19bc 100644 --- a/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.cpp +++ b/tests/3d_examples/test_3d_heart_volume_change/excitation_contraction.cpp @@ -181,7 +181,7 @@ int main(int ac, char *av[]) /** Constrain region of the inserted body. */ MuscleBaseShapeParameters muscle_base_parameters; BodyRegionByParticle muscle_base(mechanics_heart, makeShared(muscle_base_parameters, "Holder")); - SimpleDynamics constraint_holder(muscle_base); + SimpleDynamics constraint_holder(muscle_base); //---------------------------------------------------------------------- // SPH Output section //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_muscle_activation/src/muscle_activation.cpp b/tests/3d_examples/test_3d_muscle_activation/src/muscle_activation.cpp index 20ac25040b..8ed5cd9662 100644 --- a/tests/3d_examples/test_3d_muscle_activation/src/muscle_activation.cpp +++ b/tests/3d_examples/test_3d_muscle_activation/src/muscle_activation.cpp @@ -85,7 +85,7 @@ int main(int ac, char *av[]) ReduceDynamics computing_time_step_size(muscle_body); SimpleDynamics myocardium_activation(muscle_body); BodyRegionByParticle holder(muscle_body, makeShared>(Transform(translation_holder), halfsize_holder)); - SimpleDynamics constrain_holder(holder, Vecd(0.0, 1.0, 1.0)); + SimpleDynamics constrain_holder(holder, Vecd(0.0, 1.0, 1.0)); //---------------------------------------------------------------------- // Define the methods for I/O operations, observations // and regression tests of the simulation. diff --git a/tests/3d_examples/test_3d_muscle_soft_body_contact/muscle_soft_body_contact.cpp b/tests/3d_examples/test_3d_muscle_soft_body_contact/muscle_soft_body_contact.cpp index 7a20a6c98d..826f2f910e 100644 --- a/tests/3d_examples/test_3d_muscle_soft_body_contact/muscle_soft_body_contact.cpp +++ b/tests/3d_examples/test_3d_muscle_soft_body_contact/muscle_soft_body_contact.cpp @@ -101,7 +101,7 @@ int main(int ac, char *av[]) /** Constrain the holder. */ TransformShape holder_shape(Transform(translation_stationary_plate), halfsize_stationary_plate, "Holder"); BodyRegionByParticle holder(myocardium_body, holder_shape); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); /** Add spring constraint on the plate. */ SimpleDynamics spring_constraint(moving_plate, Vecd(0.2, 0, 0), 0.01); /** Damping with the solid body*/ diff --git a/tests/3d_examples/test_3d_muscle_solid_contact/muscle_solid_contact.cpp b/tests/3d_examples/test_3d_muscle_solid_contact/muscle_solid_contact.cpp index a732fbcf78..29d38a1ce1 100644 --- a/tests/3d_examples/test_3d_muscle_solid_contact/muscle_solid_contact.cpp +++ b/tests/3d_examples/test_3d_muscle_solid_contact/muscle_solid_contact.cpp @@ -97,7 +97,7 @@ int main(int ac, char *av[]) /** Constrain the holder. */ TransformShape holder_shape(Transform(translation_stationary_plate), halfsize_stationary_plate, "Holder"); BodyRegionByParticle holder(myocardium_body, holder_shape); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); /** Damping with the solid body*/ DampingWithRandomChoice>> muscle_damping(0.1, myocardium_body_inner, "Velocity", physical_viscosity); diff --git a/tests/3d_examples/test_3d_nonlinear_wave_fsi/nonlinear_wave_fsi.h b/tests/3d_examples/test_3d_nonlinear_wave_fsi/nonlinear_wave_fsi.h index 1f43427948..66a59a05fe 100644 --- a/tests/3d_examples/test_3d_nonlinear_wave_fsi/nonlinear_wave_fsi.h +++ b/tests/3d_examples/test_3d_nonlinear_wave_fsi/nonlinear_wave_fsi.h @@ -335,7 +335,7 @@ class WallBoundary : public ComplexShape //---------------------------------------------------------------------- // Boundary condition for wave_maker //---------------------------------------------------------------------- -class WaveMaking : public solid_dynamics::MotionConstraint +class WaveMaking : public BodyPartMotionConstraint { Real h; Real tf; @@ -454,8 +454,11 @@ class WaveMaking : public solid_dynamics::MotionConstraint public: WaveMaking(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part), - h(WH), tf(5), xf(4.5), fmn(0.32), fmx(0.96), a(0.0078), N(32), g(gravity_g) + : BodyPartMotionConstraint(body_part), + h(WH), tf(5), xf(4.5), fmn(0.32), fmx(0.96), a(0.0078), N(32), g(gravity_g), + mass_(*particles_->getVariableByName("Mass")), + force_(*particles_->getVariableByName("Force")) + { ComputeWaveChar(); } @@ -467,6 +470,9 @@ class WaveMaking : public solid_dynamics::MotionConstraint vel_[index_i] = getVelocity(time); force_[index_i] = mass_[index_i] * getAcceleration(time); }; + protected: + StdLargeVec &mass_; + StdLargeVec &force_; }; //---------------------------------------------------------------------- diff --git a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp index 3abb877fc7..c8c7a01644 100644 --- a/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp +++ b/tests/3d_examples/test_3d_oscillating_plate_UL/oscillating_plate_UL.cpp @@ -142,9 +142,8 @@ int main(int ac, char *av[]) Dynamics1Level plate_shear_stress_relaxation(plate_body_inner); /** Constrain the Boundary. */ BoundaryGeometry boundary_geometry(plate_body, "BoundaryGeometry"); - SimpleDynamics constrain_holder(boundary_geometry, Vecd(1.0, 1.0, 0.0)); - SimpleDynamics - constrain_mass_center(plate_body, Vecd(1.0, 1.0, 0.0)); + SimpleDynamics constrain_holder(boundary_geometry, Vecd(1.0, 1.0, 0.0)); + SimpleDynamics constrain_mass_center(plate_body, Vecd(1.0, 1.0, 0.0)); /** Output */ BodyStatesRecordingToVtp write_states(sph_system.real_bodies_); RestartIO restart_io(sph_system.real_bodies_); diff --git a/tests/3d_examples/test_3d_passive_cantilever/passive_cantilever.cpp b/tests/3d_examples/test_3d_passive_cantilever/passive_cantilever.cpp index 3590c7e1dc..fb1535b9ee 100644 --- a/tests/3d_examples/test_3d_passive_cantilever/passive_cantilever.cpp +++ b/tests/3d_examples/test_3d_passive_cantilever/passive_cantilever.cpp @@ -101,7 +101,7 @@ int main(int ac, char *av[]) /** Constrain the holder. */ TransformShape holder_shape(Transform(translation_holder), halfsize_holder, "Holder"); BodyRegionByParticle holder(cantilever_body, holder_shape); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); /** Output */ IOEnvironment io_environment(sph_system); BodyStatesRecordingToVtp write_states(sph_system.real_bodies_); diff --git a/tests/3d_examples/test_3d_passive_cantilever_neohookean/passive_cantilever_neohookean.cpp b/tests/3d_examples/test_3d_passive_cantilever_neohookean/passive_cantilever_neohookean.cpp index 4a19e996de..8f3885251c 100644 --- a/tests/3d_examples/test_3d_passive_cantilever_neohookean/passive_cantilever_neohookean.cpp +++ b/tests/3d_examples/test_3d_passive_cantilever_neohookean/passive_cantilever_neohookean.cpp @@ -101,7 +101,7 @@ int main(int ac, char *av[]) /** Constrain the holder. */ TransformShape holder_shape(Transform(translation_holder), halfsize_holder, "Holder"); BodyRegionByParticle holder(cantilever_body, holder_shape); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); DampingWithRandomChoice>> muscle_damping(0.1, cantilever_body_inner, "Velocity", physical_viscosity); /** Output */ diff --git a/tests/3d_examples/test_3d_pkj_lv_electrocontraction/pkj_lv_electrocontraction.cpp b/tests/3d_examples/test_3d_pkj_lv_electrocontraction/pkj_lv_electrocontraction.cpp index 67b257cae5..a7c58a8a13 100644 --- a/tests/3d_examples/test_3d_pkj_lv_electrocontraction/pkj_lv_electrocontraction.cpp +++ b/tests/3d_examples/test_3d_pkj_lv_electrocontraction/pkj_lv_electrocontraction.cpp @@ -213,7 +213,7 @@ int main(int ac, char *av[]) /** Constrain region of the inserted body. */ MuscleBaseShapeParameters muscle_base_parameters; BodyRegionByParticle muscle_base(mechanics_heart, makeShared(muscle_base_parameters, "Holder")); - SimpleDynamics constraint_holder(muscle_base); + SimpleDynamics constraint_holder(muscle_base); /** * Pre-simulation. */ diff --git a/tests/3d_examples/test_3d_roof/3d_roof.cpp b/tests/3d_examples/test_3d_roof/3d_roof.cpp index 2e867d17d3..74d81108ab 100644 --- a/tests/3d_examples/test_3d_roof/3d_roof.cpp +++ b/tests/3d_examples/test_3d_roof/3d_roof.cpp @@ -149,7 +149,7 @@ int main(int ac, char *av[]) Dynamics1Level stress_relaxation_second_half(cylinder_body_inner); BoundaryGeometry boundary_geometry(cylinder_body, "BoundaryGeometry"); - SimpleDynamics constrain_holder(boundary_geometry, Vecd(0.0, 1.0, 0.0)); + SimpleDynamics constrain_holder(boundary_geometry, Vecd(0.0, 1.0, 0.0)); DampingWithRandomChoice>> cylinder_position_damping(0.2, cylinder_body_inner, "Velocity", physical_viscosity); DampingWithRandomChoice>> diff --git a/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp b/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp index 4f9eb10471..471781c468 100644 --- a/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp +++ b/tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp @@ -381,7 +381,7 @@ return_data roof_under_self_weight(Real dp, bool cvt = true, int particle_number }(); constrained_edges.body_part_particles_ = constrained_edge_ids; - SimpleDynamics constrain_holder(constrained_edges, length_vec); + SimpleDynamics constrain_holder(constrained_edges, length_vec); DampingWithRandomChoice>> shell_velocity_damping(0.2, shell_body_inner, "Velocity", physical_viscosity); DampingWithRandomChoice>> diff --git a/tests/3d_examples/test_3d_twisting_column/twisting_column.cpp b/tests/3d_examples/test_3d_twisting_column/twisting_column.cpp index a795a4fc3e..d904cad1b5 100644 --- a/tests/3d_examples/test_3d_twisting_column/twisting_column.cpp +++ b/tests/3d_examples/test_3d_twisting_column/twisting_column.cpp @@ -110,7 +110,7 @@ int main(int ac, char *av[]) Dynamics1Level stress_relaxation_second_half(column_inner); TransformShape holder_shape(Transform(translation_holder), halfsize_holder, "Holder"); BodyRegionByParticle holder(column, holder_shape); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); //---------------------------------------------------------------------- // Define the methods for I/O operations, observations // and regression tests of the simulation. diff --git a/tests/extra_source_and_tests/test_2d_lid_driven_cavity/lid_driven_cavity.cpp b/tests/extra_source_and_tests/test_2d_lid_driven_cavity/lid_driven_cavity.cpp index 741a634170..c011bfcd29 100644 --- a/tests/extra_source_and_tests/test_2d_lid_driven_cavity/lid_driven_cavity.cpp +++ b/tests/extra_source_and_tests/test_2d_lid_driven_cavity/lid_driven_cavity.cpp @@ -70,11 +70,11 @@ class FluidFilling : public ComplexShape } }; -class BoundaryVelocity : public solid_dynamics::MotionConstraint +class BoundaryVelocity : public BodyPartMotionConstraint { public: BoundaryVelocity(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part) {} + : BodyPartMotionConstraint(body_part) {} void update(size_t index_i, Real dt = 0.0) { diff --git a/tests/extra_source_and_tests/test_3d_passive_cantilever_LG/passive_cantilever_LG.cpp b/tests/extra_source_and_tests/test_3d_passive_cantilever_LG/passive_cantilever_LG.cpp index 56a50e7cee..cbfd85348a 100644 --- a/tests/extra_source_and_tests/test_3d_passive_cantilever_LG/passive_cantilever_LG.cpp +++ b/tests/extra_source_and_tests/test_3d_passive_cantilever_LG/passive_cantilever_LG.cpp @@ -103,7 +103,7 @@ int main(int ac, char *av[]) /** Constrain the holder. */ TransformShape holder_shape(Transform(translation_holder), halfsize_holder, "Holder"); BodyRegionByParticle holder(cantilever_body, holder_shape); - SimpleDynamics constraint_holder(holder); + SimpleDynamics constraint_holder(holder); /** Output */ IOEnvironment io_environment(sph_system); BodyStatesRecordingToVtp write_states(sph_system.real_bodies_); diff --git a/tests/unit_tests_src/shared/particle_dynamics/fluid_dynamics/test_2d_velocity_gradient/2d_velocity_gradient.cpp b/tests/unit_tests_src/shared/particle_dynamics/fluid_dynamics/test_2d_velocity_gradient/2d_velocity_gradient.cpp index d226f3cf1b..aabe5da95f 100644 --- a/tests/unit_tests_src/shared/particle_dynamics/fluid_dynamics/test_2d_velocity_gradient/2d_velocity_gradient.cpp +++ b/tests/unit_tests_src/shared/particle_dynamics/fluid_dynamics/test_2d_velocity_gradient/2d_velocity_gradient.cpp @@ -91,11 +91,11 @@ class CouetteFlowInitialCondition } }; -class BoundaryVelocity : public solid_dynamics::MotionConstraint +class BoundaryVelocity : public BodyPartMotionConstraint { public: explicit BoundaryVelocity(BodyPartByParticle &body_part) - : solid_dynamics::MotionConstraint(body_part) {} + : BodyPartMotionConstraint(body_part) {} void update(size_t index_i, Real dt = 0.0) { From 3ad2b4f7ab3eba3c86964e1c673ee6b288b60e9b Mon Sep 17 00:00:00 2001 From: Xiangyu Hu Date: Mon, 8 Apr 2024 16:39:49 +0200 Subject: [PATCH 2/2] a compiling error --- .../particle_dynamics/general_dynamics/general_constraint.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/shared/particle_dynamics/general_dynamics/general_constraint.h b/src/shared/particle_dynamics/general_dynamics/general_constraint.h index 588f3c51ea..ac283af4ba 100644 --- a/src/shared/particle_dynamics/general_dynamics/general_constraint.h +++ b/src/shared/particle_dynamics/general_dynamics/general_constraint.h @@ -72,7 +72,7 @@ class MotionConstraint : public BaseLocalDynamics, public Ge : BaseLocalDynamics(identifier), GeneralDataDelegateSimple(identifier.getSPHBody()), pos_(this->particles_->pos_), - pos0_(*this->particles_->registerSharedVariable( + pos0_(*this->particles_->template registerSharedVariable( "InitialPosition", [&](size_t index_i) { return pos_[index_i]; })), vel_(this->particles_->vel_){};