====================== AGX Granular ====================== .. _granular_body_system: ------------- Introduction ------------- ``agx::Physics::GranularBodySystem`` is a extension of the ``agx::ParticleSystem`` class that contains spherical particles, called granular bodies, with 6 degrees of freedom (6-DOF) with Hertzian contact elasticity laws, Coulomb friction and rolling resistance. AGX Granular is based on a non-smooth dynamics formulation of the discrete element method (DEM) called non-smooth discrete element method (NDEM) [#f2]_. The granular system is solved together with the rest of AGX Dynamics which allows for a strong-coupling between articulated rigid body systems and granular bodies. .. note:: In this chapter, the term ``particle`` and ``granular body`` ( or just ``body`` ) will be used interchangeably. They will all refer to the 6-DOF particle instance mentioned in the introduction above. A ``GranularBodyPtr`` is an extension of a ``ParticlePtr`` with additional buffers for inertia, rotation and angular velocity. .. image:: images/granular_machine.png :scale: 60 % :align: center **Figure 1:** Example showcasing a granular body flow interacting with and driving a mechanical system of constrained rigid bodies. A contact situation is illustrated in Fig 1. Contact forces and velocities can be decomposed in the direction of the contact normal :math:`\pmb{\vec{n}}` and tangents :math:`\pmb{\vec{t}}`. The gap function :math:`g(\pmb{x})` measures the magnitude of the overlap between the contacting bodies. .. image:: images/GranularContactModel.png :scale: 60 % :align: center **Figure 2:** Illustration of a contact between two bodies :math:`a` and :math:`b` giving rise to a normal force :math:`f_n`, tangential force :math:`f_t` and rolling resistance :math:`\tau_r`. The contact overlap is denoted by :math:`g`, the unit normal :math:`n` and tangents :math:`t_1` and :math:`t_2`. | ``agx::Physics::GranularBodySystem`` is used to simulate the movement of dynamic soil in :ref:`agxTerrain `. ^^^^^^^^^^^^^^^^^^ Feature Overview ^^^^^^^^^^^^^^^^^^ The ``GranularBodySystem`` contains memory-aligned :ref:`attribute storages ` (position, velocity, inertia, etc... ) for increased efficiency and performance. Granular bodies have two-way interactions with all the other physical systems in AGX Dynamics, such as articulated multibody systems. :ref:`GranularBody particles ` can be created from the system object or by using :ref:`particle emitter `. :ref:`Particle distribution tables ` containing particle models with different radii and :ref:`material properties ` can be used in the creation of particles. The large contact system that often results from granular body interaction can be solved using a :ref:`parallel solver (PPGS) ` for increased performance where a warm starting option exist for increased solution fidelity. Particle states can also be written and read from files via the :ref:`GranularReaderWriter ` class. The granular API is accessible via :ref:`Python `. ------------- Basic Usage ------------- ^^^^^^^^^^^^^^^^^^^ Setup and Creation ^^^^^^^^^^^^^^^^^^^ .. _particle_creation: The following code snippet demonstrates how to create a ``agx::Physics::GranularBodySystem`` system in C++: .. code-block:: c++ // Create and configure a granular body system agx::Physics::GranularBodySystemRef granularBodySystem = new agx::Physics::GranularBodySystem(); simulation->add( granularBodySystem ); Default values for created bodies can be set in the ``GranularBodySystem`` object: .. code-block:: c++ // Properties set on the granularBodySystem are the default values of // all created particles. These can of course be changed on an individual // particle basis after creation. // Most of the methods have optional arguments for updating existing particles bool updateExistingParticles = false; // Configure default values for created particles. // When the radius is updated on particles, the mass is automatically calculated from // the volume change. granularBodySystem->setParticleRadius( 0.1, updateExistingParticles ); // Set the default agx::Material that governs body density and contact material. // When the material is updated on particles, the mass is automatically calculated from // the density in material->getBulkMaterial()->getDensity() granularBodySystem->setMaterial( particleMaterial, updateExistingParticles ); Individial bodies can be created via the ``GranularBodySystem`` object. The body initially has default values but they can be changed after creation: .. code-block:: c++ // Creation returns a agx::Physics::GranularBodyPtr representing the created particle agx::Physics::GranularBodyPtr body = granularBodySystem->createParticle(); // Set position of the particle body->setPosition( Vec3(0, 0, 4.5) ); // Set radius of particle ( updates mass ) body->setRadius( 0.1 ); // Set lifetime of particle in seconds body->setLife( 2.0 ); // Set agx::Material of particle ( updates mass ) body->setMaterial( particleMaterial->getEntity() ); We can also create batches of bodies using bounds specified by the user: .. code-block:: c++ // Create a bunch of particles in a bound. agx::Bound3 spawnBound = agx::Bound3( Vec3(0.0, 0.0, 0), agx::Vec3(4.0, 1.0, 4.0) ); // The radius of all the created particles agx::Real particleRadius = 0.1; // Spacing between the created particles in the bound. // Note that a spacing of 2 x radius will give a tight packing as // particles will be placed next to eachother agx::Vec3 spacing = agx::Vec3( radius * 2, radius * 2, radius * 2 ); // Add a random jitter on the placement of the particles. The jitterFactor is defined as // the fraction of the specified particle spacing which will be used as the bound for // the random perturbation of the particle placement. agx::Real jitterFactor = 0.1; // Place particles in an evenly spaced lattice 3D square lattice ( with jiter ) granularBodySystem->spawnParticlesInBound( bound, radius, spacing, jitterFactor ); // Place particles in a Hexagonal-Close-Packing (HCP) lattice granularBodySystem->spawnParticlesInBoundHCP( bound, radius, spacing, jitterFactor ); Unless specified in the function call, the particles spawned carry default properties set in the ``agx:::Physics::GranularBodySystem`` object. Deleting a granular body can be done in multiple ways, mainly via utiliy functions on the ``GranularBodySystem`` object: .. code-block:: c++ // This will tag the particle for removal at the end of the timestep // ensuring that it will be destroyed in a safe way. // This can be called saftely any time during the timestep and is the // recommended way of removing particles. granularBodySystem->removeParticle( particle ) // You can check if a particle has been tagged for removed by checking it's state variable // Option 1: via the particle instance particle.state().removed() // Option 2: via the granular system ( recommended ). Also checks if the particle instance // has been destroyed ( See below ) granularBodySystem->isValid( particle ) // This function immmediately destroys the particle instance. // NOT recommended to do this in the middle of a simulation step since it // will change the internal structure that can invalidate the simulation. // ONLY do this after the simulation step in post() or last() particle->destroy(); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Motion Control ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ GranularBodies have a ``MotionControl`` state similar to :ref:`RigidBodies `. This can be set on bodies via a utility function on the ``GranularBodySystem`` object. .. code-block:: c++ // Granular bodies support two MotionControls: // 1. DYNAMICS - "Affected by forces, position/velocity is integrated by the *DynamicsSystem*." // 2. KINEMATICS - "Not affected by forces, considered to have infinite mass. Velocities can be set by user." // Make particle kinematic. granularBodySystem->setMotionControl( particle, agx::Physics::GranularBodySystem::MotionControl::KINEMATICS ) // Particle will now effectivley have infinite mass. We now assign velocities to it. particle->setVelocity( agx::Vec3( 1, 0, 0 ) ) // Set fixed translational speed in x-direction particle->setAngularVelocity( agx::Vec3( 0, 0, agx::PI_2 ) ) // Rotational speed of 90 deg / s around z-axis // Returns the motion control of the particle. agx::Physics::GranularBodySystem::MotionControl particleMotionControl = granularBodySystem->getMotionControl( particle ) .. note:: In granular contacts where both granular bodies have the KINEMATICS motion control, the contact forces will not be calculated. This can be used to "freeze" areas of the simulation where particles are in rest in order to increase performance. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Access Granular Data ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. _granular_data_access: Regardless of how the granular bodies were created, individual particles can be accessed and manipulated using ``agx::Physics::GranularBodyPtrs``; the same convention as with ``agx::ParticleSystem``. A ``agx::Physics::GranularBodyPtr`` is an accessor object that gives access to the individual attributes of a granular body, stored as elements in the buffers that make up the ``EntityStorage`` that the particle in question is part of. The ``GranularBodyPtr`` exposes one method of access for each attribute of the GranularBody entities. .. csv-table:: Accessing Granular data :widths: auto "Position", "``getPosition/setPosition``", "Gets/Sets the position (``agx::Vec3``) of the body." "Rotation", "``getRotation/setRotation``", "Gets/Sets the rotation (``agx::Quat``) of the body." "Velocity", "``getVelocity/setVelocity``", "Gets/Sets the velocity (``agx::Vec3``) of the body." "Angular Velocity", "``getAngularVelocity/setAngularVelocity``", "Gets/Sets the angular velocity (``agx::Vec3``) of the body." "Force", "``getForce/setForce``", "Gets/Sets the force (``agx::Vec3``) acting on the body. Useful for using explicit force models on the bodies." "Mass", "``getMass/setMass``", "Gets/Sets the mass of the body. See note about automatic calculation below!" "Inertia", "``getInertia/setInertia``", "Gets/Sets the inertia of the body. See note about automatic calculation below!" "Life", "``getLife/setLife``", "Gets/Sets the lifetime of the particle in seconds." "Color", "``getColor/setColor``", "Gets/Sets the ``agxRender::Color (agx::Vec4f)`` of particle which is used in rendering. " "Material", "``getMaterial/setMaterial``", "Gets/Sets the ``agx::Material`` of the particles which is used for determining contact material and in calculation of mass and density." .. note:: Note that *inertia* and *mass* is automatically calculated when a ``agx::Material`` is assigned via the ``setMaterial`` function to the particle. If you want to manually edit these, be sure to do it after the ``agx::Material`` assignment. .. code-block:: c++ agx::Real averageMass = 0.0; agx::Real averageRadius = 0.0; agx::MaterialRef newMaterial = new agx::Material( "newParticleMaterial" ); newMaterial->getBulkMaterial()->setDensiy( 1700 ); // kg/m3 granularBodies = granularBodySystem->getParticles(); for ( agx::Physics::GranularBodyPtr gPtr : granularBodies ) { gPtr->setVelocity( agx::Vec3( 0, 0, 0 ) ); gPtr->setMaterial( newMaterial ); // This function updates the body mass from BulkMaterial density gPtr->setRadius( 0.05 ); // This function also updates the body mass averageMass += gPtr->getMass(); // extract the calculated mass averageRadius += gPtr->getRadius(); } averageMass /= granularBodySystem->getNumParticles(); averageRadius /= granularBodySystem->getNumParticles(); -------------------------------------------- Material Configuration and Parametrization -------------------------------------------- .. _material_parametrization: Physical properties and dynamic bulk behavior of the granular bodies are governed by the material settings and the contact material in the granular contacts. This section will cover the most important aspects of material parametrization to affect granular bulk behavior. Particle density/mass is primary derived from the set ``agx::Material`` of the particle but can also be overridden manually: .. code-block:: c++ // Set material density agx::MaterialRef particleMaterial = new agx::Material( "particleMaterial" ); particleMaterial->getBulkMaterial()->setDensiy( 1700 ); // kg/m3 // Set the material as the default of granular bodies created from this system // ( unless specified otherwise, e.g from a distribution table in an emitter ) granularBodysystem->setMaterial( particleMaterial ); // Manually create a particle and apply material to it granularBodyPtr = granularBodysystem->createParticle(); granularBodyPtr->setMaterial( particleMaterial ); // Can also override this and apply a fixed mass. Not really recommended since // material and radius assginment overries manually set mass values. granularBodyPtr->setMass( 10 ); // kg // Radius of the particle determines volume and thus resulting mass granularBodyPtr->setRadius( 0.05 ); // m In AGX Dynamics the constitutive equations in granular contacts follows Hertz contact law with Coulomb friction and rolling resistance [#f2]_. Contact material properties determines coefficients used in the contact laws where there are some unique parameters that are only active in granular contacts. The example below shows how to setup a contact material between granular bodies: .. code-block:: c++ // Create new material that will be set to particle agx::MaterialRef particleMaterial = new agx::Material( "particleMaterial" ); simulation->add( particleMaterial ); auto particleParticleCM = simulation->getMaterialManager()->getOrCreateContactMaterial( particleMaterial, particleMaterial ); // Apply the particle material to the granular body system so that it becomes the // default material for new granular bodies granularBodySystem->setMaterial( particleMaterial ); // Set the Young's Modulus of the contact material. In granular contacts, this corresponds // direcly to Young's Modulus in Hertz contact law. particleParticleCM->setYoungsModulus( 1e8 ); // The coulomb friction coefficient in the granular contacts particleParticleCM->setFrictionCoefficient( 0.3 ); // The rolling resistance coefficients in the tangential direction of a contact particleParticleCM->setRollingResistanceCoefficient( 0.3 ); // The twisting resistance coefficient in the normal direction of a contact particleParticleCM->setTwistingResistanceCoefficient(0.001); // The adhesive effects of the granular contacts. Determined by a max force in the adhesive // contact as well as the depth overlap where the adhesive force is active // Arguments are: max adhesive force ( N ), adhesive overlap ( m ) particleParticleCM->setAdhesion( 10, 0.05 ); // Set the restitution of the granular contact particleParticleCM->setRestitution( 0.0 ); // Set the numerical damping of the contact // In high timestep, it is usually some coefficient multiplied by the time step particleParticleCM->setDamping( timeStep * 2 ); .. note:: The rolling and twisting resistance coefficients of the agx::ContactMaterial is unique to the ``granular <-> granular`` and ``granular <-> geometry`` contacts. They are not used in regular geometry and rigid body contacts. ------------------ Collision Groups ------------------ Collision groups can be set to the granular body system in order to configure collision settings in between particles and other geometries. These collision groups will be applied to all particles upon creation with optional argument for applying all existing instances as well. .. code-block:: c++ // Optional argument for updating existing particles with the group bool updateExistingParticles = false; // Add collision group to system granularBodySystem->addCollisionGroup( "group1", updateExistingParticles ); // Remove collision group from system granularBodySystem->removeCollisionGroup( "group1", updateExistingParticles ); // Check if system has collision group ( expecting false ) granularBodySystem->hasCollisionGroup( "group1" ); Groups can be set on individual particles using the ``agx::Physics::GranularBodySystem`` interface: .. code-block:: c++ // Add collision group to a single particle agx::Physics::GranularBodyPtr ptr = granularBodySystem->createParticle() granularBodySystem->addCollisionGroupParticle( "group1", ptr ); // Remove collision group from GranularBody granularBodySystem->removeCollisionGroupParticle( "group1", ptr ); // Check if particle has collision group ( expecting false ) granularBodySystem->hasCollisionGroupParticle( "group1", ptr ); You can also disable collision between particles and specific geometries: .. code-block:: c++ // geometry argument is a agxCollide::Geometry granularBodySystem->setEnableCollisions( geometry, false ); ------------------------------ Particle Distribution Tables ------------------------------ .. _particle_distribution_table: Particle distribution tables are used to represent distributions of particle types with different properties such as radii and material. The ``agx::EmitterDistributionTable`` is a base class in the ``agx::Emitter`` class and can be used for particle distributions in :ref:`Particle emitters ` or for rigid body distributions in :ref:`RigidBody emitters `. Distribution tables contains a set of different particle distribution models. A distribution model consists of a ``radius``, an ``agx::Material`` and a ``weight``. The weight is the probability to create a particle of the corresponding model when the table is queried of a random model. The user specifies a ``agx::Emitter::Quantity`` when a distributions table is created. This quantity is the basis of the probability of how particle models should be selected with respect to their weight. If a model has a higher weight relative to the other models in a table it will contribute more to the total amount of created quantity ( number of particles, mass or volume ) from that table. .. csv-table:: Different quantity settings of the distribution table :widths: auto "**QUANTITY_COUNT**", "Counting number of individual particles as the quantity." "**QUANTITY_VOLUME**", "Counting the total volume of emitted particles as the quantity." "**QUANTITY_MASS**", "Counting the total mass of the emitted particles as the quantity." **Example:** The quantity is set to ``QUANTITY_MASS`` and the table contains two models, one that has weight 1 and the other weight 2. Then the probability of particle creation is set so that the second model will contribute two times more *mass* than the first model to the total mass that is created from the table. If ``QUANTITY_COUNT`` is set, then the second model will yield twice as many *discrete elements* as the first one of the total amount of particles created from the table. .. code-block:: c++ // Create a new distribution table with specified quantity agx::Emitter::DistributionTableRef table = new agx::Emitter::DistributionTable( agx::Emitter::QUANTITY_MASS ) // Arguments to a particle model is: radius, agx::Material, weight table->addModel( new agx::ParticleEmitter::DistributionModel( 0.05, pelletmat, 0.3 ) ) table->addModel( new agx::ParticleEmitter::DistributionModel( 0.09, pelletmat, 0.6 ) ) table->addModel( new agx::ParticleEmitter::DistributionModel( 0.15, pelletmat, 0.1 ) ) Random models can be queried from the table using the following call: .. code-block:: c++ // Query the table for a random model DistributionModel* model = distributionTable->getRandomModel(); // Since this will return an abstract model ( particle and rigid body models are subclasses ), // we need to cast it to a particle model agx::ParticleEmitter::DistributionModel* particleModel = model->as< agx::ParticleEmitter::DistributionModel >(); Distribution tables are used both in :ref:`particle emitters ` and when creating batches of particles. ------------------- Particle Emitter ------------------- .. _particleEmitter: .. image:: images/particle_emitter.png :scale: 60 % :align: center The ``agx::ParticleEmitter`` is used to continuously create a stream of particles in a predefined volume defined by a ``agxCollide::Geometry``. The created particles carry the default properties of the coupled granularBodySystem, with the exception of ``color`` and ``life`` which can be set individually per emitter. .. note:: In order to create emitters that emits bodies/particles with different materials, see :ref:`Distribution tables ` .. code-block:: c++ // When a particle emitter is created, it needs to be coupled to a existing ParticleSystem/granularBodySystem. agx::ParticleEmitterRef emitter = agx::ParticleEmitter( granularBodySystem ); // Set the geometry that will server as a spawning volume for the particles emitter->setGeometry( agxCollide::Geometry( agxCollide::Cylinder( 1.6, 10 ) ) ); // Add the emitter to the simulation to connect it's update task to the simulation simulation->add( emitter ); .. note:: The emitter ``agxCollide::Geometry`` acts just like any other geometry. You can add it to a body in order to create a dynamically moving emitter. A ``Quantity`` needs to be assigned to an emitter to determine the nature of the creation, similar in creation of :ref:`particle distribution tables `. This determines if the given creation rate is based on discrete number of particles, mass or volume. .. code-block:: c++ emitter->setQuantity( agx::Emitter::Quantity::COUNT ); The emitter contains a parameter to control the rate of the particle creation which is dependent on current set quantity of the emitter. .. code-block:: c++ emitter->setRate( 15 ); The maximum amount of quantity ( num particles, mass or volume ) that can be emitted can be adjusted: .. code-block:: c++ /** What the set amount means depends on the current Quantity setting of the emitter. QUANTITY_COUNT - max number of objects QUANTITY_VOLUME - max total volume of objects QUANTITY_MASS - max total mass of objects */ emitter->setMaximumEmittedQuantity( 15 ); You can get information about emitted quantity and the number of particles emitted: .. code-block:: c++ /** * Returns the quantity amount of that this emitter has emitted. */ emitter->getEmittedQuantity() const; /** \return the number of bodies that this emitter has emitted */ emitter->getNumEmitted() const; :ref:`Distribution tables ` can be used in emitters by a simple assignment: .. code-block:: c++ emitter->setDistributionTable( particle_table ); It will then query the table for a particle model when creating new elements instead of using the default settings in the coupled ``GranularBodySystem``. Collision Groups of emitted particles can also be adjusted on an emitter basis. These groups are added on top of the groups that has been added via the granular body system object: .. code-block:: c++ emitter->addCollisionGroup( "emitter_group" ); emitter->removeCollisionGroup( "emitter_group" ); bool hasGroup = emitter ->hasCollisionGroup( "emitter_group" ); The ``life`` parameter of particles can also be set on an emitter basis. .. code-block:: c++ emitter->setLife( 10 ); Initial velocity of created particles can also be set via the emitter: .. code-block:: c++ emitter->setVelocity( agx::Vec3( 1, 1, 1 ) ); --------------------------- Granular Bodies in Python --------------------------- .. _granular_python: The functionality around ``GranularBodySystem`` and particle emitters can be used via the Python API. .. note:: Due to the SWIG mapping the subclass namespace structure in C++ is compressed in the resulting Python API. This primary affects the distribution table, particle model and the ``GranularBodySystem`` classes. .. _table-granular-python: .. csv-table:: Table showing the SWIG mapping of the Granular C++ API to Python API :widths: auto :header: "C++ API", "Python API" " ``agx::Physics::GranularBodySystem`` ", " ``agx.GranularBodySystem``" " ``agx::Emitter::DistributionTable`` ", " ``agx.EmitterDistributionTable``" " ``agx::ParticleEmitter::DistributionModel`` ", " ``agx.ParticleEmitterDistributionModel``" The following code snippet demonstrates how to create a simple granular body system in Python and some example applications for creating bodies. .. code-block:: python # Create and configure a particle system. granular_body_system = agx.GranularBodySystem() sim.add( granular_body_system ) granular_body_system.setMaterial( granularMaterial ) .. code-block:: python # Create a single particle. gBodyPtr = granular_body_system.createParticle() gBodyPtr.setPosition( agx.Vec3( 0, 0, 4.5 ) ) .. code-block:: python # Create a bunch of particles in a bound. particle_radius = 0.1 spawn_bound = agx.Bound3( agx.Vec3(0, 0, 0), agx.Vec3(4,1,4) ) granular_body_system.spawnParticlesInBound( spawn_bound, particle_radius, agx.Vec3( particle_radius * 2 ) ) Sample code for setting up a ``DistributionTable``: .. code-block:: python # Create a new distribution table with specified quantity table = agx.EmitterDistributionTable( agx.Emitter.QUANTITY_MASS ) # Arguments to a particle model is: radius, agx::Material, weight table.addModel(agx.ParticleEmitterDistributionModel( 0.05, pelletmat, 0.3 ) ) table.addModel(agx.ParticleEmitterDistributionModel( 0.09, pelletmat, 0.6 ) ) table.addModel(agx.ParticleEmitterDistributionModel( 0.15, pelletmat, 0.1 ) ) Emitters can also be created in Python: .. code-block:: python # init emitter with mass quantity setting emitter = agx.ParticleEmitter( granularBodySystem, agx.Emitter.QUANTITY_MASS ) # Set the previously declared distribution table to the emitter. emitter.setDistributionTable( table ) emitter.setRate( 15 ) # 15 kg/s emitter.setMaximumEmittedQuantity( 100 ) # 100kg # Create spawn geometry for the emitter. The particles will be created inside the bound of the coupled geometry. emitterSize = agx.Vec3( 0.5, 0.5, 0.5 ) emitterGeom = agxCollide.Geometry( agxCollide.Box( emitterSize ) ) emitterGeom.setPosition( g_emitterSpecification[ "position" ] ) # Disable collisions in emitter to remove contacts with the particles created in the volume. emitterGeom.setEnableCollisions( False ) sim.add( emitterGeom ) emitter.setGeometry( emitterGeom ) # Set an initial velocity downward emitter.setVelocity( agx.Vec3( 0 , 0, -1.0 ) ) sim.add( emitter ) -------------------------------- Extract contact data in Python -------------------------------- Contact data involving granular bodies can be extracted via ``agxCollide::Space``. The following snippet showcases how to extract and read: .. warning:: Contact normals and tangents are in ``agx.Vec3f`` format. They need to be converted to ``agx.Vec3`` in order to be used in multiplication with other data in the ``agx.Vec3`` format that is returned from functions like ``getPosition()`` and ``getVelocity()``. .. code-block:: python # Extract particle <-> particle contacts particle_particle_contacts = sim.getSpace().getParticleParticleContacts() for c in particle_particle_contacts: # Get particle ids from the contact id_particle1 = c.getParticleId1() id_particle2 = c.getParticleId2() # Get ParticlePtr/GranularBodyPtr from ids by using the GranularBodySystem p1 = granularBodySystem.getParticle(id_particle1) p2 = granularBodySystem.getParticle(id_particle2) # Get normals and tangents. These are all in agx.Vec3f format # they need to be converted to agx.Vec3 in order to be multiplied with # agx.Vec3 from methods like 'getVelocity' or 'getPosition' normal = c.getNormal() tu = c.getTangentU() tv = c.getTangentV() # Extract contact force. Note that this will be zero if read in the # 'pre-step' phase of the simulation. Forces are available in the # 'post-step' phase of the simulation which is executed after the solver. f = c.getLocalForce() # Get contact point and depth point = c.getPoint() depth = c.getDepth() # Extract particle <-> geometry contacts particle_geometry_contacts = sim.getSpace().getParticleGeometryContacts() for c in particle_geometry_contacts: # Get particle id in the contact id_particle = c.getParticleId() # Get ParticlePtr/GranularBodyPtr from id by using the GranularBodySystem p = granularBodySystem.getParticle(id_particle) # Get the geometry from the contact and then it's rigid body, if it exists id_geometry = c.getGeometryId() geom = sim.getSpace().getGeometr(id_geometry) # Will return None if geometry has no rigid body rigidBody = geom.getRigidBody() # Get normals and tangents. These are all in agx.Vec3f format # they need to be converted to agx.Vec3 in order to be multiplied with # agx.Vec3 from methods like 'getVelocity' or 'getPosition' normal = c.getNormal() tangentU = c.getTangentU() tangentV = c.getTangentV() # Extract contact force. Note that this will be zero if read in the # 'pre-step' phase of the simulation. Forces are available in the # 'post-step' phase of the simulation which is executed after the solver. f = c.getLocalForce() # Get contact point and depth point = c.getPoint() depth = c.getDepth() Existsing Python tutorials goes further in showcasing ``GranularBodySystem`` functionality. ----------------- Solver Settings ----------------- .. _solver_settings: Solver settings affect the performance and fidelity of ``GranularBodySystem``. Every granular contact is solved using an iterative solver configuration; The typical large contact network that results from big granular systems quickly become unfeasible for a direct solver. .. note:: ``GranularBodySystem`` ``particle <-> particle`` and ``particle <-> geometry contacts`` are **ALWAYS** solved using an iterative solver for all the 6 contact equations *( 1 x normal, 2 x friction, 2 x tangential rolling resistance and 1 x twisting resistance )*. Contact material ``Friction model`` settings have no effect on this. A basic setup of the solver can look like this: .. code-block:: c++ // Set maximum amount of threads for usage in the parallel solver // We prefferable want to use a number of threads equal to number of physical // cores present in the machine. So account for potential HT/SMT. agx::setNumThreads( 0 ); int n = agx::getNumThreads() / 2; // Account for hyper-threading agx::setNumThreads( n ); // Apply simulation solver parameters. // Use the parallel PGS solver for best performance using multiple threads. sim->getDynamicsSystem()->getSolver()->setUseParallelPgs( true ) // Enable 32bit float solver for granular bodies to increase performance. sim.getDynamicsSystem()->getSolver()->setUse32bitGranularBodySolver( true ); // Enable granular warm starting for all solver configurations to enable better solutions. sim->getDynamicsSystem()->getSolver()->setUseGranularWarmStarting( true ); // Set number of iterations that should be used in the Parallel Gauss-Seidel Solver. // The higher number the better solution, but will also impact performance. sim->getDynamicsSystem()->getSolver()->setNumPPGSRestingIterations( 50 ); // Set the number of resting iterations in the regular serial iterative solver. // Higher iterations result in better solutions with increaded computational cost. sim->getDynamicsSystem()->getSolver()->setNumRestingIterations( 50 ); ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Parallel Projected Gauss Seidel (PPGS) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ .. _ppgs: The contact system can be solved by a Parallel-Projected Gauss-Seidel Solver (PPGS). The solver first partitions the system into different subdomains and interaction zones that are solved in parallel via multithreading. The iteration count used in the solver can be adjusted according to fidelity needs. Solver settings can be adjusted to increase or decrease simulation fidelity with respect to computational performance. .. code-block:: c++ // Use the parallel PGS solver for best performance using multiple threads sim->getDynamicsSystem()->getSolver()->setUseParallelPgs( true ) // Set number of iterations that should be used in the Parallel Gauss-Seidel Solver. // The higher number the better solution, but will also impact performance. sim->getDynamicsSystem()->getSolver()->setNumPPGSRestingIterations( 50 ) .. note:: PPGS can be used for rigid body contacts as well. When enabled every rigid body contact with a contact material that has the ``agx::IterativeProjectedConeFriction( agx::FrictionModel::ITERATIVE)`` as set ``Friction Model`` will be included in the PPGS solver. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Warm Starting Granular Contacts ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The PPGS solver can be warm-started by using the constraint multiplier lambda ( :math:`\lambda` ) from the solution of the contact constraint equations in the previous time-step. In general, this will enable contacts that persists over time to converge faster and result in a more precise solution of PPGS solver with fewer number of iterations. This option is recommended when simulating slowly evolving systems with contact networks. For rapidly evolving systems, there is little benefit from using warm-starting or even risk that it worsen the convergence. .. code-block:: c++ // Enable granular warm starting for all solver configurations to enable better solutions sim->getDynamicsSystem()->getSolver()->setUseGranularWarmStarting( true ) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 32bit Granular Solver ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ The iterative solvers that are used when solving the granular contact networks always has some error. The order of magnitude of the error is in all cases higher than the precision of 64it floats. Due to this fact there exists an option of only using 32bit floats in the granular solver to increase compute performance by up to 30-40%. .. code-block:: c++ // Enable granular warm starting for all solver configurations to enable better solutions sim->getDynamicsSystem()->getSolver()->setUseGranularWarmStarting( true ) ----------------------- GranularReaderWriter ----------------------- .. _granularReaderWriter: ``GranularBodySystem`` and particle data can be serialized and stored separately apart from the regular serialization functionality of AGX Dynamics using the ``agx::GranularReaderWriter`` class. This can be used to store granular data in ``.agx``, ``.aagx`` and ``.csv``. .. code-block:: c++ // Extension of filename determines the file format that is stored ( .agx, .aagx or .csv ) std::string filename = "particles.agx"; // Vector of GranularBodyPtr:s, can be individually picked. Can also store every body by // extracting them via the getParticles() function agx::Physics::GranularBodyPtrVector bodies = granularBodySystem.getParticles(); // Optional argument transform is applied to the position and velocity // of the bodies before storing to file. // Create simple translation matrix in this case agx::AffineMatrix4x4 transform = agx::AffineMatrix4x4::translate( 0, 1, 2 ); // Write all the existing bodies to an .agx file agx::GranularReaderWriter::writeFile( filename, bodies, granularBodySystem, transform ); We can then load stored bodies from an ``.agx`` file: .. code-block:: c++ // Specify the .agx file to load bodies from std::string filename = "particles.agx"; // Vector of GranularBodyPtr:s where the loaded bodies will be inserted agx::Physics::GranularBodyPtrVector granularBodies; // Optional argument transform is applied to the position and velocity // of the bodies before loading into the simulation. // In this case, we want to invert the previous translation matrix to bring the bodies back to the origin frame. agx::AffineMatrix4x4 transform = agx::AffineMatrix4x4::translate( 0, 1, 2 ).inverse(); // Load the bodies from the specified filename bool success = loadFileAGX( filename, granularBodies, // loaded bodies will be inserted into this vector granularBodySystem, // loaded bodies will be inserted into this granular body system transform ); .. rubric:: References .. [#f1] C. Lacoursière, Regularized, stabilized, variational methods for multibodies, in: D.F. Peter Bunus, C. Führer (Eds.), The 48th Scandinavian Conference on Simulation and Modeling (SIMS 2007), Linköping University Electronic Press 2007, pp. 40–48 .. [#f2] Wang, D. (2015). Accelerated granular matter simulation (PhD dissertation). Umeå. `pdf `_. .. _rigidBodyEmitter: ---------------------- RigidBody emitter ---------------------- An **agx::RigidBodyEmitter** creates new rigid bodies at a certain rate within a specified volume and adds them to the simulation. This is useful when you do not want to create all the rigid bodies at the start of the simulation and/or when the rigid bodies in your simulation should follow some kind of size/shape distribution. The following code snippet demonstrates how to create a rigid body emitter in c++: .. code-block:: c++ agx::RigidBodyEmitterRef emitter = new agx::RigidBodyEmitter(); emitter->setRate(30.0); emitter->setMaximumEmittedQuantity(200); emitter->setQuantity(agx::Emitter::QUANTITY_COUNT); The rate and the maximumEmittedQuantity decides how fast the emitter creates new bodies and what the maximum quantity is. There are three different ways of counting the quantity of the emitter. .. _table-rigidbody-emitter-quantities: .. csv-table:: Different quantity settings of the emitter and distribution table :widths: auto "**QUANTITY_COUNT**", "Counting number of individual rocks as the quantity." "**QUANTITY_VOLUME**", "Counting the total volume of emitted rocks as the quantity." "**QUANTITY_MASS**", "Counting the total mass of the emitted rocks as the quantity." When the emitter creates a new rigid body it randomly chooses a distributionModel from the distributionTable. A distributionModel consists of a rigidBody template and a weight. The weight is the probability to create a clone of the corresponding rigid body template and adding it to the simulation. By adding several different distributionModels to the distributionTable of the rigid body emitter, the same emitter can emit rigid bodies with different geometries, sizes and materials. A short code example on how to create a distributionTable with one distributionModel: .. code-block:: c++ // Create distribution table agx::Emitter::DistributionTableRef distributionTable = new agx::Emitter::DistributionTable(agx::Emitter::Quantity::QUANTITY_MASS); // Create template agx::RigidBodyRef bodyTemplate1 = new agx::RigidBody(); agxCollide::GeometryRef geometryTemplate1 = new agxCollide::Geometry(new agxCollide::Cylinder(0.4, 0.1)); geometryTemplate1->setMaterial(material); bodyTemplate1->add(geometryTemplate1); distributionTable->addModel(new agx::RigidBodyEmitter::DistributionModel(bodyTemplate1, 0.3)); emitter->setDistributionTable(distributionTable); Finally the rigidBodyEmitter needs a volume in which rigid bodies will spawn. It is defined as an **agxCollide::Geometry**: .. code-block:: c++ agxCollide::GeometryRef emitGeometry = new agxCollide::Geometry(new agxCollide::Sphere(3.0)); emitGeometry->setSensor(true); emitter->setGeometry(emitGeometry);