From e9d4ac65616acee87e552b679520a208783a69d7 Mon Sep 17 00:00:00 2001 From: rprat-pro Date: Sun, 24 Mar 2024 17:54:07 +0100 Subject: [PATCH 1/4] add stl mesh with polyhedra + examples --- data/config/config_move_particles.msp | 9 +- data/config/config_polyhedra.msp | 14 + data/config/config_spheres.msp | 11 + example/polyhedra/CMakeLists.txt | 2 +- example/polyhedra/stl_mesh/CMakeLists.txt | 19 + example/polyhedra/stl_mesh/alpha3.shp | 25 + example/polyhedra/stl_mesh/box.stl | 73 ++ example/polyhedra/stl_mesh/box_for_octa.stl | 73 ++ example/polyhedra/stl_mesh/octahedron.shp | 44 + .../stl_mesh/stl_mesh_box_hexapod.msp | 57 ++ .../stl_mesh/stl_mesh_box_hexapod_mini.msp | 56 + .../stl_mesh/stl_mesh_box_octahedron.msp | 57 ++ .../stl_mesh/stl_mesh_box_octahedron_mini.msp | 57 ++ src/driver/CMakeLists.txt | 2 +- src/driver/add_stl_mesh.cpp | 75 ++ src/driver/include/exaDEM/driver_base.h | 5 + src/driver/include/exaDEM/driver_stl_mesh.h | 78 ++ src/driver/include/exaDEM/drivers.h | 11 +- src/driver/update_grid_stl_mesh.cpp | 116 +++ src/meshset/CMakeLists.txt | 2 +- src/meshset/include/exaDEM/face.h | 6 +- src/polyhedra/CMakeLists.txt | 2 +- src/polyhedra/compute_hooke_interaction.cpp | 6 + .../exaDEM/compute_hooke_interaction.h | 177 +++- src/polyhedra/include/exaDEM/shape/shape.hpp | 17 +- .../include/exaDEM/shape/shape_prepro.hpp | 99 +- .../include/exaDEM/stl_mesh_to_driver.h | 110 ++ src/polyhedra/read_stl_mesh_to_shp.cpp | 67 ++ src/polyhedra/update_grid_interaction.cpp | 954 ++++++++++-------- 29 files changed, 1712 insertions(+), 512 deletions(-) create mode 100644 example/polyhedra/stl_mesh/CMakeLists.txt create mode 100644 example/polyhedra/stl_mesh/alpha3.shp create mode 100644 example/polyhedra/stl_mesh/box.stl create mode 100644 example/polyhedra/stl_mesh/box_for_octa.stl create mode 100644 example/polyhedra/stl_mesh/octahedron.shp create mode 100644 example/polyhedra/stl_mesh/stl_mesh_box_hexapod.msp create mode 100644 example/polyhedra/stl_mesh/stl_mesh_box_hexapod_mini.msp create mode 100644 example/polyhedra/stl_mesh/stl_mesh_box_octahedron.msp create mode 100644 example/polyhedra/stl_mesh/stl_mesh_box_octahedron_mini.msp create mode 100644 src/driver/add_stl_mesh.cpp create mode 100644 src/driver/include/exaDEM/driver_stl_mesh.h create mode 100644 src/driver/update_grid_stl_mesh.cpp create mode 100644 src/polyhedra/include/exaDEM/stl_mesh_to_driver.h create mode 100644 src/polyhedra/read_stl_mesh_to_shp.cpp diff --git a/data/config/config_move_particles.msp b/data/config/config_move_particles.msp index 417a810..4f4004e 100644 --- a/data/config/config_move_particles.msp +++ b/data/config/config_move_particles.msp @@ -50,18 +50,13 @@ ghost_full_update: - compute_new_vertices -################## trigger stl mesh update ######################## - -trigger_update_stl_mesh: - condition: enable_stl_mesh - body: - - build_grid_stl_mesh ################### parallel particle migration ############################ ## These operator should be define by your particle mode (spheres or polyhedra) migrate_particles: nop add_generated_particles: nop +update_stl_mesh: nop parallel_update_particles: - particle_generator: @@ -72,8 +67,8 @@ parallel_update_particles: - migrate_particles - ghost_full_update - grid_post_processing + - update_stl_mesh - update_particle_neighbors - - trigger_update_stl_mesh # define actions to initialize particles at startup, just after file read init_particles: diff --git a/data/config/config_polyhedra.msp b/data/config/config_polyhedra.msp index 1c78645..d2ea74d 100644 --- a/data/config/config_polyhedra.msp +++ b/data/config/config_polyhedra.msp @@ -1,4 +1,16 @@ ######### POLYHEDRA mode #################### + +######## STL Mesh ########################## + +trigger_update_stl_mesh: + condition: enable_stl_mesh AND trigger_load_balance + body: + - update_grid_stl_mesh + +update_stl_mesh: trigger_update_stl_mesh + +######### Move Particles ################### + trigger_move_particles: rebind: { threshold: max_displ , result: trigger_move_particles } body: @@ -8,6 +20,8 @@ trigger_move_particles: # - write_paraview_polyhedra: # basename: polyhedra + + +chunk_neighbors_impl: - update_grid_interaction - update_mutexes diff --git a/data/config/config_spheres.msp b/data/config/config_spheres.msp index 893ed04..b1086a7 100644 --- a/data/config/config_spheres.msp +++ b/data/config/config_spheres.msp @@ -1,4 +1,15 @@ ######### SPHERES mode #################### + +################## trigger stl mesh update ######################## + +trigger_update_stl_mesh: + condition: enable_stl_mesh AND trigger_load_balance + body: + - build_grid_stl_mesh + +update_stl_mesh: trigger_update_stl_mesh + +################# Move spheres #################################### migrate_particles: - migrate_cell_particles_friction - rebuild_amr diff --git a/example/polyhedra/CMakeLists.txt b/example/polyhedra/CMakeLists.txt index 23562f5..abcf104 100644 --- a/example/polyhedra/CMakeLists.txt +++ b/example/polyhedra/CMakeLists.txt @@ -1,6 +1,6 @@ enable_testing() -list(APPEND DLIST generator rotating_drum balls) +list(APPEND DLIST generator rotating_drum balls stl_mesh) foreach(DirExample IN LISTS DLIST) add_subdirectory(${DirExample}) diff --git a/example/polyhedra/stl_mesh/CMakeLists.txt b/example/polyhedra/stl_mesh/CMakeLists.txt new file mode 100644 index 0000000..46ab16a --- /dev/null +++ b/example/polyhedra/stl_mesh/CMakeLists.txt @@ -0,0 +1,19 @@ + +set(ExampleName "STLMeshHexapods") +set(FileName "stl_mesh_box_hexapod_mini.msp") + +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/alpha3.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/octahedron.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/box.stl DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/box_for_octa.stl DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + +add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) +add_test(Example_ExaDEM_${ExampleName}_omp ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 2) +add_test(Example_ExaDEM_${ExampleName}_mpi mpirun -n 2 ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) + +set(ExampleName "STLMeshOctahedra") +set(FileName "stl_mesh_box_octahedron_mini.msp") + +add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) +add_test(Example_ExaDEM_${ExampleName}_omp ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 2) +add_test(Example_ExaDEM_${ExampleName}_mpi mpirun -n 2 ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) diff --git a/example/polyhedra/stl_mesh/alpha3.shp b/example/polyhedra/stl_mesh/alpha3.shp new file mode 100644 index 0000000..54f1e1f --- /dev/null +++ b/example/polyhedra/stl_mesh/alpha3.shp @@ -0,0 +1,25 @@ + +< +name alpha3 +radius 0.237770037899799 +preCompDone y +nv 6 +0.475540075799599 0 0 +-0.475540075799599 0 0 +0 0.475540075799599 0 +0 -0.475540075799599 0 +0 0 0.475540075799599 +0 0 -0.475540075799599 +ne 3 +0 1 +2 3 +4 5 +nf 0 +obb.extent 0.713310113699398 0.713310113699398 0.713310113699398 +obb.e1 1 0 0 +obb.e2 0 1 0 +obb.e3 0 0 1 +obb.center 0 0 0 +volume 0.523598775598299 +I/m 0.13266 0.13266 0.13266 +> diff --git a/example/polyhedra/stl_mesh/box.stl b/example/polyhedra/stl_mesh/box.stl new file mode 100644 index 0000000..2d3f23c --- /dev/null +++ b/example/polyhedra/stl_mesh/box.stl @@ -0,0 +1,73 @@ +solid Created by Gmsh +facet normal -1 0 0 + outer loop + vertex -1 -1 19 + vertex -1 19 -1 + vertex -1 -1 -1 + endloop +endfacet +facet normal -1 0 0 + outer loop + vertex -1 -1 19 + vertex -1 19 19 + vertex -1 19 -1 + endloop +endfacet +facet normal 1 0 0 + outer loop + vertex 19 -1 19 + vertex 19 -1 -1 + vertex 19 19 -1 + endloop +endfacet +facet normal 1 -0 0 + outer loop + vertex 19 -1 19 + vertex 19 19 -1 + vertex 19 19 19 + endloop +endfacet +facet normal 0 -1 -0 + outer loop + vertex 19 -1 19 + vertex -1 -1 -1 + vertex 19 -1 -1 + endloop +endfacet +facet normal -0 -1 0 + outer loop + vertex 19 -1 19 + vertex -1 -1 19 + vertex -1 -1 -1 + endloop +endfacet +facet normal 0 1 0 + outer loop + vertex 19 19 19 + vertex 19 19 -1 + vertex -1 19 -1 + endloop +endfacet +facet normal 0 1 0 + outer loop + vertex 19 19 19 + vertex -1 19 -1 + vertex -1 19 19 + endloop +endfacet +facet normal -0 0 -1 + outer loop + vertex 19 19 -1 + vertex -1 -1 -1 + vertex -1 19 -1 + endloop +endfacet +facet normal 0 -0 -1 + outer loop + vertex 19 19 -1 + vertex 19 -1 -1 + vertex -1 -1 -1 + endloop +endfacet + +endsolid Created by Gmsh diff --git a/example/polyhedra/stl_mesh/box_for_octa.stl b/example/polyhedra/stl_mesh/box_for_octa.stl new file mode 100644 index 0000000..a4e0193 --- /dev/null +++ b/example/polyhedra/stl_mesh/box_for_octa.stl @@ -0,0 +1,73 @@ +solid Created by Gmsh +facet normal -1 0 0 + outer loop + vertex -1 -1 9 + vertex -1 9 -1 + vertex -1 -1 -1 + endloop +endfacet +facet normal -1 0 0 + outer loop + vertex -1 -1 9 + vertex -1 9 9 + vertex -1 9 -1 + endloop +endfacet +facet normal 1 0 0 + outer loop + vertex 9 -1 9 + vertex 9 -1 -1 + vertex 9 9 -1 + endloop +endfacet +facet normal 1 -0 0 + outer loop + vertex 9 -1 9 + vertex 9 9 -1 + vertex 9 9 9 + endloop +endfacet +facet normal 0 -1 -0 + outer loop + vertex 9 -1 9 + vertex -1 -1 -1 + vertex 9 -1 -1 + endloop +endfacet +facet normal -0 -1 0 + outer loop + vertex 9 -1 9 + vertex -1 -1 9 + vertex -1 -1 -1 + endloop +endfacet +facet normal 0 1 0 + outer loop + vertex 9 9 9 + vertex 9 9 -1 + vertex -1 9 -1 + endloop +endfacet +facet normal 0 1 0 + outer loop + vertex 9 9 9 + vertex -1 9 -1 + vertex -1 9 9 + endloop +endfacet +facet normal -0 0 -1 + outer loop + vertex 9 9 -1 + vertex -1 -1 -1 + vertex -1 9 -1 + endloop +endfacet +facet normal 0 -0 -1 + outer loop + vertex 9 9 -1 + vertex 9 -1 -1 + vertex -1 -1 -1 + endloop +endfacet + +endsolid Created by Gmsh diff --git a/example/polyhedra/stl_mesh/octahedron.shp b/example/polyhedra/stl_mesh/octahedron.shp new file mode 100644 index 0000000..82acda4 --- /dev/null +++ b/example/polyhedra/stl_mesh/octahedron.shp @@ -0,0 +1,44 @@ + +< +name Octahedron +radius 0.1 +preCompDone y +nv 6 +0.2310789034541148 -0.2310789034541148 0.0 +0.2310789034541148 0.2310789034541148 0.0 +0.0 0.0 0.32679491924311227 +-0.2310789034541148 -0.2310789034541148 0.0 +-0.2310789034541148 0.2310789034541148 0.0 +0.0 0.0 -0.32679491924311227 +ne 12 +0 1 +2 1 +2 0 +0 3 +2 3 +3 4 +4 2 +4 1 +5 0 +5 1 +5 4 +5 3 +nf 8 +3 0 1 2 +3 2 3 4 +3 1 2 4 +3 0 2 3 +3 0 5 1 +3 0 5 3 +3 3 5 4 +3 4 5 1 +obb.extent 0.33107890345411484 0.33107890345411484 0.4267949192431123 +obb.e1 1.0 0.0 0.0 +obb.e2 0.0 1.0 0.0 +obb.e3 0.0 0.0 1.0 +obb.center 0.0 0.0 0.0 +position 0.0 0.0 0.0 +orientation 1.0 0.0 0.0 0.0 +volume 0.16666666666666666 +I/m 0.04999999999999999 0.04999999999999999 0.04999999999999999 +> diff --git a/example/polyhedra/stl_mesh/stl_mesh_box_hexapod.msp b/example/polyhedra/stl_mesh/stl_mesh_box_hexapod.msp new file mode 100644 index 0000000..57a9906 --- /dev/null +++ b/example/polyhedra/stl_mesh/stl_mesh_box_hexapod.msp @@ -0,0 +1,57 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: SC + types: [ 0 ] + size: [ 1.5 , 1.5 , 1.5 ] + repeats: [ 12 , 12 , 50 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: alpha3.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_stl_mesh: + id: 0 + filename: box.stl + minskowski: 0.01 + - update_grid_stl_mesh + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,false,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 50000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: 1000 + dt: 0.0001 s + rcut_inc: 0.1 m + enable_stl_mesh: true diff --git a/example/polyhedra/stl_mesh/stl_mesh_box_hexapod_mini.msp b/example/polyhedra/stl_mesh/stl_mesh_box_hexapod_mini.msp new file mode 100644 index 0000000..8adcbd3 --- /dev/null +++ b/example/polyhedra/stl_mesh/stl_mesh_box_hexapod_mini.msp @@ -0,0 +1,56 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: SC + types: [ 0 ] + size: [ 1.5 , 1.5 , 1.5 ] + repeats: [ 3 , 3 , 3 ] + - read_shape_file: + filename: alpha3.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_stl_mesh: + id: 0 + filename: box.stl + minskowski: 0.01 + - update_grid_stl_mesh + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,false,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 20000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: -1 + dt: 0.0001 s + rcut_inc: 0.1 m + enable_stl_mesh: true diff --git a/example/polyhedra/stl_mesh/stl_mesh_box_octahedron.msp b/example/polyhedra/stl_mesh/stl_mesh_box_octahedron.msp new file mode 100644 index 0000000..657e3b2 --- /dev/null +++ b/example/polyhedra/stl_mesh/stl_mesh_box_octahedron.msp @@ -0,0 +1,57 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: SC + types: [ 0 ] + size: [ 0.8 , 0.8 , 1.0 ] + repeats: [ 11 , 11 , 30 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: octahedron.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_stl_mesh: + id: 0 + filename: box_for_octa.stl + minskowski: 0.01 + - update_grid_stl_mesh + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,false,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 50000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: 1000 + dt: 0.0001 s + rcut_inc: 0.1 m + enable_stl_mesh: true diff --git a/example/polyhedra/stl_mesh/stl_mesh_box_octahedron_mini.msp b/example/polyhedra/stl_mesh/stl_mesh_box_octahedron_mini.msp new file mode 100644 index 0000000..49e7b1d --- /dev/null +++ b/example/polyhedra/stl_mesh/stl_mesh_box_octahedron_mini.msp @@ -0,0 +1,57 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: SC + types: [ 0 ] + size: [ 0.8 , 0.8 , 1.0 ] + repeats: [ 3 , 3 , 3 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: octahedron.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_stl_mesh: + id: 0 + filename: box_for_octa.stl + minskowski: 0.01 + - update_grid_stl_mesh + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,false,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 20000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: -1 + dt: 0.0001 s + rcut_inc: 0.1 m + enable_stl_mesh: true diff --git a/src/driver/CMakeLists.txt b/src/driver/CMakeLists.txt index 0273d21..d005e85 100644 --- a/src/driver/CMakeLists.txt +++ b/src/driver/CMakeLists.txt @@ -1,2 +1,2 @@ -set(exadem_driver_LINK_LIBRARIES exanbIO exanbDefBox exanbParticleNeighbors exadem_numerical_scheme exadem_friction exadem_force_field exadem_io) +set(exadem_driver_LINK_LIBRARIES exanbIO exanbDefBox exanbParticleNeighbors exadem_numerical_scheme exadem_friction exadem_force_field exadem_io exadem_mesh_set exadem_polyhedra) xstamp_add_plugin(exadem_driver ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/driver/add_stl_mesh.cpp b/src/driver/add_stl_mesh.cpp new file mode 100644 index 0000000..0a2664b --- /dev/null +++ b/src/driver/add_stl_mesh.cpp @@ -0,0 +1,75 @@ +//#pragma xstamp_cuda_enable //! DO NOT REMOVE THIS LINE +#include "exanb/core/operator.h" +#include "exanb/core/operator_slot.h" +#include "exanb/core/operator_factory.h" +#include "exanb/core/make_grid_variant_operator.h" +#include "exanb/core/parallel_grid_algorithm.h" +#include "exanb/core/grid.h" +#include "exanb/core/domain.h" +#include "exanb/compute/compute_cell_particles.h" +#include +#include +#include +#include +#include +#include +#include +#include + +namespace exaDEM +{ + + using namespace exanb; + + template + class AddSTLMesh : public OperatorNode + { + static constexpr Vec3d null= { 0.0, 0.0, 0.0 }; + + ADD_SLOT( Drivers , drivers , INPUT_OUTPUT, DocString{"List of Drivers"}); + ADD_SLOT( int , id , INPUT , REQUIRED , DocString{"Driver index"}); + ADD_SLOT( std::string , filename , INPUT , REQUIRED , DocString{"Input filename"}); + ADD_SLOT( Vec3d , center , INPUT , null , DocString{"Defined but not used"}); + ADD_SLOT( Vec3d , angular_velocity, INPUT , null , DocString{"Defined but not used"}); + ADD_SLOT( Vec3d , velocity , INPUT , null , DocString{"Defined but not used"}); + ADD_SLOT( double , minskowski , INPUT , REQUIRED , DocString{"Minskowski radius value"} ); + ADD_SLOT( double , rcut_inc , INPUT , DocString{"value added to the search distance to update neighbor list less frequently. in physical space"} ); + + public: + + inline std::string documentation() const override final + { + return R"EOF( + This operator add a stl mesh to the drivers list. + )EOF"; + } + + inline void execute () override final + { + stl_mesh mesh; + mesh.read_stl(*filename); + std::string output_name_vtk = *filename; + std::string old_extension = ".stl"; + + std::string::size_type it = output_name_vtk.find(old_extension); + if(it != std::string::npos) + { + output_name_vtk.erase(it, old_extension.length()); + } + shape shp = build_shape(mesh, output_name_vtk); + shp.m_radius = *minskowski; + shp.increase_obb (*rcut_inc); + exaDEM::Stl_mesh driver= {*center, *velocity, *angular_velocity, shp}; + drivers->add_driver(*id, driver); + } + }; + + template using AddSTLMeshTmpl = AddSTLMesh; + + // === register factories === + CONSTRUCTOR_FUNCTION + { + OperatorNodeFactory::instance()->register_factory( "add_stl_mesh", make_grid_variant_operator< AddSTLMeshTmpl > ); + } +} + diff --git a/src/driver/include/exaDEM/driver_base.h b/src/driver/include/exaDEM/driver_base.h index cfe1af4..6495708 100644 --- a/src/driver/include/exaDEM/driver_base.h +++ b/src/driver/include/exaDEM/driver_base.h @@ -11,6 +11,7 @@ namespace exaDEM CYLINDER, SURFACE, BALL, + STL_MESH, UNDEFINED }; @@ -23,6 +24,7 @@ namespace exaDEM case DRIVER_TYPE::CYLINDER: return "Cylinder"; case DRIVER_TYPE::SURFACE: return "Surface"; case DRIVER_TYPE::BALL: return "Ball"; + case DRIVER_TYPE::STL_MESH: return "Stl_mesh"; case DRIVER_TYPE::UNDEFINED: return "Undefined Driver"; default: return "Undefined Driver"; } @@ -40,6 +42,7 @@ namespace exaDEM case str2int("CYLINDER"): return DRIVER_TYPE::CYLINDER; case str2int("SURFACE"): return DRIVER_TYPE::SURFACE; case str2int("BALL"): return DRIVER_TYPE::BALL; + case str2int("STL_MESH"): return DRIVER_TYPE::STL_MESH; default: std::cout << "error, no driver " << driver_name << " found" << std::endl; std::cout << "Use: CYLINDER, SURFACE, or BALL" << std::endl; std::abort(); @@ -49,6 +52,7 @@ namespace exaDEM struct Cylinder; struct Surface; struct Ball; + struct Stl_mesh; struct UndefinedDriver; template @@ -56,6 +60,7 @@ namespace exaDEM template<> constexpr DRIVER_TYPE get_type () { return DRIVER_TYPE::CYLINDER; } template<> constexpr DRIVER_TYPE get_type () { return DRIVER_TYPE::SURFACE; } template<> constexpr DRIVER_TYPE get_type () { return DRIVER_TYPE::BALL; } + template<> constexpr DRIVER_TYPE get_type () { return DRIVER_TYPE::STL_MESH; } template<> constexpr DRIVER_TYPE get_type () { return DRIVER_TYPE::UNDEFINED; } struct Driver diff --git a/src/driver/include/exaDEM/driver_stl_mesh.h b/src/driver/include/exaDEM/driver_stl_mesh.h new file mode 100644 index 0000000..ddf7c2e --- /dev/null +++ b/src/driver/include/exaDEM/driver_stl_mesh.h @@ -0,0 +1,78 @@ +#pragma once +#include +#include +#include + +namespace exaDEM +{ + using namespace exanb; + + template using vector_t = std::vector; + + struct list_of_elements + { + vector_t vertices; + vector_t edges; + vector_t faces; + }; + + struct Stl_mesh + { + exanb::Vec3d center; // normal * offset + exanb::Vec3d vel; // 0,0,0 + exanb::Vec3d vrot; + shape shp; + vector_t grid_indexes; + + + constexpr DRIVER_TYPE get_type() {return DRIVER_TYPE::STL_MESH;} + + void print() + { + std::cout << "Driver Type: MESH STL" << std::endl; + std::cout << "center : " << center << std::endl; + std::cout << "Vel : " << vel << std::endl; + std::cout << "AngVel : " << vrot << std::endl; + std::cout << "Number of faces : " << shp.get_number_of_faces() << std::endl; + std::cout << "Number of edges : " << shp.get_number_of_edges() << std::endl; + std::cout << "Number of vertices : " << shp.get_number_of_vertices() << std::endl; + } + + inline void initialize () + { + // checks + } + + inline void update_position ( const double t ) + { + center = center + t * vel; + } + + inline void grid_indexes_summary() + { + const size_t size = grid_indexes.size(); + size_t nb_fill_cells(0), nb_v (0), nb_e(0), nb_f(0), max_v(0), max_e(0), max_f(0); + +#pragma omp parallel for reduction(+: nb_fill_cells, nb_v, nb_e, nb_f) reduction(max: max_v, max_e, max_f) + for(size_t i = 0 ; i < size ; i++) + { + auto& list = grid_indexes[i]; + if ( list.vertices.size() == 0 && list.edges.size() == 0 && list.faces.size()) continue; + nb_fill_cells++; + nb_v += list.vertices.size(); + nb_e += list.edges.size(); + nb_f += list.faces.size(); + max_v = std::max( max_v, list.vertices.size() ); + max_e = std::max( max_e, list.edges.size() ); + max_f = std::max( max_f, list.faces.size() ); + } + + lout << "========= STL Grid summary ======" << std::endl; + lout << "Number of emplty cells = " << nb_fill_cells << " / " << size << std::endl; + lout << "Vertices (Total/Max) = " << nb_v << " / " << max_v << std::endl; + lout << "Edges (Total/Max) = " << nb_e << " / " << max_e << std::endl; + lout << "Faces (Total/Max) = " << nb_f << " / " << max_f << std::endl; + lout << "=================================" << std::endl; + } + }; +} diff --git a/src/driver/include/exaDEM/drivers.h b/src/driver/include/exaDEM/drivers.h index 8288ab3..10b8ab0 100644 --- a/src/driver/include/exaDEM/drivers.h +++ b/src/driver/include/exaDEM/drivers.h @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -16,7 +17,7 @@ namespace exaDEM struct Drivers { template using vector_t = onika::memory::CudaMMVector; - using data_t = std::variant; + using data_t = std::variant; using driver_t = DRIVER_TYPE; vector_t m_type; vector_t m_data; @@ -64,27 +65,27 @@ namespace exaDEM // Accessors ONIKA_HOST_DEVICE_FUNC - inline DRIVER_TYPE type(int idx) + inline DRIVER_TYPE type(size_t idx) { assert( idx < m_type.size()); return m_type[idx]; } ONIKA_HOST_DEVICE_FUNC - inline const data_t& data(int idx) const + inline const data_t& data(size_t idx) const { assert( idx < m_data.size()); return m_data[idx]; } ONIKA_HOST_DEVICE_FUNC - inline data_t& data(int idx) + inline data_t& data(size_t idx) { assert( idx < m_data.size()); return m_data[idx]; } - bool well_former() + bool well_defined() { for( auto& it : m_type ) { diff --git a/src/driver/update_grid_stl_mesh.cpp b/src/driver/update_grid_stl_mesh.cpp new file mode 100644 index 0000000..2a40264 --- /dev/null +++ b/src/driver/update_grid_stl_mesh.cpp @@ -0,0 +1,116 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for ONIKA_ASSUME_ALIGNED macro +#include +#include +#include + +#include +#include +#include +#include + +#include + +namespace exaDEM +{ + using namespace exanb; + template< class GridT, class = AssertGridHasFields< GridT >> class UpdateGridSTLMeshOperator : public OperatorNode + { + using ComputeFields = FieldSet< field::_rx ,field::_ry ,field::_rz>; + static constexpr ComputeFields compute_field_set {}; + ADD_SLOT( MPI_Comm , mpi , INPUT , MPI_COMM_WORLD , DocString{"MPI communicator for parallel processing."}); + ADD_SLOT( GridT , grid , INPUT_OUTPUT , DocString{"Grid used for computations."} ); + ADD_SLOT( Drivers , drivers , INPUT_OUTPUT, DocString{"List of Drivers"}); + + public: + inline std::string documentation() const override final + { + return R"EOF( + )EOF"; + } + + + // could be speedup by an OBB TREE + void update_indexes ( list_of_elements& list, OBB& bx, shape& shp) + { + // first vertices + auto& obb_v = shp.m_obb_vertices; + auto& obb_e = shp.m_obb_edges; + auto& obb_f = shp.m_obb_faces; + list.vertices.clear(); + list.faces.clear(); + list.edges.clear(); + for(size_t it = 0; it < obb_v.size() ; it++) + { + if ( bx.intersect( obb_v[it] ) ) list.vertices.push_back( it ); + //list.vertices.push_back( it ); + } + + for(size_t it = 0; it < obb_e.size() ; it++) + { + if ( bx.intersect( obb_e[it] ) ) list.edges.push_back( it ); + //list.edges.push_back( it ); + } + + for(size_t it = 0; it < obb_f.size() ; it++) + { + if ( bx.intersect( obb_f[it] ) ) list.faces.push_back( it ); + //list.faces.push_back( it ); + } + } + + inline void execute () override final + { + const size_t n_cells = grid->number_of_cells(); // nbh.size(); + const IJK dims = grid->dimension(); + const int gl = grid->ghost_layers(); + const double csize = grid->cell_size(); + + for(size_t id = 0 ; id < drivers->get_size() ; id++) + { + if ( drivers->type(id) == DRIVER_TYPE::STL_MESH) + { + exaDEM::Stl_mesh& mesh = std::get (drivers->data(id)); + auto& grid_stl = mesh.grid_indexes; + grid_stl.clear(); + grid_stl.resize(n_cells); + +# pragma omp parallel + { + OBB bx; + bx.extent = {0.5 * csize, 0.5 * csize, 0.5 * csize}; + GRID_OMP_FOR_BEGIN(dims-2*gl,_,block_loc, schedule(dynamic)) + { + IJK loc_a = block_loc + gl; + size_t cell_a = grid_ijk_to_index( dims , loc_a ); + auto cb = grid->cell_bounds(loc_a); + auto center = (cb.bmin + cb.bmax) / 2; + bx.center = { center.x , center.y, center.z}; + update_indexes ( grid_stl[cell_a], bx, mesh.shp); + } + GRID_OMP_FOR_END + } + //mesh.grid_indexes_summary(); for debug + } + } + }; + }; + + // this helps older versions of gcc handle the unnamed default second template parameter + template using UpdateGridSTLMeshOperatorTemplate = UpdateGridSTLMeshOperator; + + // === register factories === + CONSTRUCTOR_FUNCTION + { + OperatorNodeFactory::instance()->register_factory( "update_grid_stl_mesh", make_grid_variant_operator< UpdateGridSTLMeshOperatorTemplate > ); + } +} diff --git a/src/meshset/CMakeLists.txt b/src/meshset/CMakeLists.txt index becc263..bc2948e 100644 --- a/src/meshset/CMakeLists.txt +++ b/src/meshset/CMakeLists.txt @@ -1,2 +1,2 @@ -set(exadem_mesh_set_LINK_LIBRARIES exanbDefBox exanbCompute exanbIO exanbMPI exanbParticleNeighbors exanbGridCellParticles exadem_friction exadem_force_field) +set(exadem_mesh_set_LINK_LIBRARIES exanbDefBox exanbCompute exanbIO exanbMPI exanbParticleNeighbors exanbGridCellParticles exadem_friction exadem_force_field exadem_driver) xstamp_add_plugin(exadem_mesh_set ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/meshset/include/exaDEM/face.h b/src/meshset/include/exaDEM/face.h index 14b3735..cfb5a8c 100644 --- a/src/meshset/include/exaDEM/face.h +++ b/src/meshset/include/exaDEM/face.h @@ -44,7 +44,7 @@ namespace exaDEM * @brief Normalize a 3D vector. * @param v The input vector to be normalized. */ - inline void normalize(Vec3d& v) { v = v / exanb::norm(v); } + inline void _normalize(Vec3d& v) { v = v / exanb::norm(v); } /** @@ -106,9 +106,9 @@ namespace exaDEM const Vec3d& pc = vertices[nb_vertices-1]; Vec3d v1 = pb - pa; Vec3d v2 = pc - pa; - normalize(v1); + _normalize(v1); Vec3d n = exanb::cross(v1,v2); - normalize(n); + _normalize(n); Vec3d iv = center;// - pa; double dist = exanb::dot(iv,n); if(dist < 0.0) diff --git a/src/polyhedra/CMakeLists.txt b/src/polyhedra/CMakeLists.txt index 0a78605..dfcd228 100644 --- a/src/polyhedra/CMakeLists.txt +++ b/src/polyhedra/CMakeLists.txt @@ -1,2 +1,2 @@ -set(exadem_polyhedra_LINK_LIBRARIES exanbDefBox exanbCompute exanbIO exanbMPI exanbParticleNeighbors exanbGridCellParticles exanbExtraStorage exadem_driver exadem_force_field exadem_io) +set(exadem_polyhedra_LINK_LIBRARIES exanbDefBox exanbCompute exanbIO exanbMPI exanbParticleNeighbors exanbGridCellParticles exanbExtraStorage exadem_driver exadem_force_field exadem_io exadem_mesh_set) xstamp_add_plugin(exadem_polyhedra ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/src/polyhedra/compute_hooke_interaction.cpp b/src/polyhedra/compute_hooke_interaction.cpp index 0ea349d..40c4b4f 100644 --- a/src/polyhedra/compute_hooke_interaction.cpp +++ b/src/polyhedra/compute_hooke_interaction.cpp @@ -20,6 +20,7 @@ #include #include #include +#include namespace exaDEM { @@ -261,6 +262,11 @@ size_t current_cell = indexes[ci]; Ball& driver = std::get(drvs.data(driver_idx)) ; compute_driver(cells[item.cell_i], driver, hkp_drvs, shp, item, time, locker); } + else if ( item.type >= 7 && item.type <= 12) + { + Stl_mesh& driver = std::get(drvs.data(driver_idx)) ; + compute_driver_stl(cells[item.cell_i], driver, hkp_drvs, shp, item, time, locker); + } } } } diff --git a/src/polyhedra/include/exaDEM/compute_hooke_interaction.h b/src/polyhedra/include/exaDEM/compute_hooke_interaction.h index 7df9b95..a451243 100644 --- a/src/polyhedra/include/exaDEM/compute_hooke_interaction.h +++ b/src/polyhedra/include/exaDEM/compute_hooke_interaction.h @@ -1,47 +1,168 @@ #pragma once -#include +#include +#include +#include namespace exaDEM { - template - struct HookeForceInteractionFunctor - { - GridT& cells; - mutexes& locker; - const shapes & shps; - const HookeParams params; - const double time; - inline const Vec3d get_r(const int cell_id, const int p_id) + // C for cell and D for driver + template + ONIKA_HOST_DEVICE_FUNC void compute_driver(C& cell, D& driver, const HookeParams& hkp, const shape* shp, Interaction& I, const double time, mutexes& locker) { - const Vec3d res = { - cells[cell_id][field::rx][p_id], - cells[cell_id][field::ry][p_id], - cells[cell_id][field::rz][p_id]}; - return res; - }; + const size_t p = I.p_i; + const size_t sub = I.sub_i; + // === positions + const Vec3d r = { cell[field::rx][p], cell[field::ry][p], cell[field::rz][p] }; + // === vrot + const Vec3d& vrot = cell[field::vrot][p]; + // === vertex array + const auto& vertices = cell[field::vertices][p]; - inline const Vec3d get_v(const int cell_id, const int p_id) - { - const Vec3d res = { - cells[cell_id][field::vx][p_id], - cells[cell_id][field::vy][p_id], - cells[cell_id][field::vz][p_id]}; - return res; - }; + auto [contact, dn, n, contact_position] = exaDEM::detector_vertex_driver(driver, vertices, sub, shp); + + if(contact) + { + constexpr Vec3d null = {0,0,0}; + auto& mom = cell[field::mom][p]; + const Vec3d v = { cell[field::vx][p], cell[field::vy][p], cell[field::vz][p] }; + const double meff = cell[field::mass][p]; + Vec3d f = null; + hooke_force_core(dn, n, time, hkp.m_kn, hkp.m_kt, hkp.m_kr, + hkp.m_mu, hkp.m_damp_rate, meff, + I.friction, contact_position, + r, v, f, I.moment, vrot, // particle i + driver.center, driver.vel, driver.vrot // particle j + ); - ONIKA_HOST_DEVICE_FUNC inline void compute_force_interaction ( ) + // === update informations + locker.lock(I.cell_i, p); + mom += compute_moments(contact_position, r, f, I.moment); + cell[field::fx][p] += f.x; + cell[field::fy][p] += f.y; + cell[field::fz][p] += f.z; + locker.unlock(I.cell_i, p); + } + else { + I.reset(); } + } - ONIKA_HOST_DEVICE_FUNC inline void operator () ( Interaction * __restrict__ ptr, const size_t size) const + struct stl_mesh_dectector + { + ONIKA_HOST_DEVICE_FUNC std::tuple operator() ( + const uint16_t type, + const Vec3d& pi, const int i, const shape* shpi, const exanb::Quaternion& oi, + const Vec3d& pj, const int j, const shape* shpj, const exanb::Quaternion& oj) const { + #define __params__ pi, i, shpi, oi, pj, j, shpj, oj + #define __inv_params__ pj, j, shpj, oj, pi, i, shpi, oi + assert( type >= 7 && type <= 12 ); + switch (type) + { + case 7: return exaDEM::detection_vertex_vertex ( __params__ ); + case 8: return exaDEM::detection_vertex_edge ( __params__ ); + case 9: return exaDEM::detection_vertex_face ( __params__ ); + case 10: return exaDEM::detection_edge_edge ( __params__ ); + case 11: return exaDEM::detection_vertex_edge ( __params__ ); + case 12: return exaDEM::detection_vertex_face ( __inv_params__ ); + } + #undef __params__ + #undef __inv_params__ + return std::tuple(); } - ONIKA_HOST_DEVICE_FUNC inline void operator () ( Interaction& it) const + }; + + template + ONIKA_HOST_DEVICE_FUNC void compute_driver_stl(C& cell, Stl_mesh& driver, const HookeParams& hkp, const shape* shp_i, Interaction& I, const double time, mutexes& locker) { + const stl_mesh_dectector func; + const size_t p_i = I.p_i; + const size_t sub_i = I.sub_i; + const size_t sub_j = I.sub_j; + // === positions + const Vec3d r_i = { cell[field::rx][p_i], cell[field::ry][p_i], cell[field::rz][p_i] }; + // === vrot + const Vec3d& vrot_i = cell[field::vrot][p_i]; + const auto& t_i = cell[field::type][p_i]; + const Quaternion& orient_i = cell[field::orient][p_i]; + const auto& shp_j = driver.shp; + + // WARNING + const Quaternion orient_j = {1.0,0.0,0.0,0.0}; + auto [contact, dn, n, contact_position] = func(I.type, r_i, sub_i, &shp_i[t_i], orient_i, driver.center, sub_j, &shp_j, orient_j); + + if(contact) + { + constexpr Vec3d null = {0,0,0}; + auto& mom = cell[field::mom][p_i]; + const Vec3d v_i = { cell[field::vx][p_i], cell[field::vy][p_i], cell[field::vz][p_i] }; + const double meff = cell[field::mass][p_i]; + Vec3d f = null; + hooke_force_core(dn, n, time, hkp.m_kn, hkp.m_kt, hkp.m_kr, + hkp.m_mu, hkp.m_damp_rate, meff, + I.friction, contact_position, + r_i, v_i, f, I.moment, vrot_i, // particle i + driver.center, driver.vel, driver.vrot // particle j + ); + + // === update informations + locker.lock(I.cell_i, p_i); + mom += compute_moments(contact_position, r_i, f, I.moment); + cell[field::fx][p_i] += f.x; + cell[field::fy][p_i] += f.y; + cell[field::fz][p_i] += f.z; + locker.unlock(I.cell_i, p_i); + } + else + { + I.reset(); + } } - }; + + + + template + struct HookeForceInteractionFunctor + { + GridT& cells; + mutexes& locker; + const shapes & shps; + const HookeParams params; + const double time; + + inline const Vec3d get_r(const int cell_id, const int p_id) + { + const Vec3d res = { + cells[cell_id][field::rx][p_id], + cells[cell_id][field::ry][p_id], + cells[cell_id][field::rz][p_id]}; + return res; + }; + + inline const Vec3d get_v(const int cell_id, const int p_id) + { + const Vec3d res = { + cells[cell_id][field::vx][p_id], + cells[cell_id][field::vy][p_id], + cells[cell_id][field::vz][p_id]}; + return res; + }; + + ONIKA_HOST_DEVICE_FUNC inline void compute_force_interaction ( ) + { + } + + ONIKA_HOST_DEVICE_FUNC inline void operator () ( Interaction * __restrict__ ptr, const size_t size) const + { + } + + ONIKA_HOST_DEVICE_FUNC inline void operator () ( Interaction& it) const + { + } + }; } diff --git a/src/polyhedra/include/exaDEM/shape/shape.hpp b/src/polyhedra/include/exaDEM/shape/shape.hpp index dfc1961..19060de 100644 --- a/src/polyhedra/include/exaDEM/shape/shape.hpp +++ b/src/polyhedra/include/exaDEM/shape/shape.hpp @@ -34,6 +34,7 @@ namespace exaDEM VectorT m_vertices; ///< exanb::Vec3d m_inertia_on_mass; + VectorT m_obb_vertices; ///< only used for stl meshes VectorT m_obb_edges; VectorT m_obb_faces; OBB obb; @@ -44,8 +45,10 @@ namespace exaDEM std::string m_name = "undefined"; // functions in shape_prepro.hpp + inline void pre_compute_obb_vertices(); inline void pre_compute_obb_edges(); inline void pre_compute_obb_faces(); + inline void increase_obb ( const double value ); ONIKA_HOST_DEVICE_FUNC @@ -306,7 +309,7 @@ namespace exaDEM lout << "Vertex[" << idx++ << "]: [" << vertex.x << "," << vertex.y << "," << vertex.z << "]" << std::endl; }; - lout << "Number of vertices: " << this->get_number_of_vertices() << std::endl; + lout << "Number of vertices= " << this->get_number_of_vertices() << std::endl; for_all_vertices(printer); } @@ -323,7 +326,7 @@ namespace exaDEM } else { - lout << "Number of edges: " << this->get_number_of_edges() << std::endl; + lout << "Number of edge = " << this->get_number_of_edges() << std::endl; for_all_edges(printer); } } @@ -345,7 +348,7 @@ namespace exaDEM } else { - lout << "Number of faces: " << this->get_number_of_faces() << std::endl; + lout << "Number of faces = " << this->get_number_of_faces() << std::endl; for_all_faces(printer); } } @@ -353,10 +356,10 @@ namespace exaDEM inline void print() { lout << "======= Shape Configuration =====" << std::endl; - lout << "Shape Name: " << this->m_name << std::endl; - lout << "Shape Radius: " << this->m_radius << std::endl; - lout << "Shape I/m: [" << this->m_inertia_on_mass << "]" << std::endl; - lout << "Shape Volume: " << this->m_volume << std::endl; + lout << "Shape Name = " << this->m_name << std::endl; + lout << "Shape Radius = " << this->m_radius << std::endl; + lout << "Shape I/m = [" << this->m_inertia_on_mass << "]" << std::endl; + lout << "Shape Volume = " << this->m_volume << std::endl; print_vertices(); print_edges(); print_faces(); diff --git a/src/polyhedra/include/exaDEM/shape/shape_prepro.hpp b/src/polyhedra/include/exaDEM/shape/shape_prepro.hpp index 5f30e73..c0cf466 100644 --- a/src/polyhedra/include/exaDEM/shape/shape_prepro.hpp +++ b/src/polyhedra/include/exaDEM/shape/shape_prepro.hpp @@ -85,6 +85,15 @@ namespace exaDEM //#define OLD_VERSION // general functon; + inline OBB build_obb_vertex(const Vec3d& position, const int index, const shape* shp, const exanb::Quaternion& orientation) + { + const double ext = shp->m_radius; + auto vertex = shp->get_vertex(index, position, orientation); + std::vector v = {conv_to_vec3r(vertex)}; + OBB res = build_OBB (v, ext); + return res; + } + inline OBB build_obb_edge(const Vec3d& position, const int index, const shape* shp, const exanb::Quaternion& orientation) { const double ext = shp->m_radius; @@ -94,26 +103,6 @@ namespace exaDEM std::vector v = {conv_to_vec3r(vf), conv_to_vec3r(vs)}; OBB res = build_OBB (v, ext); return res; -#ifdef OLD_VERSION - // TODO use orientation - OBB res; - const Vec3d center = (vf + vs) / 2; - res.center = { center.x , center.y, center.z }; - - constexpr int DIM = 3; - // TODO: add obb orientation HERE - // Use distance between vf and vs and apply rotation - - const auto v = conv_to_vec3r (vf); - - for ( int dim = 0 ; dim < DIM ; dim++ ) - { - res.extent[dim] = std::abs ( res.center[dim] - v [dim] ); - } - - res.enlarge( ext ); - return res; -#endif } inline OBB build_obb_face(const Vec3d& position, const int index, const shape* shp, const exanb::Quaternion& orientation) @@ -127,39 +116,23 @@ namespace exaDEM } OBB res = build_OBB (v, ext); return res; -#ifdef OLD_VERSION - OBB obb_face; - assert ( nf > 0 ); + } - constexpr int DIM = 3; - for ( int dim = 0 ; dim < DIM ; dim++ ) - { - double acc = 0.0; - for (int f = 0 ; f < nf ; f++) - { - vec3r v = conv_to_vec3r (shp->get_vertex(data[f], position, orientation)); - acc += v[dim]; - } - obb_face.center[dim] = acc / nf; - } + inline void shape::pre_compute_obb_vertices() + { + // This function could be optimized by avoiding to use `position` and `orientation` in `build_obb_face` + const size_t size = this->get_number_of_vertices(); + m_obb_vertices.resize(size); + const exanb::Vec3d center = conv_to_Vec3d (this->obb.center); + const exanb::Quaternion onull = {1,0,0,0}; - // TODO: add obb orientation HERE - // Use distance between vf and vs and apply rotation - // TODO : we could do it with only one loop - for ( int dim = 0 ; dim < DIM ; dim++ ) + lout << "build " << size << " obb [vertices]" << std::endl; + +#pragma omp parallel for schedule (static) + for ( size_t i = 0 ; i < size ; i++) { - double acc = 0.0; - for (int f = 0 ; f < nf ; f++) - { - vec3r v = conv_to_vec3r (shp->get_vertex(data[f], position, orientation)); - acc = std::max ( acc , std::abs( v[dim] - obb_face.center[dim] )); - } - obb_face.extent[dim] = acc; + m_obb_vertices[i] = build_obb_vertex (center, i, this, onull); } - - obb_face.enlarge( ext ); - return obb_face; -#endif } inline void shape::pre_compute_obb_edges() @@ -171,7 +144,8 @@ namespace exaDEM //const exanb::Vec3d vnull = {0,0,0}; const exanb::Quaternion onull = {1,0,0,0}; const exanb::Vec3d center = conv_to_Vec3d (this->obb.center); -#pragma omp parallel for + lout << "build " << size << " obb [edges]" << std::endl; +#pragma omp parallel for schedule(static) for ( size_t i = 0 ; i < size ; i++) { m_obb_edges[i] = build_obb_edge (center, i, this, onull); @@ -186,11 +160,34 @@ namespace exaDEM //const exanb::Vec3d vnull = {0,0,0}; const exanb::Quaternion onull = {1,0,0,0}; const exanb::Vec3d center = conv_to_Vec3d (this->obb.center); + lout << "build " << size << " obb [faces]" << std::endl; -#pragma omp parallel for +#pragma omp parallel for schedule(static) for ( size_t i = 0 ; i < size ; i++) { m_obb_faces[i] = build_obb_face (center, i, this, onull); } } + + inline void shape::increase_obb ( const double value ) + { +#pragma omp parallel + { +#pragma omp for schedule(static) nowait + for (size_t i = 0 ; i < m_obb_vertices.size() ; i++) + { + m_obb_vertices[i].enlarge(value); + } +#pragma omp for schedule(static) nowait + for (size_t i = 0 ; i < m_obb_edges.size() ; i++) + { + m_obb_edges[i].enlarge(value); + } +#pragma omp for schedule(static) + for (size_t i = 0 ; i < m_obb_faces.size() ; i++) + { + m_obb_faces[i].enlarge(value); + } + } + } } diff --git a/src/polyhedra/include/exaDEM/stl_mesh_to_driver.h b/src/polyhedra/include/exaDEM/stl_mesh_to_driver.h new file mode 100644 index 0000000..1d214f2 --- /dev/null +++ b/src/polyhedra/include/exaDEM/stl_mesh_to_driver.h @@ -0,0 +1,110 @@ +#pragma once + +#include +#include +#include + +namespace exaDEM +{ + using namespace exanb; + shape build_shape(stl_mesh& mesh, std::string name) + { + shape shp; + const int n_faces = mesh.m_data.size(); + std::vector idx; + auto& shp_v = shp.m_vertices; + auto& shp_e = shp.m_edges; + auto& shp_f = shp.m_faces; + + // tmp edges + std::vector> tmp_e; + + shp_f.resize(1); + shp_f[0] = n_faces; + +#pragma omp parallel + { + // first get vertices + std::vector v; +#pragma omp for + for (int i = 0 ; i < n_faces ; i++ ) + { + Face& face = mesh.get_data( i ); + v.insert( std::end(v), std::begin( face.vertices ), std::end ( face.vertices ) ); + } + + auto last = std::unique(v.begin(), v.end()); + v.erase ( last, v.end() ); + std::sort(v.begin() , v.end()); + +#pragma omp critical + { + shp_v.insert( std::end(shp_v), std::begin(v), std::end(v)); + } + } + + { + auto last = std::unique(shp_v.begin(), shp_v.end()); + shp_v.erase ( last, shp_v.end() ); + std::sort(shp_v.begin() , shp_v.end()); + } + +#pragma omp parallel + { + // first get vertices + std::vector idxs; + std::vector> edges; + std::vector faces; +#pragma omp for + for (int i = 0 ; i < n_faces ; i++ ) + { + Face& face = mesh.get_data( i ); + auto& v = face.vertices; + idxs.resize ( v.size() ); + for( size_t idx = 0 ; idx < v.size() ; idx++ ) + { + auto it = std::lower_bound ( std::begin(shp_v), std::end(shp_v) , v[idx]); + idxs[idx] = std::distance ( shp_v.begin() , it ); + assert ( v[idx] == shp_v[idxs[idx]]); + } + // add edges and faces + faces.push_back( v.size() ); + for( size_t idx = 0 ; idx < v.size() ; idx++ ) + { + faces.push_back( idxs[idx] ); + size_t second ; + if( idx + 1 == v.size() ) second = idxs[0]; + else second = idxs[idx+1]; + edges.push_back( {idxs[idx], second} ); + } + } +#pragma omp critical + { + tmp_e.insert( std::end(tmp_e), std::begin(edges), std::end(edges)); + } +#pragma omp critical + { + shp_f.insert( std::end(shp_f), std::begin(faces), std::end(faces)); + } + } + + { + auto last = std::unique(tmp_e.begin(), tmp_e.end()); + tmp_e.erase ( last, tmp_e.end() ); + shp_e.resize( tmp_e.size() * 2 ); +#pragma omp parallel for schedule(static) + for(size_t i = 0 ; i< tmp_e.size() ; i++) + { + shp_e[2*i] = tmp_e[i].first; + shp_e[2*i+1] = tmp_e[i].second; + } + } + + shp.m_name = name; + shp.pre_compute_obb_vertices(); + shp.pre_compute_obb_edges(); + shp.pre_compute_obb_faces(); + shp.write_paraview(); + return shp; + } +} diff --git a/src/polyhedra/read_stl_mesh_to_shp.cpp b/src/polyhedra/read_stl_mesh_to_shp.cpp new file mode 100644 index 0000000..7f7bdcb --- /dev/null +++ b/src/polyhedra/read_stl_mesh_to_shp.cpp @@ -0,0 +1,67 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for ONIKA_ASSUME_ALIGNED macro +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +namespace exaDEM +{ + using namespace exanb; + template< class GridT, class = AssertGridHasFields< GridT >> class ReadSTLMeshOperator : public OperatorNode + { + using ComputeFields = FieldSet< field::_rx ,field::_ry ,field::_rz>; + static constexpr ComputeFields compute_field_set {}; + ADD_SLOT( MPI_Comm , mpi , INPUT , MPI_COMM_WORLD); + ADD_SLOT( GridT , grid , INPUT_OUTPUT ); + ADD_SLOT( Domain , domain , INPUT , REQUIRED ); + ADD_SLOT( std::string , filename, INPUT , REQUIRED , DocString{"Input filename"}); + + public: + inline std::string documentation() const override final + { + return R"EOF( This operator initialize a mesh composed of faces from an stl input file. + )EOF"; + } + + inline void execute () override final + { + stl_mesh mesh; + mesh.read_stl(*filename); + std::string output_name_vtk = *filename; + std::string old_extension = ".stl"; + + std::string::size_type it = output_name_vtk.find(old_extension); + if(it != std::string::npos) + { + output_name_vtk.erase(it, old_extension.length()); + } + shape shp = build_shape(mesh, output_name_vtk); + }; + }; + + // this helps older versions of gcc handle the unnamed default second template parameter + template using ReadSTLMeshOperatorTemplate = ReadSTLMeshOperator; + + // === register factories === + CONSTRUCTOR_FUNCTION + { + OperatorNodeFactory::instance()->register_factory( "read_stl_mesh", make_grid_variant_operator< ReadSTLMeshOperatorTemplate > ); + } +} diff --git a/src/polyhedra/update_grid_interaction.cpp b/src/polyhedra/update_grid_interaction.cpp index c4cfba2..039326c 100644 --- a/src/polyhedra/update_grid_interaction.cpp +++ b/src/polyhedra/update_grid_interaction.cpp @@ -20,429 +20,569 @@ namespace exaDEM { - using namespace exanb; - using VertexArray = ::onika::oarray_t<::exanb::Vec3d, 8>; + using namespace exanb; + using VertexArray = ::onika::oarray_t<::exanb::Vec3d, 8>; - template - > - class UpdateGridCellInteraction : public OperatorNode - { - using ComputeFields = FieldSet<>; - static constexpr ComputeFields compute_field_set {}; + template + > + class UpdateGridCellInteraction : public OperatorNode + { + using ComputeFields = FieldSet<>; + static constexpr ComputeFields compute_field_set {}; - ADD_SLOT( GridT , grid , INPUT_OUTPUT , REQUIRED ); - ADD_SLOT( exanb::GridChunkNeighbors , chunk_neighbors , INPUT , OPTIONAL , DocString{"Neighbor list"} ); - ADD_SLOT( GridCellParticleInteraction , ges , INPUT_OUTPUT , DocString{"Interaction list"} ); - ADD_SLOT( shapes , shapes_collection , INPUT , DocString{"Collection of shapes"}); - ADD_SLOT(double , rcut_inc , INPUT_OUTPUT , DocString{"value added to the search distance to update neighbor list less frequently. in physical space"} ); - ADD_SLOT( Drivers , drivers , INPUT , DocString{"List of Drivers"}); - ADD_SLOT( std::vector , idxs , INPUT_OUTPUT , DocString{"List of non empty cells"}); + ADD_SLOT( GridT , grid , INPUT_OUTPUT , REQUIRED ); + ADD_SLOT( exanb::GridChunkNeighbors , chunk_neighbors , INPUT , OPTIONAL , DocString{"Neighbor list"} ); + ADD_SLOT( GridCellParticleInteraction , ges , INPUT_OUTPUT , DocString{"Interaction list"} ); + ADD_SLOT( shapes , shapes_collection , INPUT , DocString{"Collection of shapes"}); + ADD_SLOT(double , rcut_inc , INPUT_OUTPUT , DocString{"value added to the search distance to update neighbor list less frequently. in physical space"} ); + ADD_SLOT( Drivers , drivers , INPUT , DocString{"List of Drivers"}); + ADD_SLOT( std::vector , idxs , INPUT_OUTPUT , DocString{"List of non empty cells"}); - public: + public: - inline std::string documentation() const override final - { - return R"EOF( + inline std::string documentation() const override final + { + return R"EOF( )EOF"; - } - - - template - void add_driver_interaction( D& driver, std::vector& driver_data, std::vector& driver_count, - Interaction& item, const size_t n_particles, const double rVerlet, - const uint8_t* __restrict__ type, const uint64_t* __restrict__ id, const VertexArray* __restrict__ vertices, shapes& shps) - { - for(size_t p = 0 ; p < n_particles ; p++) - { - const auto va = vertices[p]; - const shape* shp = shps[type[p]]; - int nv = shp->get_number_of_vertices(); - for(int sub = 0 ; sub < nv ; sub++) - { - bool contact = exaDEM::filter_vertex_driver( - driver, rVerlet, - va, sub, shp); - if(contact) - { - item.p_i = p; - item.id_i = id[p]; - driver_count[p]++; - item.sub_i = sub; - item.sub_j = -1; - driver_data.push_back(item); - } - } - } - } - - inline void execute () override final - { - auto& g = *grid; - const auto cells = g.cells(); - const size_t n_cells = g.number_of_cells(); // nbh.size(); - const IJK dims = g.dimension(); - const int gl = g.ghost_layers(); - auto & interactions = ges->m_data; - - // if grid structure (dimensions) changed, we invalidate thie whole data - if( interactions.size() != n_cells ) - { - ldbg << "number of cells changed, reset friction data" << std::endl; - interactions.clear(); - interactions.resize( n_cells ); - } - assert( interactions.size() == n_cells ); - - if( ! chunk_neighbors.has_value() ) - { + } + + void add_driver_interaction( Stl_mesh& mesh, size_t cell_a, std::vector& driver_data, std::vector& driver_count, + Interaction& item, const size_t n_particles, const double rVerlet, + const uint8_t* __restrict__ type, + const uint64_t* __restrict__ id, + const double* __restrict__ rx, + const double* __restrict__ ry, + const double* __restrict__ rz, + const VertexArray* __restrict__ vertices, + const exanb::Quaternion* __restrict__ orient, + shapes& shps) + { + assert ( cell_a < mesh.grid_indexes.size() ); + auto& list = mesh.grid_indexes[cell_a]; + const size_t stl_nv = list.vertices.size(); + const size_t stl_ne = list.edges.size(); + const size_t stl_nf = list.faces.size(); + auto& stl_shp = mesh.shp; + OBB* __restrict__ stl_obb_vertices = onika::cuda::vector_data( stl_shp.m_obb_vertices ); + OBB* __restrict__ stl_obb_edges = onika::cuda::vector_data( stl_shp.m_obb_edges ); + OBB* __restrict__ stl_obb_faces = onika::cuda::vector_data( stl_shp.m_obb_faces ); + + if( stl_nv == 0 && stl_ne == 0 && stl_nf == 0) return; + + for( size_t p = 0 ; p < n_particles ; p++ ) + { + Vec3d r = {rx[p], ry[p], rz[p]}; // position + item.p_i = p; + item.id_i = id[p]; + auto ti = type[p]; + const shape* shpi = shps[ti]; + const size_t nv = shpi->get_number_of_vertices(); + const size_t ne = shpi->get_number_of_edges(); + const size_t nf = shpi->get_number_of_faces(); + + for( size_t i = 0 ; i < nv ; i++ ) + { + vec3r v = conv_to_vec3r(vertices[p][i]); + // vertex - vertex + item.type = 7; + item.sub_i = i; + for( size_t j = 0 ; j < stl_nv ; j++ ) + { + size_t idx = list.vertices[j]; + OBB& obb = stl_obb_vertices[idx]; + obb.enlarge( shpi->m_radius ); + if ( obb.intersect( v )) + { + driver_count[p]++; + item.sub_j = idx; + driver_data.push_back(item); + } + obb.enlarge( -shpi->m_radius ); + } + // vertex - edge + item.type = 8; + for( size_t j = 0 ; j < stl_ne ; j++ ) + { + size_t idx = list.edges[j]; + OBB& obb = stl_obb_edges[idx]; + obb.enlarge( shpi->m_radius ); + if ( obb.intersect(v) ) + { + driver_count[p]++; + item.sub_j = idx; + driver_data.push_back(item); + } + obb.enlarge( -shpi->m_radius ); + } + // vertex - face + item.type = 9; + for( size_t j = 0 ; j < stl_nf ; j++ ) + { + size_t idx = list.faces[j]; + OBB& obb = stl_obb_faces[idx]; + obb.enlarge( shpi->m_radius ); + if ( obb.intersect(v) ) + { + driver_count[p]++; + item.sub_j = idx; + driver_data.push_back(item); + } + obb.enlarge( -shpi->m_radius ); + } + } + + for( size_t i = 0 ; i < ne ; i++ ) + { + item.type = 10; + item.sub_i = i; + OBB obb_edge_i = shpi->get_obb_edge(r, i, orient[i]); + // edge - edge + for(size_t j = 0; j < stl_ne ; j++) + { + const size_t idx = list.edges[j]; + OBB& stl_obb_edge = stl_obb_edges[idx]; + if( obb_edge_i.intersect(stl_obb_edge)) + { + driver_count[p]++; + item.sub_j = idx; + driver_data.push_back(item); + } + } + + // edge - vertex + item.type = 11; + for(size_t j = 0; j < stl_nv ; j++) + { + const size_t idx = list.vertices[j]; + OBB& obb = stl_obb_vertices[idx]; + if( obb_edge_i.intersect(obb)) + { + driver_count[p]++; + item.sub_j = idx; + driver_data.push_back(item); + } + } + } + + // face vertex + item.type = 12; + for( size_t i = 0 ; i < nf ; i++ ) + { + item.sub_i = i; + OBB obb_face_i = shpi->get_obb_face(r, i, orient[i]); + for(size_t j = 0; j < stl_nv ; j++) + { + size_t idx = list.vertices[j]; + OBB& obb = stl_obb_vertices[idx]; + if( obb_face_i.intersect(obb)) + { + driver_count[p]++; + item.sub_j = idx; + driver_data.push_back(item); + } + } + } + } + } + + template + void add_driver_interaction( D& driver, std::vector& driver_data, std::vector& driver_count, + Interaction& item, const size_t n_particles, const double rVerlet, + const uint8_t* __restrict__ type, const uint64_t* __restrict__ id, const VertexArray* __restrict__ vertices, shapes& shps) + { + for(size_t p = 0 ; p < n_particles ; p++) + { + const auto va = vertices[p]; + const shape* shp = shps[type[p]]; + int nv = shp->get_number_of_vertices(); + for(int sub = 0 ; sub < nv ; sub++) + { + bool contact = exaDEM::filter_vertex_driver( + driver, rVerlet, + va, sub, shp); + if(contact) + { + item.p_i = p; + item.id_i = id[p]; + driver_count[p]++; + item.sub_i = sub; + item.sub_j = -1; + driver_data.push_back(item); + } + } + } + } + + inline void execute () override final + { + auto& g = *grid; + const auto cells = g.cells(); + const size_t n_cells = g.number_of_cells(); // nbh.size(); + const IJK dims = g.dimension(); + const int gl = g.ghost_layers(); + auto & interactions = ges->m_data; + + // if grid structure (dimensions) changed, we invalidate thie whole data + if( interactions.size() != n_cells ) + { + ldbg << "number of cells has changed, reset friction data" << std::endl; + interactions.clear(); + interactions.resize( n_cells ); + } + assert( interactions.size() == n_cells ); + + if( ! chunk_neighbors.has_value() ) + { #pragma omp parallel for schedule(static) - for(size_t i = 0 ; i < n_cells ; i++) interactions[i].initialize(0); - return; - } + for(size_t i = 0 ; i < n_cells ; i++) interactions[i].initialize(0); + return; + } - auto & shps = *shapes_collection; - double rVerlet = *rcut_inc; + auto & shps = *shapes_collection; + double rVerlet = *rcut_inc; - // use OBB for vertex/edge and vertex/faces - constexpr bool skip_obb = false; + // use OBB for vertex/edge and vertex/faces + constexpr bool skip_obb = false; - auto& indexes = *idxs; - indexes.clear(); - GRID_FOR_BEGIN(dims-2*gl,_,loc) - { - IJK loc_a = loc + gl; - size_t cell_a = grid_ijk_to_index( dims , loc_a ); - const unsigned int n_particles = cells[cell_a].size(); - if( n_particles > 0 ) indexes.push_back(cell_a); - } - GRID_FOR_END + auto& indexes = *idxs; + indexes.clear(); + GRID_FOR_BEGIN(dims-2*gl,_,loc) + { + IJK loc_a = loc + gl; + size_t cell_a = grid_ijk_to_index( dims , loc_a ); + const unsigned int n_particles = cells[cell_a].size(); + if( n_particles > 0 ) indexes.push_back(cell_a); + } + GRID_FOR_END # pragma omp parallel - { - Interaction item; - std::vector local_history; - std::vector driver_data; - std::vector poly_data; - std::vector driver_count; - std::vector poly_count; - /* - GRID_OMP_FOR_BEGIN(dims-2*gl,_,block_loc, schedule(dynamic) ) - { - IJK loc_a = block_loc + gl; - size_t cell_a = grid_ijk_to_index( dims , loc_a ); - */ + { + Interaction item; + std::vector local_history; + std::vector driver_data; + std::vector poly_data; + std::vector driver_count; + std::vector poly_count; #pragma omp for schedule(dynamic) - for(size_t ci = 0 ; ci < indexes.size() ; ci++) - { - size_t cell_a = indexes[ci]; - IJK loc_a = grid_index_to_ijk( dims , cell_a ); - - const unsigned int n_particles = cells[cell_a].size(); - CellExtraDynamicDataStorageT& storage = interactions[cell_a]; - - assert ( - interaction_test::check_extra_interaction_storage_consistency( - storage.number_of_particles(), - storage.m_info.data(), - storage.m_data.data() - )); - - if(n_particles == 0) continue; - - // extract history before reset it - const size_t data_size = storage.m_data.size(); - - Interaction* __restrict__ data_ptr = storage.m_data.data(); - extract_history(local_history, data_ptr, data_size); - std::sort ( local_history.begin(), local_history.end() ); - - driver_data.clear(); - poly_data.clear(); - driver_count.assign(n_particles, 0); - poly_count.assign(n_particles, 0); - - const uint64_t* __restrict__ id_a = cells[cell_a][ field::id ]; ONIKA_ASSUME_ALIGNED(id_a); - const auto* __restrict__ rx_a = cells[cell_a][ field::rx ]; ONIKA_ASSUME_ALIGNED(rx_a); - const auto* __restrict__ ry_a = cells[cell_a][ field::ry ]; ONIKA_ASSUME_ALIGNED(ry_a); - const auto* __restrict__ rz_a = cells[cell_a][ field::rz ]; ONIKA_ASSUME_ALIGNED(rz_a); - const auto* __restrict__ t_a = cells[cell_a][ field::type ]; ONIKA_ASSUME_ALIGNED(t_a); - const auto* __restrict__ orient_a = cells[cell_a][ field::orient ]; ONIKA_ASSUME_ALIGNED(orient_a); - const auto* __restrict__ vertices_a = cells[cell_a][ field::vertices ]; ONIKA_ASSUME_ALIGNED(vertices_a); - - storage.initialize(n_particles); - auto& info_particles = storage.m_info; - - auto add_contact = []( std::vector& list, Interaction& item, int sub_i, int sub_j) -> void - { - item.sub_i = sub_i; - item.sub_j = sub_j; - list.push_back(item); - }; - - auto incr_particle_interactions = [] (std::vector& count, int p_a) - { - count[p_a]++; - }; - - // get particle id - for( size_t it = 0 ; it < n_particles ; it++) - { - std::get<2> (info_particles[it]) = id_a[it]; - } - - // First drivers - if ( drivers.has_value() ) - { - auto& drvs = *drivers; - item.cell_i = cell_a; - item.id_j = -1; - item.cell_j = -1; - item.p_j = -1; - item.moment = Vec3d{0,0,0}; - item.friction = Vec3d{0,0,0}; - for( size_t drvs_idx = 0 ; drvs_idx < drvs.get_size() ; drvs_idx++ ) - { - item.id_j = drvs_idx; // we store the driver idx - if (drvs.type(drvs_idx) == DRIVER_TYPE::CYLINDER) - { - item.type = 4; - Cylinder driver = std::get(drvs.data(drvs_idx)) ; - add_driver_interaction( driver, driver_data, driver_count, - item, n_particles, rVerlet, - t_a, id_a, vertices_a, shps); - } - else if ( drvs.type(drvs_idx) == DRIVER_TYPE::SURFACE) - { - item.type = 5; - Surface driver = std::get(drvs.data(drvs_idx)) ; - add_driver_interaction( driver, driver_data, driver_count, - item, n_particles, rVerlet, - t_a, id_a, vertices_a, shps); - } - else if ( drvs.type(drvs_idx) == DRIVER_TYPE::BALL) - { - item.type = 6; - Ball driver = std::get(drvs.data(drvs_idx)) ; - add_driver_interaction( driver, driver_data, driver_count, - item, n_particles, rVerlet, - t_a, id_a, vertices_a, shps); - } - } - } - - // Second polyhedra - apply_cell_particle_neighbors(*grid, *chunk_neighbors, cell_a, loc_a, std::false_type() /* not symetric */, - [&g , cells, &info_particles, cell_a, &poly_data, &poly_count, &item, &shps, rVerlet, id_a, rx_a, ry_a, rz_a, t_a, orient_a, vertices_a, &add_contact, &incr_particle_interactions] - ( int p_a, size_t cell_b, unsigned int p_b , size_t p_nbh_index ){ - // default value of the interaction studied (A or i -> B or j) - const uint64_t id_nbh = cells[cell_b][field::id][p_b]; - if( id_a[p_a] >= id_nbh) - { - if ( !g.is_ghost_cell(cell_b) ) return; - } - - const uint8_t type_nbh = cells[cell_b][field::type][p_b]; - const Quaternion orient_nbh = cells[cell_b][field::orient][p_b]; - const double rx_nbh = cells[cell_b][field::rx][p_b]; - const double ry_nbh = cells[cell_b][field::ry][p_b]; - const double rz_nbh = cells[cell_b][field::rz][p_b]; - const auto& vertices_b = cells[cell_b][field::vertices][p_b]; - - // prev - const shape* shp = shps[t_a[p_a]]; - const shape* shp_nbh = shps[type_nbh]; - OBB obb_i = shp->obb; - OBB obb_j = shp_nbh->obb; - const Quaternion& orient = orient_a[p_a]; - const double rx = rx_a[p_a]; - const double ry = ry_a[p_a]; - const double rz = rz_a[p_a]; - quat conv_orient_i = quat{vec3r{orient.x, orient.y, orient.z}, orient.w}; - quat conv_orient_j = quat{vec3r{orient_nbh.x, orient_nbh.y, orient_nbh.z}, orient_nbh.w}; - obb_i.rotate(conv_orient_i); - obb_j.rotate(conv_orient_j); - obb_i.translate(vec3r{rx, ry, rz}); - obb_j.translate(vec3r{rx_nbh, ry_nbh, rz_nbh}); - - obb_i.enlarge(rVerlet); - obb_j.enlarge(rVerlet); - - if ( ! obb_i.intersect(obb_j) ) return; - - // reset rVerlet - obb_i . enlarge(-rVerlet); - obb_j . enlarge(-rVerlet); - - // Add interactions - item.id_i = id_a[p_a]; - item.p_i = p_a; - - item.cell_i = cell_a; - item.p_j = p_b; - item.cell_j = cell_b; - - const Vec3d r = {rx, ry, rz}; - const Vec3d r_nbh = {rx_nbh, ry_nbh, rz_nbh}; - - // get particle j data. - const int nv = shp->get_number_of_vertices(); - const int ne = shp->get_number_of_edges(); - const int nf = shp->get_number_of_faces(); - const int nv_nbh = shp_nbh->get_number_of_vertices(); - const int ne_nbh = shp_nbh->get_number_of_edges(); - const int nf_nbh = shp_nbh->get_number_of_faces(); - - item.id_j = id_nbh; - // exclude possibilities with obb - for( int i = 0 ; i < nv ; i++) - { - auto vi = shp -> get_vertex (i, r, orient); - OBB obbvi; - obbvi.center = {vi.x, vi.y, vi.z}; - obbvi.enlarge(shp->m_radius + rVerlet); - if (obb_j.intersect( obbvi ) ) - { - item.type = 0; // === Vertex - Vertex - for(int j = 0; j < nv_nbh ; j++) - { - bool contact = exaDEM::filter_vertex_vertex(rVerlet, vertices_a[p_a], i, shp, vertices_b, j, shp_nbh); - if ( contact ) - { - incr_particle_interactions(poly_count, p_a); - add_contact(poly_data, item, i, j); - } - } - - item.type = 1; // === vertex edge - for(int j = 0; j < ne_nbh ; j++) - { - bool contact = exaDEM::filter_vertex_edge (obbvi, r_nbh, j, shp_nbh, orient_nbh); - if(contact) - { - incr_particle_interactions(poly_count, p_a); - add_contact(poly_data, item, i, j); - } - } - - item.type = 2; // === vertex face - for(int j = 0; j < nf_nbh ; j++) - { - bool contact = exaDEM::filter_vertex_face (obbvi, r_nbh, j, shp_nbh, orient_nbh); - if(contact) - { - incr_particle_interactions(poly_count, p_a); - add_contact(poly_data, item, i, j); - } - } - } - } - - item.type = 3; // === edge edge - for( int i = 0 ; i < ne ; i++) - { - OBB obb_edge_i = shp->get_obb_edge(r, i, orient); - if ( obb_j.intersect (obb_edge_i) ) - { - obb_edge_i.enlarge(rVerlet); - for(int j = 0; j < ne_nbh ; j++) - { - OBB obb_edge_j = shp_nbh->get_obb_edge(r_nbh, j, orient_nbh); - if( obb_edge_i.intersect(obb_edge_j)) - { - incr_particle_interactions(poly_count, p_a); - add_contact(poly_data, item, i, j); - } - } - } - } - - item.cell_j = cell_a; - item.id_j = id_a[p_a]; - item.p_j = p_a; - - item.cell_i = cell_b; - item.p_i = p_b; - item.id_i = id_nbh; - - for( int j = 0 ; j < nv_nbh ; j++) - { - auto vj = shp -> get_vertex (j, r_nbh, orient_nbh); - OBB obbvj; - obbvj.center = {vj.x, vj.y, vj.z}; - obbvj.enlarge(shp_nbh->m_radius + rVerlet); - - if (obb_i.intersect( obbvj ) ) - { - item.type = 1; // === vertex edge - for(int i = 0; i < ne ; i++) - { - bool contact = exaDEM::filter_vertex_edge (obbvj, r, i, shp, orient); - if(contact) - { - incr_particle_interactions(poly_count, p_a); - add_contact(poly_data, item, j, i); - } - } - - item.type = 2; // === vertex face - for(int i = 0; i < nf ; i++) - { - bool contact = exaDEM::filter_vertex_face (obbvj, r, i, shp, orient); - if(contact) - { - incr_particle_interactions(poly_count, p_a); - add_contact(poly_data, item, j, i); - } - } - } - } - }); - - // - update_friction_moment(driver_data, local_history); - update_friction_moment(poly_data, local_history); - - // build storage - size_t offset = 0 ; - size_t offset_driver = 0 ; - size_t offset_poly = 0 ; - storage.m_data.resize( poly_data.size () + driver_data.size() ); - data_ptr = storage.m_data.data(); - for(size_t p = 0 ; p < n_particles ; p++) - { - std::get<0> (info_particles[p]) = offset; - for(size_t it = 0 ; it < driver_count[p] ; it++ ) data_ptr[offset ++ ] = driver_data[offset_driver++]; - for(size_t it = 0 ; it < poly_count[p] ; it++ ) data_ptr[offset ++ ] = poly_data[ offset_poly++]; - std::get<1> (info_particles[p]) = driver_count[p] + poly_count[p]; - } - assert ( offset_driver == driver_data.size() ); - assert ( offset_poly == poly_data.size() ); - - // add history, local history and local have to be sorted. local is sorted by construction. - assert ( - interaction_test::check_extra_interaction_storage_consistency( - storage.number_of_particles(), - storage.m_info.data(), - storage.m_data.data() - )); - } -// GRID_OMP_FOR_END - } + for(size_t ci = 0 ; ci < indexes.size() ; ci++) + { + size_t cell_a = indexes[ci]; + IJK loc_a = grid_index_to_ijk( dims , cell_a ); + + const unsigned int n_particles = cells[cell_a].size(); + CellExtraDynamicDataStorageT& storage = interactions[cell_a]; + + assert ( + interaction_test::check_extra_interaction_storage_consistency( + storage.number_of_particles(), + storage.m_info.data(), + storage.m_data.data() + )); + + if(n_particles == 0) continue; + + // extract history before reset it + const size_t data_size = storage.m_data.size(); + + Interaction* __restrict__ data_ptr = storage.m_data.data(); + extract_history(local_history, data_ptr, data_size); + std::sort ( local_history.begin(), local_history.end() ); + + driver_data.clear(); + poly_data.clear(); + driver_count.assign(n_particles, 0); + poly_count.assign(n_particles, 0); + + const uint64_t* __restrict__ id_a = cells[cell_a][ field::id ]; ONIKA_ASSUME_ALIGNED(id_a); + const auto* __restrict__ rx_a = cells[cell_a][ field::rx ]; ONIKA_ASSUME_ALIGNED(rx_a); + const auto* __restrict__ ry_a = cells[cell_a][ field::ry ]; ONIKA_ASSUME_ALIGNED(ry_a); + const auto* __restrict__ rz_a = cells[cell_a][ field::rz ]; ONIKA_ASSUME_ALIGNED(rz_a); + const auto* __restrict__ t_a = cells[cell_a][ field::type ]; ONIKA_ASSUME_ALIGNED(t_a); + const auto* __restrict__ orient_a = cells[cell_a][ field::orient ]; ONIKA_ASSUME_ALIGNED(orient_a); + const auto* __restrict__ vertices_a = cells[cell_a][ field::vertices ]; ONIKA_ASSUME_ALIGNED(vertices_a); + + storage.initialize(n_particles); + auto& info_particles = storage.m_info; + + auto add_contact = []( std::vector& list, Interaction& item, int sub_i, int sub_j) -> void + { + item.sub_i = sub_i; + item.sub_j = sub_j; + list.push_back(item); + }; + + auto incr_particle_interactions = [] (std::vector& count, int p_a) + { + count[p_a]++; + }; + + // get particle id + for( size_t it = 0 ; it < n_particles ; it++) + { + std::get<2> (info_particles[it]) = id_a[it]; + } + + // First drivers + if ( drivers.has_value() ) + { + auto& drvs = *drivers; + item.cell_i = cell_a; + item.id_j = -1; + item.cell_j = -1; + item.p_j = -1; + item.moment = Vec3d{0,0,0}; + item.friction = Vec3d{0,0,0}; + for( size_t drvs_idx = 0 ; drvs_idx < drvs.get_size() ; drvs_idx++ ) + { + item.id_j = drvs_idx; // we store the driver idx + if (drvs.type(drvs_idx) == DRIVER_TYPE::CYLINDER) + { + item.type = 4; + Cylinder& driver = std::get(drvs.data(drvs_idx)) ; + add_driver_interaction( driver, driver_data, driver_count, + item, n_particles, rVerlet, + t_a, id_a, vertices_a, shps); + } + else if ( drvs.type(drvs_idx) == DRIVER_TYPE::SURFACE) + { + item.type = 5; + Surface& driver = std::get(drvs.data(drvs_idx)) ; + add_driver_interaction( driver, driver_data, driver_count, + item, n_particles, rVerlet, + t_a, id_a, vertices_a, shps); + } + else if ( drvs.type(drvs_idx) == DRIVER_TYPE::BALL) + { + item.type = 6; + Ball& driver = std::get(drvs.data(drvs_idx)) ; + add_driver_interaction( driver, driver_data, driver_count, + item, n_particles, rVerlet, + t_a, id_a, vertices_a, shps); + } + else if ( drvs.type(drvs_idx) == DRIVER_TYPE::STL_MESH) + { + Stl_mesh& driver = std::get(drvs.data(drvs_idx)) ; + //driver.grid_indexes_summary(); + add_driver_interaction ( driver, cell_a, driver_data, driver_count, + item, n_particles, rVerlet, + t_a, id_a, rx_a, ry_a, rz_a, vertices_a, orient_a, shps); + } + } + } + + // Second polyhedra + apply_cell_particle_neighbors(*grid, *chunk_neighbors, cell_a, loc_a, std::false_type() /* not symetric */, + [&g , cells, &info_particles, cell_a, &poly_data, &poly_count, &item, &shps, rVerlet, id_a, rx_a, ry_a, rz_a, t_a, orient_a, vertices_a, &add_contact, &incr_particle_interactions] + ( int p_a, size_t cell_b, unsigned int p_b , size_t p_nbh_index ){ + // default value of the interaction studied (A or i -> B or j) + const uint64_t id_nbh = cells[cell_b][field::id][p_b]; + if( id_a[p_a] >= id_nbh) + { + if ( !g.is_ghost_cell(cell_b) ) return; + } + + const uint8_t type_nbh = cells[cell_b][field::type][p_b]; + const Quaternion orient_nbh = cells[cell_b][field::orient][p_b]; + const double rx_nbh = cells[cell_b][field::rx][p_b]; + const double ry_nbh = cells[cell_b][field::ry][p_b]; + const double rz_nbh = cells[cell_b][field::rz][p_b]; + const auto& vertices_b = cells[cell_b][field::vertices][p_b]; + + // prev + const shape* shp = shps[t_a[p_a]]; + const shape* shp_nbh = shps[type_nbh]; + OBB obb_i = shp->obb; + OBB obb_j = shp_nbh->obb; + const Quaternion& orient = orient_a[p_a]; + const double rx = rx_a[p_a]; + const double ry = ry_a[p_a]; + const double rz = rz_a[p_a]; + quat conv_orient_i = quat{vec3r{orient.x, orient.y, orient.z}, orient.w}; + quat conv_orient_j = quat{vec3r{orient_nbh.x, orient_nbh.y, orient_nbh.z}, orient_nbh.w}; + obb_i.rotate(conv_orient_i); + obb_j.rotate(conv_orient_j); + obb_i.translate(vec3r{rx, ry, rz}); + obb_j.translate(vec3r{rx_nbh, ry_nbh, rz_nbh}); + + obb_i.enlarge(rVerlet); + obb_j.enlarge(rVerlet); + + if ( ! obb_i.intersect(obb_j) ) return; + + // reset rVerlet + obb_i . enlarge(-rVerlet); + obb_j . enlarge(-rVerlet); + + // Add interactions + item.id_i = id_a[p_a]; + item.p_i = p_a; + + item.cell_i = cell_a; + item.p_j = p_b; + item.cell_j = cell_b; + + const Vec3d r = {rx, ry, rz}; + const Vec3d r_nbh = {rx_nbh, ry_nbh, rz_nbh}; + + // get particle j data. + const int nv = shp->get_number_of_vertices(); + const int ne = shp->get_number_of_edges(); + const int nf = shp->get_number_of_faces(); + const int nv_nbh = shp_nbh->get_number_of_vertices(); + const int ne_nbh = shp_nbh->get_number_of_edges(); + const int nf_nbh = shp_nbh->get_number_of_faces(); + + item.id_j = id_nbh; + // exclude possibilities with obb + for( int i = 0 ; i < nv ; i++) + { + auto vi = shp -> get_vertex (i, r, orient); + OBB obbvi; + obbvi.center = {vi.x, vi.y, vi.z}; + obbvi.enlarge(shp->m_radius + rVerlet); + if (obb_j.intersect( obbvi ) ) + { + item.type = 0; // === Vertex - Vertex + for(int j = 0; j < nv_nbh ; j++) + { + bool contact = exaDEM::filter_vertex_vertex(rVerlet, vertices_a[p_a], i, shp, vertices_b, j, shp_nbh); + if ( contact ) + { + incr_particle_interactions(poly_count, p_a); + add_contact(poly_data, item, i, j); + } + } + + item.type = 1; // === vertex edge + for(int j = 0; j < ne_nbh ; j++) + { + bool contact = exaDEM::filter_vertex_edge (obbvi, r_nbh, j, shp_nbh, orient_nbh); + if(contact) + { + incr_particle_interactions(poly_count, p_a); + add_contact(poly_data, item, i, j); + } + } + + item.type = 2; // === vertex face + for(int j = 0; j < nf_nbh ; j++) + { + bool contact = exaDEM::filter_vertex_face (obbvi, r_nbh, j, shp_nbh, orient_nbh); + if(contact) + { + incr_particle_interactions(poly_count, p_a); + add_contact(poly_data, item, i, j); + } + } + } + } + + item.type = 3; // === edge edge + for( int i = 0 ; i < ne ; i++) + { + OBB obb_edge_i = shp->get_obb_edge(r, i, orient); + if ( obb_j.intersect (obb_edge_i) ) + { + obb_edge_i.enlarge(rVerlet); + for(int j = 0; j < ne_nbh ; j++) + { + OBB obb_edge_j = shp_nbh->get_obb_edge(r_nbh, j, orient_nbh); + if( obb_edge_i.intersect(obb_edge_j)) + { + incr_particle_interactions(poly_count, p_a); + add_contact(poly_data, item, i, j); + } + } + } + } + + item.cell_j = cell_a; + item.id_j = id_a[p_a]; + item.p_j = p_a; + + item.cell_i = cell_b; + item.p_i = p_b; + item.id_i = id_nbh; + + for( int j = 0 ; j < nv_nbh ; j++) + { + auto vj = shp -> get_vertex (j, r_nbh, orient_nbh); + OBB obbvj; + obbvj.center = {vj.x, vj.y, vj.z}; + obbvj.enlarge(shp_nbh->m_radius + rVerlet); + + if (obb_i.intersect( obbvj ) ) + { + item.type = 1; // === vertex edge + for(int i = 0; i < ne ; i++) + { + bool contact = exaDEM::filter_vertex_edge (obbvj, r, i, shp, orient); + if(contact) + { + incr_particle_interactions(poly_count, p_a); + add_contact(poly_data, item, j, i); + } + } + + item.type = 2; // === vertex face + for(int i = 0; i < nf ; i++) + { + bool contact = exaDEM::filter_vertex_face (obbvj, r, i, shp, orient); + if(contact) + { + incr_particle_interactions(poly_count, p_a); + add_contact(poly_data, item, j, i); + } + } + } + } + }); + + // + update_friction_moment(driver_data, local_history); + update_friction_moment(poly_data, local_history); + + // build storage + size_t offset = 0 ; + size_t offset_driver = 0 ; + size_t offset_poly = 0 ; + storage.m_data.resize( poly_data.size () + driver_data.size() ); + data_ptr = storage.m_data.data(); + for(size_t p = 0 ; p < n_particles ; p++) + { + std::get<0> (info_particles[p]) = offset; + for(size_t it = 0 ; it < driver_count[p] ; it++ ) data_ptr[offset ++ ] = driver_data[offset_driver++]; + for(size_t it = 0 ; it < poly_count[p] ; it++ ) data_ptr[offset ++ ] = poly_data[ offset_poly++]; + std::get<1> (info_particles[p]) = driver_count[p] + poly_count[p]; + } + assert ( offset_driver == driver_data.size() ); + assert ( offset_poly == poly_data.size() ); + + // add history, local history and local have to be sorted. local is sorted by construction. + assert ( + interaction_test::check_extra_interaction_storage_consistency( + storage.number_of_particles(), + storage.m_info.data(), + storage.m_data.data() + )); + } + // GRID_OMP_FOR_END + } + } + }; + + template using UpdateGridCellInteractionTmpl = UpdateGridCellInteraction; + + // === register factories === + CONSTRUCTOR_FUNCTION + { + OperatorNodeFactory::instance()->register_factory( "update_grid_interaction", make_grid_variant_operator< UpdateGridCellInteraction > ); } - }; - - template using UpdateGridCellInteractionTmpl = UpdateGridCellInteraction; - - // === register factories === - CONSTRUCTOR_FUNCTION - { - OperatorNodeFactory::instance()->register_factory( "update_grid_interaction", make_grid_variant_operator< UpdateGridCellInteraction > ); - } - } +} From 54fd454a968ce03dec82f16e73a278dc0991dcf2 Mon Sep 17 00:00:00 2001 From: rprat-pro Date: Thu, 28 Mar 2024 09:36:04 +0100 Subject: [PATCH 2/4] fix issue polyhedron face vs stl vertex --- src/polyhedra/include/exaDEM/compute_hooke_interaction.h | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/polyhedra/include/exaDEM/compute_hooke_interaction.h b/src/polyhedra/include/exaDEM/compute_hooke_interaction.h index a451243..77681ac 100644 --- a/src/polyhedra/include/exaDEM/compute_hooke_interaction.h +++ b/src/polyhedra/include/exaDEM/compute_hooke_interaction.h @@ -87,13 +87,12 @@ namespace exaDEM const Vec3d r_i = { cell[field::rx][p_i], cell[field::ry][p_i], cell[field::rz][p_i] }; // === vrot const Vec3d& vrot_i = cell[field::vrot][p_i]; - const auto& t_i = cell[field::type][p_i]; const Quaternion& orient_i = cell[field::orient][p_i]; const auto& shp_j = driver.shp; // WARNING const Quaternion orient_j = {1.0,0.0,0.0,0.0}; - auto [contact, dn, n, contact_position] = func(I.type, r_i, sub_i, &shp_i[t_i], orient_i, driver.center, sub_j, &shp_j, orient_j); + auto [contact, dn, n, contact_position] = func(I.type, r_i, sub_i, shp_i, orient_i, driver.center, sub_j, &shp_j, orient_j); if(contact) { From aebf70c3d24eb1421b2b98c4bdbe979f05a28470 Mon Sep 17 00:00:00 2001 From: rprat-pro Date: Thu, 28 Mar 2024 09:36:43 +0100 Subject: [PATCH 3/4] add type field in paraview outputs --- data/config/config_iteration_dump.msp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/data/config/config_iteration_dump.msp b/data/config/config_iteration_dump.msp index 7c682cc..01a350a 100644 --- a/data/config/config_iteration_dump.msp +++ b/data/config/config_iteration_dump.msp @@ -23,7 +23,7 @@ write_paraview_generic: binary: false write_ghost: false write_box: true - fields: ["vx","vy","vz","id","orient"] + fields: ["type","vx","vy","vz","id","orient"] # fields: ["vx","vy","vz","id","radius"] dump_data_paraview: From a508f333273a82bf6f0d7d4f580d01c8a05cb6f4 Mon Sep 17 00:00:00 2001 From: rprat-pro Date: Thu, 28 Mar 2024 10:41:05 +0100 Subject: [PATCH 4/4] add example for polyhedra + mesh or cylinder --- .../polyhedra/rotating_drum/CMakeLists.txt | 9 +++ example/polyhedra/rotating_drum/alpha3.shp | 25 +++++++ .../polyhedra/rotating_drum/octahedron.shp | 44 +++++++++++ .../rotating-drum-mixte-mini.msp | 60 +++++++++++++++ .../rotating_drum/rotating-drum-mixte.msp | 60 +++++++++++++++ example/polyhedra/stl_mesh/CMakeLists.txt | 8 ++ example/polyhedra/stl_mesh/box_mixte.stl | 73 +++++++++++++++++++ .../polyhedra/stl_mesh/stl_mesh_box_mixte.msp | 59 +++++++++++++++ .../stl_mesh/stl_mesh_box_mixte_mini.msp | 59 +++++++++++++++ 9 files changed, 397 insertions(+) create mode 100644 example/polyhedra/rotating_drum/alpha3.shp create mode 100644 example/polyhedra/rotating_drum/octahedron.shp create mode 100644 example/polyhedra/rotating_drum/rotating-drum-mixte-mini.msp create mode 100644 example/polyhedra/rotating_drum/rotating-drum-mixte.msp create mode 100644 example/polyhedra/stl_mesh/box_mixte.stl create mode 100644 example/polyhedra/stl_mesh/stl_mesh_box_mixte.msp create mode 100644 example/polyhedra/stl_mesh/stl_mesh_box_mixte_mini.msp diff --git a/example/polyhedra/rotating_drum/CMakeLists.txt b/example/polyhedra/rotating_drum/CMakeLists.txt index 79f786d..ff598a2 100644 --- a/example/polyhedra/rotating_drum/CMakeLists.txt +++ b/example/polyhedra/rotating_drum/CMakeLists.txt @@ -3,6 +3,15 @@ set(ExampleName "OctahedraRotatingDrum") set(FileName "rotating-drum.msp") file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/shapes.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/alpha3.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/octahedron.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) + +add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) +add_test(Example_ExaDEM_${ExampleName}_omp ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 2) +add_test(Example_ExaDEM_${ExampleName}_mpi mpirun -n 2 ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) + +set(ExampleName "OctahedraHexapodsRotatingDrum") +set(FileName "rotating-drum-mixte-mini.msp") add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) add_test(Example_ExaDEM_${ExampleName}_omp ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 2) diff --git a/example/polyhedra/rotating_drum/alpha3.shp b/example/polyhedra/rotating_drum/alpha3.shp new file mode 100644 index 0000000..54f1e1f --- /dev/null +++ b/example/polyhedra/rotating_drum/alpha3.shp @@ -0,0 +1,25 @@ + +< +name alpha3 +radius 0.237770037899799 +preCompDone y +nv 6 +0.475540075799599 0 0 +-0.475540075799599 0 0 +0 0.475540075799599 0 +0 -0.475540075799599 0 +0 0 0.475540075799599 +0 0 -0.475540075799599 +ne 3 +0 1 +2 3 +4 5 +nf 0 +obb.extent 0.713310113699398 0.713310113699398 0.713310113699398 +obb.e1 1 0 0 +obb.e2 0 1 0 +obb.e3 0 0 1 +obb.center 0 0 0 +volume 0.523598775598299 +I/m 0.13266 0.13266 0.13266 +> diff --git a/example/polyhedra/rotating_drum/octahedron.shp b/example/polyhedra/rotating_drum/octahedron.shp new file mode 100644 index 0000000..82acda4 --- /dev/null +++ b/example/polyhedra/rotating_drum/octahedron.shp @@ -0,0 +1,44 @@ + +< +name Octahedron +radius 0.1 +preCompDone y +nv 6 +0.2310789034541148 -0.2310789034541148 0.0 +0.2310789034541148 0.2310789034541148 0.0 +0.0 0.0 0.32679491924311227 +-0.2310789034541148 -0.2310789034541148 0.0 +-0.2310789034541148 0.2310789034541148 0.0 +0.0 0.0 -0.32679491924311227 +ne 12 +0 1 +2 1 +2 0 +0 3 +2 3 +3 4 +4 2 +4 1 +5 0 +5 1 +5 4 +5 3 +nf 8 +3 0 1 2 +3 2 3 4 +3 1 2 4 +3 0 2 3 +3 0 5 1 +3 0 5 3 +3 3 5 4 +3 4 5 1 +obb.extent 0.33107890345411484 0.33107890345411484 0.4267949192431123 +obb.e1 1.0 0.0 0.0 +obb.e2 0.0 1.0 0.0 +obb.e3 0.0 0.0 1.0 +obb.center 0.0 0.0 0.0 +position 0.0 0.0 0.0 +orientation 1.0 0.0 0.0 0.0 +volume 0.16666666666666666 +I/m 0.04999999999999999 0.04999999999999999 0.04999999999999999 +> diff --git a/example/polyhedra/rotating_drum/rotating-drum-mixte-mini.msp b/example/polyhedra/rotating_drum/rotating-drum-mixte-mini.msp new file mode 100644 index 0000000..89817c5 --- /dev/null +++ b/example/polyhedra/rotating_drum/rotating-drum-mixte-mini.msp @@ -0,0 +1,60 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: BCC + types: [ 0 , 1 ] + size: [ 1.5 , 1.5 , 1.5 ] + repeats: [ 5 , 5 , 5 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: alpha3.shp + - read_shape_file: + filename: octahedron.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_cylinder: + id: 0 + center: [7.5, 7.5, 7.5] + axis: [1, 0, 1] + radius: 13 + angular_velocity: [0,0.2,0] + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.5, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,true,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 10000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: -1 + dt: 0.00005 s + rcut_inc: 0.1 m + enable_stl_mesh: true diff --git a/example/polyhedra/rotating_drum/rotating-drum-mixte.msp b/example/polyhedra/rotating_drum/rotating-drum-mixte.msp new file mode 100644 index 0000000..a7907f7 --- /dev/null +++ b/example/polyhedra/rotating_drum/rotating-drum-mixte.msp @@ -0,0 +1,60 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: BCC + types: [ 0 , 1 ] + size: [ 1.5 , 1.5 , 1.5 ] + repeats: [ 30 , 30 , 30 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: alpha3.shp + - read_shape_file: + filename: octahedron.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_cylinder: + id: 0 + center: [22.5, 22.5, 22.5] + axis: [1, 0, 1] + radius: 39 + angular_velocity: [0,0.05,0] + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.5, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,true,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 1000000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: 5000 + dt: 0.0001 s + rcut_inc: 0.05 m + enable_stl_mesh: true diff --git a/example/polyhedra/stl_mesh/CMakeLists.txt b/example/polyhedra/stl_mesh/CMakeLists.txt index 46ab16a..fcd55d1 100644 --- a/example/polyhedra/stl_mesh/CMakeLists.txt +++ b/example/polyhedra/stl_mesh/CMakeLists.txt @@ -5,6 +5,7 @@ set(FileName "stl_mesh_box_hexapod_mini.msp") file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/alpha3.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/octahedron.shp DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/box.stl DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/box_mixte.stl DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/box_for_octa.stl DESTINATION ${CMAKE_CURRENT_BINARY_DIR}) add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) @@ -17,3 +18,10 @@ set(FileName "stl_mesh_box_octahedron_mini.msp") add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) add_test(Example_ExaDEM_${ExampleName}_omp ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 2) add_test(Example_ExaDEM_${ExampleName}_mpi mpirun -n 2 ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) + +set(ExampleName "STLMeshOctahedraHexapods") +set(FileName "stl_mesh_box_mixte_mini.msp") + +add_test(Example_ExaDEM_${ExampleName} ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) +add_test(Example_ExaDEM_${ExampleName}_omp ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 2) +add_test(Example_ExaDEM_${ExampleName}_mpi mpirun -n 2 ${CMAKE_BINARY_DIR}/exaDEM ${CMAKE_CURRENT_SOURCE_DIR}/${FileName} --omp_num_threads 1) diff --git a/example/polyhedra/stl_mesh/box_mixte.stl b/example/polyhedra/stl_mesh/box_mixte.stl new file mode 100644 index 0000000..2d3f23c --- /dev/null +++ b/example/polyhedra/stl_mesh/box_mixte.stl @@ -0,0 +1,73 @@ +solid Created by Gmsh +facet normal -1 0 0 + outer loop + vertex -1 -1 19 + vertex -1 19 -1 + vertex -1 -1 -1 + endloop +endfacet +facet normal -1 0 0 + outer loop + vertex -1 -1 19 + vertex -1 19 19 + vertex -1 19 -1 + endloop +endfacet +facet normal 1 0 0 + outer loop + vertex 19 -1 19 + vertex 19 -1 -1 + vertex 19 19 -1 + endloop +endfacet +facet normal 1 -0 0 + outer loop + vertex 19 -1 19 + vertex 19 19 -1 + vertex 19 19 19 + endloop +endfacet +facet normal 0 -1 -0 + outer loop + vertex 19 -1 19 + vertex -1 -1 -1 + vertex 19 -1 -1 + endloop +endfacet +facet normal -0 -1 0 + outer loop + vertex 19 -1 19 + vertex -1 -1 19 + vertex -1 -1 -1 + endloop +endfacet +facet normal 0 1 0 + outer loop + vertex 19 19 19 + vertex 19 19 -1 + vertex -1 19 -1 + endloop +endfacet +facet normal 0 1 0 + outer loop + vertex 19 19 19 + vertex -1 19 -1 + vertex -1 19 19 + endloop +endfacet +facet normal -0 0 -1 + outer loop + vertex 19 19 -1 + vertex -1 -1 -1 + vertex -1 19 -1 + endloop +endfacet +facet normal 0 -0 -1 + outer loop + vertex 19 19 -1 + vertex 19 -1 -1 + vertex -1 -1 -1 + endloop +endfacet + +endsolid Created by Gmsh diff --git a/example/polyhedra/stl_mesh/stl_mesh_box_mixte.msp b/example/polyhedra/stl_mesh/stl_mesh_box_mixte.msp new file mode 100644 index 0000000..0cc04ec --- /dev/null +++ b/example/polyhedra/stl_mesh/stl_mesh_box_mixte.msp @@ -0,0 +1,59 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: BCC + types: [ 0 , 1 ] + size: [ 1.5 , 1.5 , 1.5 ] + repeats: [ 12 , 12 , 40 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: alpha3.shp + - read_shape_file: + filename: octahedron.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_stl_mesh: + id: 0 + filename: box_mixte.stl + minskowski: 0.01 + - update_grid_stl_mesh + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.5, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,false,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 100000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: 1000 + dt: 0.0001 s + rcut_inc: 0.1 m + enable_stl_mesh: true diff --git a/example/polyhedra/stl_mesh/stl_mesh_box_mixte_mini.msp b/example/polyhedra/stl_mesh/stl_mesh_box_mixte_mini.msp new file mode 100644 index 0000000..ca686f9 --- /dev/null +++ b/example/polyhedra/stl_mesh/stl_mesh_box_mixte_mini.msp @@ -0,0 +1,59 @@ +grid_flavor: grid_flavor_dem + +includes: + - config_polyhedra.msp + +input_data: + - lattice: + structure: BCC + types: [ 0 , 1 ] + size: [ 1.5 , 1.5 , 1.5 ] + repeats: [ 5 , 5 , 5 ] + enlarge_bounds: 0.0 m + - read_shape_file: + filename: alpha3.shp + - read_shape_file: + filename: octahedron.shp + - polyhedra_define_radius + - polyhedra_set_density + - set_rand_velocity: + var: 0.1 + mean: [0.0,0.0,0.0] + - set_rand_vrot_arot + - set_quaternion + - polyhedra_update_inertia + + +input_data: + - reader1 + +setup_drivers: + - add_stl_mesh: + id: 0 + filename: box_mixte.stl + minskowski: 0.01 + - update_grid_stl_mesh + + +compute_force: + - gravity_force + - compute_hooke_interaction: + config: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.1, damp_rate: 0.999} + config_driver: { rcut: 0.0 m , dncut: 0.0 m, kn: 10000, kt: 8000, kr: 0.0, fc: 0.0, mu: 0.5, damp_rate: 0.999} + +domain: + cell_size: 2 m + periodic: [false,false,false] + +write_vtklegacy: + ascii: true + ghost: false + +global: + simulation_dump_frequency: -1 + simulation_end_iteration: 12000 + simulation_log_frequency: 1000 + simulation_paraview_frequency: -1 + dt: 0.0001 s + rcut_inc: 0.05 m + enable_stl_mesh: true