Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Xiangyu/split shape from body #668

Draft
wants to merge 13 commits into
base: master
Choose a base branch
from
Draft
25 changes: 11 additions & 14 deletions modules/opencascade/opencascade/relax_dynamics_surface.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,22 +10,21 @@ namespace SPH
namespace relax_dynamics
{
//=================================================================================================//
ShapeSurfaceBounding2::ShapeSurfaceBounding2(RealBody &real_body_)
ShapeSurfaceConstraint::ShapeSurfaceConstraint(RealBody &real_body_, Shape &shape)
: LocalDynamics(real_body_),
pos_(particles_->getVariableDataByName<Vecd>("Position"))
{
shape_ = &real_body_.getInitialShape();
}
pos_(particles_->getVariableDataByName<Vecd>("Position")),
shape_(&shape) {}
//=================================================================================================//
void ShapeSurfaceBounding2::update(size_t index_i, Real dt)
void ShapeSurfaceConstraint::update(size_t index_i, Real dt)
{
pos_[index_i] = shape_->findClosestPoint(pos_[index_i]);
}
//=================================================================================================//
RelaxationStepInnerFirstHalf::
RelaxationStepInnerFirstHalf(BaseInnerRelation &inner_relation)
RelaxationStepInnerFirstHalf(BaseInnerRelation &inner_relation, Shape &shape)
: BaseDynamics<void>(), real_body_(inner_relation.real_body_),
inner_relation_(inner_relation), relaxation_acceleration_inner_(inner_relation) {}
inner_relation_(inner_relation),
relaxation_acceleration_inner_(inner_relation, shape) {}
//=================================================================================================//
void RelaxationStepInnerFirstHalf::exec(Real dt)
{
Expand All @@ -36,12 +35,10 @@ void RelaxationStepInnerFirstHalf::exec(Real dt)

//=================================================================================================//
RelaxationStepInnerSecondHalf::
RelaxationStepInnerSecondHalf(BaseInnerRelation &inner_relation)
RelaxationStepInnerSecondHalf(BaseInnerRelation &inner_relation, Shape &shape)
: BaseDynamics<void>(), real_body_(inner_relation.real_body_),
get_time_step_square_(*real_body_), update_particle_position_(*real_body_),
surface_bounding_(*real_body_)
{
}
surface_bounding_(*real_body_, shape) {}
//=================================================================================================//
void RelaxationStepInnerSecondHalf::exec(Real dt)
{
Expand All @@ -51,9 +48,9 @@ void RelaxationStepInnerSecondHalf::exec(Real dt)
}

//=================================================================================================//
SurfaceNormalDirection::SurfaceNormalDirection(SPHBody &sph_body)
SurfaceNormalDirection::SurfaceNormalDirection(SPHBody &sph_body, Shape &shape)
: LocalDynamics(sph_body),
surface_shape_(DynamicCast<SurfaceShape>(this, &sph_body.getInitialShape())),
surface_shape_(DynamicCast<SurfaceShape>(this, &shape)),
pos_(particles_->getVariableDataByName<Vecd>("Position")),
n_(particles_->registerStateVariable<Vecd>("NormalDirection")) {}

Expand Down
16 changes: 8 additions & 8 deletions modules/opencascade/opencascade/relax_dynamics_surface.h
Original file line number Diff line number Diff line change
Expand Up @@ -41,11 +41,11 @@ class SurfaceShape;

namespace relax_dynamics
{
class ShapeSurfaceBounding2 : public LocalDynamics
class ShapeSurfaceConstraint : public LocalDynamics
{
public:
ShapeSurfaceBounding2(RealBody &real_body_);
virtual ~ShapeSurfaceBounding2(){};
ShapeSurfaceConstraint(RealBody &real_body, Shape &shape);
virtual ~ShapeSurfaceConstraint(){};
void update(size_t index_i, Real dt = 0.0);

protected:
Expand All @@ -56,7 +56,7 @@ class ShapeSurfaceBounding2 : public LocalDynamics
class RelaxationStepInnerFirstHalf : public BaseDynamics<void>
{
public:
explicit RelaxationStepInnerFirstHalf(BaseInnerRelation &inner_relation);
explicit RelaxationStepInnerFirstHalf(BaseInnerRelation &inner_relation, Shape &shape);
virtual ~RelaxationStepInnerFirstHalf(){};
virtual void exec(Real dt = 0.0) override;

Expand All @@ -69,16 +69,16 @@ class RelaxationStepInnerFirstHalf : public BaseDynamics<void>
class RelaxationStepInnerSecondHalf : public BaseDynamics<void>
{
public:
explicit RelaxationStepInnerSecondHalf(BaseInnerRelation &inner_relation);
explicit RelaxationStepInnerSecondHalf(BaseInnerRelation &inner_relation, Shape &shape);
virtual ~RelaxationStepInnerSecondHalf(){};
SimpleDynamics<ShapeSurfaceBounding2> &SurfaceBounding() { return surface_bounding_; };
SimpleDynamics<ShapeSurfaceConstraint> &SurfaceBounding() { return surface_bounding_; };
virtual void exec(Real dt = 0.0) override;

protected:
RealBody *real_body_;
ReduceDynamics<RelaxationScaling> get_time_step_square_;
SimpleDynamics<PositionRelaxation> update_particle_position_;
SimpleDynamics<ShapeSurfaceBounding2> surface_bounding_;
SimpleDynamics<ShapeSurfaceConstraint> surface_bounding_;
};

/**
Expand All @@ -88,7 +88,7 @@ class RelaxationStepInnerSecondHalf : public BaseDynamics<void>
class SurfaceNormalDirection : public LocalDynamics
{
public:
explicit SurfaceNormalDirection(SPHBody &sph_body);
explicit SurfaceNormalDirection(SPHBody &sph_body, Shape &shape);
virtual ~SurfaceNormalDirection(){};
void update(size_t index_i, Real dt = 0.0);

Expand Down
32 changes: 17 additions & 15 deletions modules/opencascade/tests/test_3d_aortic_valve/aortic_valve.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -59,12 +59,11 @@ template <>
class ParticleGenerator<SurfaceParticles, Leaflet> : public ParticleGenerator<SurfaceParticles>
{
public:
explicit ParticleGenerator(SPHBody &sph_body, SurfaceParticles &surface_particles)
: ParticleGenerator<SurfaceParticles>(sph_body, surface_particles), sph_body_(sph_body){};
explicit ParticleGenerator(SPHBody &sph_body, SurfaceParticles &surface_particles, Shape &shape)
: ParticleGenerator<SurfaceParticles>(sph_body, surface_particles),
surface_shape_(DynamicCast<SurfaceShape>(this, &shape)){};
virtual void prepareGeometricData() override
{
SurfaceShape *a = DynamicCast<SurfaceShape>(this, &sph_body_.getInitialShape());

Standard_Real u1 = 0;
Standard_Real v1 = DELTA1;
Standard_Real u2 = DELTA2;
Expand All @@ -74,24 +73,24 @@ class ParticleGenerator<SurfaceParticles, Leaflet> : public ParticleGenerator<Su
for (size_t k = 0; k <= 1 / DELTA1; k++)
{
Standard_Real u = u1 + k * DELTA1;
points.push_back(a->getCartesianPoint(u, 0));
points.push_back(surface_shape_->getCartesianPoint(u, 0));
}
for (size_t k = 0; k <= 1 / DELTA1 - 1; k++)
{
Standard_Real v = v1 + k * DELTA1;
points.push_back(a->getCartesianPoint(0, v));
points.push_back(surface_shape_->getCartesianPoint(0, v));
}

for (size_t n = 0; n <= 1 / DELTA2 - 2; n++)
{
Standard_Real u = u2 + n * DELTA2;
points.push_back(a->getCartesianPoint(u, 1));
points.push_back(surface_shape_->getCartesianPoint(u, 1));
}

for (size_t n = 0; n <= 1 / DELTA2 - 1; n++)
{
Standard_Real v = v2 + n * DELTA2;
points.push_back(a->getCartesianPoint(1, v));
points.push_back(surface_shape_->getCartesianPoint(1, v));
}

for (size_t k = 0; k <= 8; k++)
Expand All @@ -104,7 +103,7 @@ class ParticleGenerator<SurfaceParticles, Leaflet> : public ParticleGenerator<Su
for (size_t n = 0; n <= 20; n++)
{
Standard_Real u = u1 + n * U_DELTA;
points.push_back(a->getCartesianPoint(u, v));
points.push_back(surface_shape_->getCartesianPoint(u, v));
}
}

Expand All @@ -115,7 +114,9 @@ class ParticleGenerator<SurfaceParticles, Leaflet> : public ParticleGenerator<Su
addSurfaceProperties(n_0, thickness);
}
}
SPHBody &sph_body_;

protected:
SurfaceShape *surface_shape_;
};
} // namespace SPH

Expand All @@ -132,10 +133,11 @@ int main(int ac, char *av[])
//----------------------------------------------------------------------
// Creating body, materials and particles.
//----------------------------------------------------------------------
SolidBody leaflet(sph_system, makeShared<SurfaceShapeSTEP>(full_path_to_geometry, "Leaflet"));
SurfaceShapeSTEP leaflet_shape(full_path_to_geometry, "Leaflet");
SolidBody leaflet(sph_system, leaflet_shape.getName());
// here dummy linear elastic solid is use because no solid dynamics in particle relaxation
leaflet.defineMaterial<Solid>();
leaflet.generateParticles<SurfaceParticles, Leaflet>();
leaflet.generateParticles<SurfaceParticles, Leaflet>(leaflet_shape);
//----------------------------------------------------------------------
// Define simple file input and outputs functions.
//----------------------------------------------------------------------
Expand All @@ -154,11 +156,11 @@ int main(int ac, char *av[])
// Methods used for particle relaxation.
//----------------------------------------------------------------------
using namespace relax_dynamics;
RelaxationStepInnerFirstHalf leaflet_relaxation_first_half(leaflet_inner);
RelaxationStepInnerSecondHalf leaflet_relaxation_second_half(leaflet_inner);
RelaxationStepInnerFirstHalf leaflet_relaxation_first_half(leaflet_inner, leaflet_shape);
RelaxationStepInnerSecondHalf leaflet_relaxation_second_half(leaflet_inner, leaflet_shape);
/** Constrain the boundary. */
BoundaryGeometry boundary_geometry(leaflet, "BoundaryGeometry");
SimpleDynamics<SurfaceNormalDirection> surface_normal_direction(leaflet);
SimpleDynamics<SurfaceNormalDirection> surface_normal_direction(leaflet, leaflet_shape);
//----------------------------------------------------------------------
// Particle relaxation starts here.
//----------------------------------------------------------------------
Expand Down
44 changes: 26 additions & 18 deletions modules/structural_simulation/structural_simulation_class.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -23,20 +23,20 @@ BodyPartFromMesh::BodyPartFromMesh(SPHBody &body, SharedPtr<TriangleMeshShape> t
SolidBodyFromMesh::SolidBodyFromMesh(
SPHSystem &system, SharedPtr<TriangleMeshShape> triangle_mesh_shape, Real resolution,
SharedPtr<SaintVenantKirchhoffSolid> material_model, Vecd *pos_0, Real *volume)
: SolidBody(system, triangle_mesh_shape)
: SolidBody(system, triangle_mesh_shape->getName())
{
defineAdaptationRatios(1.15, system.resolution_ref_ / resolution);
defineBodyLevelSetShape()->cleanLevelSet();
LevelSetShape level_set_shape(*this, *triangle_mesh_shape, 1.0);
assignMaterial(material_model.get());
generateParticles<BaseParticles, Lattice>();
generateParticles<BaseParticles, Lattice>(*level_set_shape.cleanLevelSet());
}

SolidBodyForSimulation::SolidBodyForSimulation(
SPHSystem &system, SharedPtr<TriangleMeshShape> triangle_mesh_shape, Real resolution,
Real physical_viscosity, SharedPtr<SaintVenantKirchhoffSolid> material_model, Vecd *pos_0, Real *volume)
: solid_body_from_mesh_(system, triangle_mesh_shape, resolution, material_model, pos_0, volume),
inner_body_relation_(solid_body_from_mesh_),
initial_normal_direction_(SimpleDynamics<NormalDirectionFromBodyShape>(solid_body_from_mesh_)),
initial_normal_direction_(SimpleDynamics<NormalDirectionFromBodyShape>(solid_body_from_mesh_, *triangle_mesh_shape)),
correct_configuration_(inner_body_relation_),
stress_relaxation_first_half_(inner_body_relation_),
stress_relaxation_second_half_(inner_body_relation_),
Expand All @@ -62,17 +62,18 @@ void expandBoundingBox(BoundingBox *original, BoundingBox *additional)
}

void relaxParticlesSingleResolution(bool write_particle_relaxation_data,
SolidBody &solid_body_from_mesh,
InnerRelation &solid_body_from_mesh_inner)
InnerRelation &solid_body_from_mesh_inner,
LevelSetShape *level_set_shape)
{
RealBody &solid_body_from_mesh = dynamic_cast<RealBody &>(solid_body_from_mesh_inner.getSPHBody());
BodyStatesRecordingToVtp write_solid_body_from_mesh_to_vtp(solid_body_from_mesh);

//----------------------------------------------------------------------
// Methods used for particle relaxation.
//----------------------------------------------------------------------
SimpleDynamics<relax_dynamics::RandomizeParticlePosition> random_solid_body_from_mesh_particles(solid_body_from_mesh);
/** A Physics relaxation step. */
relax_dynamics::RelaxationStepLevelSetCorrectionInner relaxation_step_inner(solid_body_from_mesh_inner);
relax_dynamics::RelaxationStepLevelSetCorrectionInner relaxation_step_inner(solid_body_from_mesh_inner, level_set_shape);
//----------------------------------------------------------------------
// Particle relaxation starts here.
//----------------------------------------------------------------------
Expand Down Expand Up @@ -108,16 +109,16 @@ std::tuple<Vecd *, Real *> generateAndRelaxParticlesFromMesh(
{
BoundingBox bb = triangle_mesh_shape->getBounds();
SPHSystem system(bb, resolution);
SolidBody model(system, triangle_mesh_shape);
model.defineBodyLevelSetShape()->cleanLevelSet();
SolidBody model(system, triangle_mesh_shape->getName());
LevelSetShape level_set_shape(model, *triangle_mesh_shape, 1.0);
model.defineMaterial<Solid>();
model.generateParticles<BaseParticles, Lattice>();
model.generateParticles<BaseParticles, Lattice>(*level_set_shape.cleanLevelSet());

if (particle_relaxation)
{
system.setIOEnvironment();
InnerRelation inner_relation(model);
relaxParticlesSingleResolution(write_particle_relaxation_data, model, inner_relation);
relaxParticlesSingleResolution(write_particle_relaxation_data, inner_relation, &level_set_shape);
}

return std::tuple<Vecd *, Real *>(model.getBaseParticles().ParticlePositions(), model.getBaseParticles().VolumetricMeasures());
Expand Down Expand Up @@ -345,10 +346,12 @@ void StructuralSimulation::initializeElasticSolidBodies()
void StructuralSimulation::initializeContactBetweenTwoBodies(int first, int second)
{
SolidBodyFromMesh *first_body = solid_body_list_[first]->getSolidBodyFromMesh();
Shape &first_shape = *body_mesh_list_[first];
SolidBodyFromMesh *second_body = solid_body_list_[second]->getSolidBodyFromMesh();
Shape &second_shape = *body_mesh_list_[second];

contact_list_.emplace_back(makeShared<SurfaceContactRelation>(*first_body, RealBodyVector({second_body})));
contact_list_.emplace_back(makeShared<SurfaceContactRelation>(*second_body, RealBodyVector({first_body})));
contact_list_.emplace_back(makeShared<SurfaceContactRelation>(*first_body, first_shape, RealBodyVector({second_body})));
contact_list_.emplace_back(makeShared<SurfaceContactRelation>(*second_body, second_shape, RealBodyVector({first_body})));

int last = contact_list_.size() - 1;
contact_density_list_.push_back(makeShared<InteractionDynamics<solid_dynamics::ContactFactorSummation>>(*contact_list_[last - 1]));
Expand All @@ -374,7 +377,7 @@ void StructuralSimulation::initializeAllContacts()
target_list.emplace_back(solid_body_list_[target_i]->getSolidBodyFromMesh());
}

contact_list_.emplace_back(makeShared<SurfaceContactRelation>(*contact_body, target_list));
contact_list_.emplace_back(makeShared<SurfaceContactRelation>(*contact_body, *body_mesh_list_[i], target_list));
int last = contact_list_.size() - 1;
contact_density_list_.emplace_back(makeShared<InteractionDynamics<solid_dynamics::ContactFactorSummation>>(*contact_list_[last]));
contact_force_list_.emplace_back(makeShared<InteractionWithUpdate<solid_dynamics::ContactForce>>(*contact_list_[last]));
Expand Down Expand Up @@ -461,7 +464,8 @@ void StructuralSimulation::initializeSurfacePressure()
StdVec<std::array<Real, 2>> pressure_over_time = std::get<3>(surface_pressure_tuple_[i]);

BodyPartByParticle *bp = body_part_tri_mesh_ptr_keeper_.createPtr<BodyPartFromMesh>(*solid_body_list_[body_index]->getSolidBodyFromMesh(), tri_mesh);
surface_pressure_.emplace_back(makeShared<SimpleDynamics<solid_dynamics::SurfacePressureFromSource>>(*bp, point, pressure_over_time));
surface_pressure_.emplace_back(makeShared<SimpleDynamics<solid_dynamics::SurfacePressureFromSource>>(
*bp, *body_mesh_list_[i], point, pressure_over_time));
}
}

Expand Down Expand Up @@ -489,7 +493,7 @@ void StructuralSimulation::initializeSpringNormalOnSurfaceParticles()
Real damping_coefficient = std::get<5>(surface_spring_tuple_[i]);

surface_spring_.emplace_back(makeShared<SimpleDynamics<solid_dynamics::SpringNormalOnSurfaceParticles>>(
*solid_body_list_[body_index]->getSolidBodyFromMesh(),
*solid_body_list_[body_index]->getSolidBodyFromMesh(), *body_mesh_list_[body_index],
inner_outer, point, spring_stiffness, damping_coefficient));
}
}
Expand Down Expand Up @@ -538,7 +542,9 @@ void StructuralSimulation::initializePositionSolidBody()
Real end_time = std::get<2>(position_solid_body_tuple_[i]);
Vecd pos_end_center = std::get<3>(position_solid_body_tuple_[i]);
position_solid_body_.emplace_back(makeShared<SimpleDynamics<solid_dynamics::PositionSolidBody>>(
*solid_body_list_[body_index]->getSolidBodyFromMesh(), start_time, end_time, pos_end_center));
*solid_body_list_[body_index]->getSolidBodyFromMesh(),
body_mesh_list_[body_index]->getBounds(),
start_time, end_time, pos_end_center));
}
}

Expand All @@ -552,7 +558,9 @@ void StructuralSimulation::initializePositionScaleSolidBody()
Real end_time = std::get<2>(position_scale_solid_body_tuple_[i]);
Real scale = std::get<3>(position_scale_solid_body_tuple_[i]);
position_scale_solid_body_.emplace_back(makeShared<SimpleDynamics<solid_dynamics::PositionScaleSolidBody>>(
*solid_body_list_[body_index]->getSolidBodyFromMesh(), start_time, end_time, scale));
*solid_body_list_[body_index]->getSolidBodyFromMesh(),
body_mesh_list_[body_index]->getBounds(),
start_time, end_time, scale));
}
}

Expand Down
5 changes: 0 additions & 5 deletions modules/structural_simulation/structural_simulation_class.h
Original file line number Diff line number Diff line change
Expand Up @@ -116,11 +116,6 @@ class SolidBodyForSimulation

void expandBoundingBox(BoundingBox *original, BoundingBox *additional);

void relaxParticlesSingleResolution(bool write_particles_to_file,
SolidBodyFromMesh &solid_body_from_mesh,
BaseParticles &solid_body_from_mesh_particles,
InnerRelation &solid_body_from_mesh_inner);

static inline Real getPhysicalViscosityGeneral(Real rho, Real youngs_modulus, Real length_scale, Real shape_constant = 1.0)
{
// the physical viscosity is defined in the paper pf prof. Hu
Expand Down
Loading
Loading