all compiled, to test

tests/2d_examples/test_2d_two_phase_dambreak/two_phase_dambreak.cpp

@@ -60,8 +60,8 @@ int main(int ac, char *av[])
     //	Define the main numerical methods used in the simulation.
     //	Note that there may be data dependence on the constructors of these methods.
     //----------------------------------------------------------------------
-    /** Initialize particle acceleration. */
-    SimpleDynamics<NormalDirectionFromSubShapeAndOp> inner_normal_direction(wall_boundary, "InnerWall");
+    SubShapeAndOp *inner_wall_shape_and_op = wall_boundary_shape.getSubShapeAndOpByName("InnerWall");
+    SimpleDynamics<NormalDirectionFromSubShapeAndOp> inner_normal_direction(wall_boundary, inner_wall_shape_and_op);
 
     Gravity gravity(Vecd(0.0, -gravity_g));
     SimpleDynamics<GravityForce<Gravity>> constant_gravity_to_water(water_block, gravity);

tests/3d_examples/test_3d_incompressible_channel_flow/test_3d_incompressible_channel_flow.cpp

@@ -22,7 +22,7 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     //	Creating body, materials and particles.
     //----------------------------------------------------------------------
-    FluidBody air_block(sph_system, makeShared<AirBody>("AirBody"));
+    FluidBody air_block(sph_system, "AirBody");
     air_block.defineMaterial<WeaklyCompressibleFluid>(rho0_f, c_f, mu_f);
     Ghost<ReserveSizeFactor> ghost_boundary(0.5);
     air_block.generateParticlesWithReserve<BaseParticles, UnstructuredMesh>(ghost_boundary, read_mesh_data);

tests/3d_examples/test_3d_incompressible_channel_flow/test_3d_incompressible_channel_flow.h

@@ -33,20 +33,6 @@ Real mu_f = 0.0; /**< Dynamics viscosity. */
 //----------------------------------------------------------------------
 std::string mesh_fullpath = "./input/Channel_ICEM.msh";
 //----------------------------------------------------------------------
-//	Define geometries and body shapes
-//----------------------------------------------------------------------
-
-class AirBody : public ComplexShape
-{
-  public:
-    explicit AirBody(const std::string &shape_name) : ComplexShape(shape_name)
-    {
-        Vecd halfsize_wave(0.5 * DH, 0.5 * DL, 0.5 * DW);
-        Transform translation_wave(halfsize_wave);
-        add<TransformShape<GeometricShapeBox>>(Transform(translation_wave), halfsize_wave);
-    }
-};
-///----------------------------------------------------------------------
 //	Initialization
 //----------------------------------------------------------------------
 class InvCFInitialCondition

tests/3d_examples/test_3d_particle_relaxation_single_resolution/particle_relaxation_single_resolution_3D.cpp

@@ -98,7 +98,7 @@ int main(int ac, char *av[])
     // level set shape is used for particle relaxation
     LevelSetShape level_set_shape(imported_model, import_model_shape);
     level_set_shape.correctLevelSetSign()->writeLevelSet(sph_system);
-    imported_model.generateParticles<BaseParticles, Lattice>();
+    imported_model.generateParticles<BaseParticles, Lattice>(level_set_shape);
     //----------------------------------------------------------------------
     //	Define simple file input and outputs functions.
     //----------------------------------------------------------------------

tests/3d_examples/test_3d_pkj_lv_electrocontraction/pkj_lv_electrocontraction.cpp

@@ -34,23 +34,22 @@ int main(int ac, char *av[])
     IOEnvironment io_environment(sph_system);
 
     Heart triangle_mesh_heart_model("HeartModel");
-    SharedPtr<LevelSetShape> level_set_heart_model =
-        makeShared<LevelSetShape>(triangle_mesh_heart_model, makeShared<SPHAdaptation>(dp_0));
-    level_set_heart_model->correctLevelSetSign()->writeLevelSet(sph_system);
+    LevelSetShape level_set_heart_model(triangle_mesh_heart_model, makeShared<SPHAdaptation>(dp_0));
+    level_set_heart_model.correctLevelSetSign()->writeLevelSet(sph_system);
     //----------------------------------------------------------------------
     //	SPH Particle relaxation section
     //----------------------------------------------------------------------
     /** check whether run particle relaxation for body fitted particle distribution. */
     if (sph_system.RunParticleRelaxation() && !sph_system.ReloadParticles())
     {
-        SolidBody herat_model(sph_system, level_set_heart_model);
+        SolidBody herat_model(sph_system, level_set_heart_model.getName());
         herat_model.defineMaterial<LocallyOrthotropicMuscle>(rho0_s, bulk_modulus, fiber_direction, sheet_direction, a0, b0);
-        herat_model.generateParticles<BaseParticles, Lattice>();
+        herat_model.generateParticles<BaseParticles, Lattice>(level_set_heart_model);
         /** topology */
         InnerRelation herat_model_inner(herat_model);
         using namespace relax_dynamics;
         SimpleDynamics<RandomizeParticlePosition> random_particles(herat_model);
-        RelaxationStepInner relaxation_step_inner(herat_model_inner);
+        RelaxationStepInner relaxation_step_inner(herat_model_inner, level_set_heart_model);
         BodyStatesRecordingToVtp write_herat_model_state_to_vtp({herat_model});
         ReloadParticleIO write_particle_reload_files(herat_model);
         write_particle_reload_files.addToReload<Vecd>(herat_model, "Fiber");
@@ -82,9 +81,9 @@ int main(int ac, char *av[])
         IsotropicDiffusion diffusion("Phi", "Phi", diffusion_coeff);
         GetDiffusionTimeStepSize get_time_step_size(herat_model, diffusion);
         FiberDirectionDiffusionRelaxation diffusion_relaxation(herat_model_inner, &diffusion);
-        SimpleDynamics<ComputeFiberAndSheetDirections> compute_fiber_sheet(herat_model, "Phi");
-        BodySurface surface_part(herat_model);
-        SimpleDynamics<DiffusionBCs> impose_diffusion_bc(surface_part, "Phi");
+        SimpleDynamics<ComputeFiberAndSheetDirections> compute_fiber_sheet(herat_model, level_set_heart_model, "Phi");
+        BodySurface surface_part(herat_model, level_set_heart_model);
+        SimpleDynamics<DiffusionBCs> impose_diffusion_bc(surface_part, level_set_heart_model, "Phi");
         impose_diffusion_bc.exec();
         write_herat_model_state_to_vtp.addToWrite<Real>(herat_model, "Phi");
         write_herat_model_state_to_vtp.writeToFile(ite);
@@ -114,28 +113,28 @@ int main(int ac, char *av[])
     //	SPH simulation section
     //----------------------------------------------------------------------
     /** create a SPH body, material and particles */
-    SolidBody physiology_heart(sph_system, level_set_heart_model, "PhysiologyHeart");
+    SolidBody physiology_heart(sph_system, "PhysiologyHeart");
     AlievPanfilowModel aliev_panfilow_model(k_a, c_m, k, a, b, mu_1, mu_2, epsilon);
     MonoFieldElectroPhysiology<LocalDirectionalDiffusion> *myocardium_physiology =
         physiology_heart.defineMaterial<MonoFieldElectroPhysiology<LocalDirectionalDiffusion>>(
             aliev_panfilow_model, diffusion_coeff, bias_coeff, fiber_direction);
     (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles())
         ? physiology_heart.generateParticles<BaseParticles, Reload>("HeartModel")
-        : physiology_heart.generateParticles<BaseParticles, Lattice>();
+        : physiology_heart.generateParticles<BaseParticles, Lattice>(level_set_heart_model);
 
     /** create a SPH body, material and particles */
-    SolidBody mechanics_heart(sph_system, level_set_heart_model, "MechanicalHeart");
+    SolidBody mechanics_heart(sph_system, "MechanicalHeart");
     mechanics_heart.defineMaterial<ActiveMuscle<LocallyOrthotropicMuscle>>(rho0_s, bulk_modulus, fiber_direction, sheet_direction, a0, b0);
     (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles())
         ? mechanics_heart.generateParticles<BaseParticles, Reload>("HeartModel")
-        : mechanics_heart.generateParticles<BaseParticles, Lattice>();
+        : mechanics_heart.generateParticles<BaseParticles, Lattice>(level_set_heart_model);
 
     /** Creat a Purkinje network for fast diffusion, material and particles */
-    TreeBody pkj_body(sph_system, level_set_heart_model, "Purkinje");
+    TreeBody pkj_body(sph_system, "Purkinje");
     MonoFieldElectroPhysiology<IsotropicDiffusion> *pkj_physiology =
         pkj_body.defineMaterial<MonoFieldElectroPhysiology<IsotropicDiffusion>>(
             aliev_panfilow_model, diffusion_coeff * acceleration_factor);
-    pkj_body.generateParticles<BaseParticles, NetworkWithExtraCheck>(starting_point, second_point, 50, 1.0);
+    pkj_body.generateParticles<BaseParticles, NetworkWithExtraCheck>(level_set_heart_model, starting_point, second_point, 50, 1.0);
     TreeTerminates pkj_leaves(pkj_body);
     //----------------------------------------------------------------------
     //	SPH Observation section

tests/3d_examples/test_3d_pkj_lv_electrocontraction/pkj_lv_electrocontraction.h

@@ -290,9 +290,9 @@ class ParticleGenerator<BaseParticles, NetworkWithExtraCheck>
     };
 
   public:
-    ParticleGenerator(SPHBody &sph_body, BaseParticles &base_particles,
+    ParticleGenerator(SPHBody &sph_body, BaseParticles &base_particles, LevelSetShape &level_set_shape,
                       Vecd starting_pnt, Vecd second_pnt, int iterator, Real grad_factor)
         : ParticleGenerator<BaseParticles, Network>(
-              sph_body, base_particles, starting_pnt, second_pnt, iterator, grad_factor){};
+              sph_body, base_particles, level_set_shape, starting_pnt, second_pnt, iterator, grad_factor){};
 };
 } // namespace SPH

tests/3d_examples/test_3d_poiseuille_flow_shell/poiseuille_flow_shell.cpp

@@ -160,10 +160,9 @@ void poiseuille_flow(const Real resolution_ref, const Real resolution_shell, con
     //----------------------------------------------------------------------
     //  Define water shape
     //----------------------------------------------------------------------
-    auto water_block_shape = makeShared<ComplexShape>("WaterBody");
-    water_block_shape->add<TriangleMeshShapeCylinder>(SimTK::UnitVec3(0., 1., 0.), fluid_radius,
-                                                      full_length * 0.5, SimTK_resolution,
-                                                      translation_fluid);
+    TriangleMeshShapeCylinder water_block_shape(SimTK::UnitVec3(0., 1., 0.), fluid_radius,
+                                                full_length * 0.5, SimTK_resolution,
+                                                translation_fluid);
 
     //----------------------------------------------------------------------
     //  Build up -- a SPHSystem --
@@ -174,7 +173,7 @@ void poiseuille_flow(const Real resolution_ref, const Real resolution_shell, con
     //----------------------------------------------------------------------
     //	Creating bodies with corresponding materials and particles.
     //----------------------------------------------------------------------
-    FluidBody water_block(system, water_block_shape);
+    FluidBody water_block(system, "WaterBody");
     water_block.defineMaterial<WeaklyCompressibleFluid>(rho0_f, c_f, mu_f);
     ParticleBuffer<ReserveSizeFactor> inlet_particle_buffer(0.5);
     water_block.generateParticlesWithReserve<BaseParticles, Lattice>(inlet_particle_buffer, water_block_shape);

tests/3d_examples/test_3d_repose_angle/repose_angle.cpp

@@ -84,12 +84,14 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     //	Creating bodies with corresponding materials and particles.
     //----------------------------------------------------------------------
-    RealBody soil_block(sph_system, makeShared<SoilBlock>("GranularBody"));
-    soil_block.defineBodyLevelSetShape()->writeLevelSet(sph_system);
+    SoilBlock soil_block_shape("GranularBody");
+    RealBody soil_block(sph_system, soil_block_shape.getName());
+    LevelSetShape level_set_shape(soil_block, soil_block_shape);
+    level_set_shape.writeLevelSet(sph_system);
     soil_block.defineMaterial<PlasticContinuum>(rho0_s, c_s, Youngs_modulus, poisson, friction_angle);
     (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles())
         ? soil_block.generateParticles<BaseParticles, Reload>(soil_block.getName())
-        : soil_block.generateParticles<BaseParticles, Lattice>();
+        : soil_block.generateParticles<BaseParticles, Lattice>(level_set_shape);
 
     WallBoundary wall_boundary_shape("WallBoundary");
     SolidBody wall_boundary(sph_system, wall_boundary_shape.getName());
@@ -120,7 +122,7 @@ int main(int ac, char *av[])
         //----------------------------------------------------------------------
         using namespace relax_dynamics;
         SimpleDynamics<RandomizeParticlePosition> random_column_particles(soil_block);
-        RelaxationStepInner relaxation_step_inner(soil_block_inner);
+        RelaxationStepInner relaxation_step_inner(soil_block_inner, level_set_shape);
         //----------------------------------------------------------------------
         //	Output for particle relaxation.
         //----------------------------------------------------------------------

tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp

@@ -13,37 +13,6 @@ using namespace SPH;
 
 Real to_rad(Real angle) { return angle * Pi / 180; }
 
-void relax_shell(RealBody &plate_body, Real thickness)
-{
-    // BUG: apparently only works if dp > thickness, otherwise ShellNormalDirectionPrediction::correctNormalDirection() throws error
-    using namespace relax_dynamics;
-    InnerRelation imported_model_inner(plate_body);
-    SimpleDynamics<RandomizeParticlePosition> random_imported_model_particles(plate_body);
-    ShellRelaxationStep relaxation_step_inner(imported_model_inner);
-    ShellNormalDirectionPrediction shell_normal_prediction(imported_model_inner, thickness);
-    //----------------------------------------------------------------------
-    //	Particle relaxation starts here.
-    //----------------------------------------------------------------------
-    random_imported_model_particles.exec(0.25);
-    relaxation_step_inner.MidSurfaceBounding().exec();
-    plate_body.updateCellLinkedList();
-    //----------------------------------------------------------------------
-    //	Particle relaxation time stepping start here.
-    //----------------------------------------------------------------------
-    int ite_p = 0;
-    while (ite_p < 1000)
-    {
-        if (ite_p % 100 == 0)
-        {
-            std::cout << std::fixed << std::setprecision(9) << "Relaxation steps for the inserted body N = " << ite_p << "\n";
-        }
-        relaxation_step_inner.exec();
-        ite_p += 1;
-    }
-    shell_normal_prediction.exec();
-    std::cout << "The physics relaxation process of imported model finish !" << std::endl;
-}
-
 namespace SPH
 {
 class ShellRoof;
@@ -331,13 +300,10 @@ return_data roof_under_self_weight(Real dp, bool cvt = true, int particle_number
         bb_system = get_particles_bounding_box(obj_vertices); // store this
     }
 
-    // shell
-    auto shell_shape = makeShared<ComplexShape>("shell_shape" + std::to_string(dp_cm)); // keep all data for parameter study
-
     // starting the actual simulation
     SPHSystem system(bb_system, dp);
     system.setIOEnvironment(false);
-    SolidBody shell_body(system, shell_shape);
+    SolidBody shell_body(system, "shell_shape" + std::to_string(dp_cm));
     shell_body.defineMaterial<LinearElasticSolid>(rho, E, mu);
     if (cvt)
     {
@@ -354,10 +320,6 @@ return_data roof_under_self_weight(Real dp, bool cvt = true, int particle_number
     system.system_domain_bounds_ = bb_system;
     std::cout << "bb_system.first_: " << bb_system.first_ << std::endl;
     std::cout << "bb_system.second_: " << bb_system.second_ << std::endl;
-    { // recalculate the volume/area after knowing the particle positions
-      // for (auto& vol: shell_particles->Vol_) vol = total_area / shell_particles->TotalRealParticles();
-      // for (auto& mass: shell_particles->mass_) mass = total_area*rho / shell_particles->TotalRealParticles();
-    }
 
     // methods
     InnerRelation shell_body_inner(shell_body);

tests/3d_examples/test_3d_self_contact/3d_self_contact.cpp

@@ -69,16 +69,19 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     //	Creating body, materials and particles.
     //----------------------------------------------------------------------
-    SolidBody coil(sph_system, makeShared<Coil>("Coil"));
-    coil.defineBodyLevelSetShape()->writeLevelSet(sph_system);
+    Coil coil_shape("Coil");
+    SolidBody coil(sph_system, coil_shape.getName());
+    LevelSetShape level_set_shape(coil, coil_shape);
+    level_set_shape.writeLevelSet(sph_system);
     coil.defineMaterial<NeoHookeanSolid>(rho0_s, Youngs_modulus, poisson);
     (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles())
         ? coil.generateParticles<BaseParticles, Reload>(coil.getName())
-        : coil.generateParticles<BaseParticles, Lattice>();
+        : coil.generateParticles<BaseParticles, Lattice>(level_set_shape);
 
-    SolidBody stationary_plate(sph_system, makeShared<StationaryPlate>("StationaryPlate"));
+    StationaryPlate stationary_plate_shape("StationaryPlate");
+    SolidBody stationary_plate(sph_system, stationary_plate_shape.getName());
     stationary_plate.defineMaterial<SaintVenantKirchhoffSolid>(rho0_s, Youngs_modulus, poisson);
-    stationary_plate.generateParticles<BaseParticles, Lattice>();
+    stationary_plate.generateParticles<BaseParticles, Lattice>(stationary_plate_shape);
     //----------------------------------------------------------------------
     //	Define body relation map.
     //	The contact map gives the topological connections between the bodies.
@@ -88,8 +91,8 @@ int main(int ac, char *av[])
     //  inner and contact relations.
     //----------------------------------------------------------------------
     InnerRelation coil_inner(coil);
-    SelfSurfaceContactRelation coil_self_contact(coil);
-    SurfaceContactRelation coil_contact(coil_self_contact, {&stationary_plate});
+    SelfSurfaceContactRelation coil_self_contact(coil, level_set_shape);
+    SurfaceContactRelation coil_contact(coil_self_contact, level_set_shape, {&stationary_plate});
     //----------------------------------------------------------------------
     //	Define simple file input and outputs functions.
     //----------------------------------------------------------------------
@@ -107,7 +110,7 @@ int main(int ac, char *av[])
         // Write the particle reload files.
         ReloadParticleIO write_particle_reload_files(coil);
         // A  Physics relaxation step.
-        RelaxationStepInner relaxation_step_inner(coil_inner);
+        RelaxationStepInner relaxation_step_inner(coil_inner, level_set_shape);
         //----------------------------------------------------------------------
         //	Particle relaxation starts here.
         //----------------------------------------------------------------------

tests/3d_examples/test_3d_shell_particle_relaxation/test_3d_shell_particle_relaxation.cpp

@@ -48,9 +48,11 @@ int main(int ac, char *av[])
     //----------------------------------------------------------------------
     //	Creating body, materials and particles.
     //----------------------------------------------------------------------
-    RealBody imported_model(sph_system, makeShared<ImportedShellModel>("ImportedShellModel"));
-    imported_model.defineBodyLevelSetShape(level_set_refinement_ratio)->correctLevelSetSign()->writeLevelSet(sph_system);
-    imported_model.generateParticles<SurfaceParticles, Lattice>(thickness);
+    ImportedShellModel imported_model_shape("ImportedShellModel");
+    RealBody imported_model(sph_system, imported_model_shape.getName());
+    LevelSetShape level_set_shape(imported_model, imported_model_shape, level_set_refinement_ratio);
+    level_set_shape.correctLevelSetSign()->writeLevelSet(sph_system);
+    imported_model.generateParticles<SurfaceParticles, Lattice>(level_set_shape, thickness);
     //----------------------------------------------------------------------
     //	Define simple file input and outputs functions.
     //----------------------------------------------------------------------
@@ -72,8 +74,8 @@ int main(int ac, char *av[])
     using namespace relax_dynamics;
     SimpleDynamics<RandomizeParticlePosition> random_imported_model_particles(imported_model);
     /** A  Physics relaxation step. */
-    ShellRelaxationStep relaxation_step_inner(imported_model_inner);
-    ShellNormalDirectionPrediction shell_normal_prediction(imported_model_inner, thickness);
+    ShellRelaxationStep relaxation_step_inner(imported_model_inner, level_set_shape);
+    ShellNormalDirectionPrediction shell_normal_prediction(imported_model_inner, level_set_shape, thickness);
     //----------------------------------------------------------------------
     //	Particle relaxation starts here.
     //----------------------------------------------------------------------