diff options
Diffstat (limited to 'particles/builtin_particle_forces.cpp')
| -rw-r--r-- | particles/builtin_particle_forces.cpp | 544 |
1 files changed, 544 insertions, 0 deletions
diff --git a/particles/builtin_particle_forces.cpp b/particles/builtin_particle_forces.cpp new file mode 100644 index 0000000..1495526 --- /dev/null +++ b/particles/builtin_particle_forces.cpp @@ -0,0 +1,544 @@ +//========= Copyright Valve Corporation, All rights reserved. ============// +// +// Purpose: particle system code +// +//===========================================================================// + +#include "tier0/platform.h" +#include "particles/particles.h" +#include "filesystem.h" +#include "tier2/tier2.h" +#include "tier2/fileutils.h" +#include "tier1/UtlStringMap.h" +#include "tier1/strtools.h" + +#ifdef USE_BLOBULATOR +// TODO: These should be in public by the time the SDK ships +#include "../common/blobulator/Physics/PhysParticle.h" +#include "../common/blobulator/Physics/PhysParticleCache_inl.h" +#include "../common/blobulator/Physics/PhysTiler.h" +#endif + +// memdbgon must be the last include file in a .cpp file!!! +#include "tier0/memdbgon.h" + +class C_OP_RandomForce : public CParticleOperatorInstance +{ + DECLARE_PARTICLE_OPERATOR( C_OP_RandomForce ); + + uint32 GetWrittenAttributes( void ) const + { + return 0; + } + + uint32 GetReadAttributes( void ) const + { + return 0; + } + + + virtual void AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const; + + Vector m_MinForce; + Vector m_MaxForce; +}; + +void C_OP_RandomForce::AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const +{ + FourVectors box_min,box_max; + box_min.DuplicateVector( m_MinForce * flStrength ); + box_max.DuplicateVector( m_MaxForce * flStrength); + box_max -= box_min; + int nContext = GetSIMDRandContext(); + for(int i=0;i<nBlocks;i++) + { + pAccumulatedForces->x = AddSIMD( + pAccumulatedForces->x, AddSIMD( box_min.x, MulSIMD( box_max.x, RandSIMD( nContext) ) ) ); + pAccumulatedForces->y = AddSIMD( + pAccumulatedForces->y, AddSIMD( box_min.y, MulSIMD( box_max.y, RandSIMD( nContext) ) ) ); + pAccumulatedForces->z = AddSIMD( + pAccumulatedForces->z, AddSIMD( box_min.z, MulSIMD( box_max.z, RandSIMD( nContext) ) ) ); + pAccumulatedForces++; + } + ReleaseSIMDRandContext( nContext ); +} + +DEFINE_PARTICLE_OPERATOR( C_OP_RandomForce, "random force", OPERATOR_GENERIC ); + +BEGIN_PARTICLE_OPERATOR_UNPACK( C_OP_RandomForce ) + DMXELEMENT_UNPACK_FIELD( "min force", "0 0 0", Vector, m_MinForce ) + DMXELEMENT_UNPACK_FIELD( "max force", "0 0 0", Vector, m_MaxForce ) +END_PARTICLE_OPERATOR_UNPACK( C_OP_RandomForce ) + +class C_OP_TwistAroundAxis : public CParticleOperatorInstance +{ + DECLARE_PARTICLE_OPERATOR( C_OP_TwistAroundAxis ); + + uint32 GetWrittenAttributes( void ) const + { + return 0; + } + + uint32 GetReadAttributes( void ) const + { + return PARTICLE_ATTRIBUTE_XYZ_MASK; + } + + + virtual void AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const; + + float m_fForceAmount; + Vector m_TwistAxis; + bool m_bLocalSpace; +}; + +void C_OP_TwistAroundAxis::AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const +{ + FourVectors Twist_AxisInWorldSpace; + Twist_AxisInWorldSpace.DuplicateVector( pParticles->TransformAxis( m_TwistAxis, m_bLocalSpace ) ); + Twist_AxisInWorldSpace.VectorNormalize(); + + Vector vecCenter; + pParticles->GetControlPointAtTime( 0, pParticles->m_flCurTime, &vecCenter ); + FourVectors Center; + Center.DuplicateVector( vecCenter ); + size_t nPosStride; + fltx4 ForceScale = ReplicateX4( m_fForceAmount * flStrength ); + const FourVectors *pPos=pParticles->Get4VAttributePtr( PARTICLE_ATTRIBUTE_XYZ, &nPosStride ); + for(int i=0;i<nBlocks;i++) + { + FourVectors ofs=*pPos; + ofs -= Center; + fltx4 bGoodLen = CmpGtSIMD( ofs*ofs, Four_Epsilons ); + ofs.VectorNormalize(); + FourVectors parallel_comp=ofs; + parallel_comp *= ( ofs*Twist_AxisInWorldSpace ); + ofs-=parallel_comp; + bGoodLen = AndSIMD( bGoodLen, CmpGtSIMD( ofs*ofs, Four_Epsilons ) ); + ofs.VectorNormalize(); + FourVectors TangentialForce = ofs ^ Twist_AxisInWorldSpace; + TangentialForce *= ForceScale; + TangentialForce.x = AndSIMD( TangentialForce.x, bGoodLen ); + TangentialForce.y = AndSIMD( TangentialForce.y, bGoodLen ); + TangentialForce.z = AndSIMD( TangentialForce.z, bGoodLen ); + + *(pAccumulatedForces++) += TangentialForce; + pPos += nPosStride; + } + +} + +DEFINE_PARTICLE_OPERATOR( C_OP_TwistAroundAxis, "twist around axis", OPERATOR_GENERIC ); + +BEGIN_PARTICLE_OPERATOR_UNPACK( C_OP_TwistAroundAxis ) + DMXELEMENT_UNPACK_FIELD( "amount of force", "0", float, m_fForceAmount ) + DMXELEMENT_UNPACK_FIELD( "twist axis", "0 0 1", Vector, m_TwistAxis ) + DMXELEMENT_UNPACK_FIELD( "object local space axis 0/1","0", bool, m_bLocalSpace ) +END_PARTICLE_OPERATOR_UNPACK( C_OP_TwistAroundAxis ) + +class C_OP_AttractToControlPoint : public CParticleOperatorInstance +{ + DECLARE_PARTICLE_OPERATOR( C_OP_AttractToControlPoint ); + + uint32 GetWrittenAttributes( void ) const + { + return 0; + } + + uint32 GetReadAttributes( void ) const + { + return 0; + } + + + virtual uint64 GetReadControlPointMask() const + { + return 1ULL << m_nControlPointNumber; + } + + + virtual void AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const; + + float m_fForceAmount; + float m_fFalloffPower; + int m_nControlPointNumber; +}; + +void C_OP_AttractToControlPoint::AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const +{ + int power_frac=-4.0*m_fFalloffPower; // convert to what pow_fixedpoint_exponent_simd wants + fltx4 fForceScale=ReplicateX4( -m_fForceAmount * flStrength ); + + Vector vecCenter; + pParticles->GetControlPointAtTime( m_nControlPointNumber, pParticles->m_flCurTime, &vecCenter ); + FourVectors Center; + Center.DuplicateVector( vecCenter ); + size_t nPosStride; + const FourVectors *pPos=pParticles->Get4VAttributePtr( PARTICLE_ATTRIBUTE_XYZ, &nPosStride ); + + for(int i=0;i<nBlocks;i++) + { + FourVectors ofs=*pPos; + ofs -= Center; + fltx4 len = ofs.length(); + ofs *= MulSIMD( fForceScale, ReciprocalSaturateSIMD( len )); // normalize and scale + ofs *= Pow_FixedPoint_Exponent_SIMD( len, power_frac ); // * 1/pow(dist, exponent) + fltx4 bGood = CmpGtSIMD( len, Four_Epsilons ); + ofs.x = AndSIMD( bGood, ofs.x ); + ofs.y = AndSIMD( bGood, ofs.y ); + ofs.z = AndSIMD( bGood, ofs.z ); + *(pAccumulatedForces++) += ofs; + pPos += nPosStride; + } + +} + +DEFINE_PARTICLE_OPERATOR( C_OP_AttractToControlPoint, "Pull towards control point", OPERATOR_GENERIC ); + +BEGIN_PARTICLE_OPERATOR_UNPACK( C_OP_AttractToControlPoint ) + DMXELEMENT_UNPACK_FIELD( "amount of force", "0", float, m_fForceAmount ) + DMXELEMENT_UNPACK_FIELD( "falloff power", "2", float, m_fFalloffPower ) + DMXELEMENT_UNPACK_FIELD( "control point number", "0", int, m_nControlPointNumber ) +END_PARTICLE_OPERATOR_UNPACK( C_OP_AttractToControlPoint ) + + +#undef USE_BLOBULATOR // TODO (Ilya): Must fix this code + +#ifdef USE_BLOBULATOR + +class C_OP_LennardJonesForce : public CParticleOperatorInstance +{ + DECLARE_PARTICLE_OPERATOR( C_OP_LennardJonesForce ); + + uint32 GetWrittenAttributes( void ) const + { + return 0; + } + + uint32 GetReadAttributes( void ) const + { + return 0; + } + + void InitParams( CParticleSystemDefinition *pDef, CDmxElement *pElement ) + { + //m_pParticleCache = new ParticleCache(m_fInteractionRadius); + m_pPhysTiler = new PhysTiler(m_fInteractionRadius); + } + + virtual void AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const; + + // TODO: Have to destroy PhysTiler in destructor somewhere!!!! + + //ParticleCache* m_pParticleCache; + PhysTiler* m_pPhysTiler; + float m_fInteractionRadius; + float m_fSurfaceTension; + float m_fLennardJonesRepulsion; + float m_fLennardJonesAttraction; + float m_fMaxRepulsion; + float m_fMaxAttraction; + +private: + //virtual void addParticleForce(PhysParticle* a, PhysParticleCacheNode* bcn, float flStrength, float ts) const; + virtual void addParticleForce(PhysParticle* a, PhysParticle* b, float distSq, float flStrength, float ts) const; +}; + + + + +// TODO: I should make sure I don't have divide by zero errors. +// TODO: ts is not used +void C_OP_LennardJonesForce::addParticleForce(PhysParticle* a, PhysParticle* b, float distSq, float flStrength, float ts) const +{ + float d = sqrtf(distSq); + + //======================================================== + // based on equation of force between two molecules which is + // factor * ((distance/bond_length)^-7 - (distance/bond_length)^-13) + + float f; + if(a->group == b->group) // In the same group + { + float p = a->radius * 2.0f / (d+FLT_EPSILON); + float p2 = p * p; + float p4 = p2 * p2; + + + // Surface tension: + + //Notes: + // Can average the neighbor count between the two particles... + // I tried this, and discovered that rather than averaging, I can take maybe take the + // larger of the two neighbor counts, so the attraction between two particles on the surface will be strong, but + // the attraction between a particle inside and a particle on the surface will be weak. I can also try + // taking the min so that the attraction between a particle on the surface and a particle inside the fluid will + // be strong, but the attraction between two particles completely on the inside will be weak. + // + // int symmetric_neighbor_count = min(a->neighbor_count, b->neighbor_count); + // + // Can try having neighbors only cause stronger attraction (no repulsion) + // Can try lower exponents for the LennardJones forces. + + // This is a trick to prevent single particles from floating off... the less neighbors a particle has.. the more it sticks + // This also tends to simulate surface tension + float surface_tension_modifier = ((24.0f * m_fSurfaceTension) / (a->neighbor_count + b->neighbor_count + 0.1f)) + 1.0f; + //float lennard_jones_force = fLennardJones * 2.0f * (p2 - (p4 * p4)); + float lennard_jones_force = m_fLennardJonesAttraction * p2 - m_fLennardJonesRepulsion*p4; + f = surface_tension_modifier * lennard_jones_force; + + // This is some older code: + //f = ((35.0f * LampScene::simulationSurfaceTension) / (a->neighbor_count + 0.1f)) * (p2 - (p4 * p4)); + // used to be 68' + + + //float factor = (b->neighbor_count < 13 && neighbor_count < 13 ? 4.0f : 0.5f); + //f = factor * (p2 - (p2 * p2 * p2 * p2)); + } + else + { + // This was 3.5 ... made 3.0 so particles get closer when they collide + if(d > a->radius * 3.0f) return; + + float p = a->radius * 4.0f / d; + f = -1.0f * p * p; + } + + // These checks are great to have, but are they really necessary? + // It might also be good to have a limit on velocity + + // Attraction is a positive value. + // Repulsion is negative. + if(f < -m_fMaxRepulsion) f = -m_fMaxRepulsion; + if(f > m_fMaxAttraction) f = m_fMaxAttraction; + + Point3D scaledr = (b->center - a->center) * (f/(d+FLT_EPSILON)); // Dividing by d scales distance down to a unit vector + a->force.add(scaledr); + b->force.subtract(scaledr); +} + +void C_OP_LennardJonesForce::AddForces( FourVectors *pAccumulatedForces, + CParticleCollection *pParticles, + int nBlocks, + float flStrength, + void *pContext ) const +{ + int nParticles = pParticles->m_nActiveParticles; // Not sure if this is correct! + + size_t nPosStride; + const FourVectors *pPos=pParticles->Get4VAttributePtr( PARTICLE_ATTRIBUTE_XYZ, &nPosStride ); + + // The +4 is because particles are stored by PET in blocks of 4 + // However, not every block is full. Thus, nParticles may be + // less than nBlocks*4. Could get rid of this if the swizzling/unswizzling + // loop were better written. + static SmartArray<PhysParticle> imp_particles_sa; // This doesn't specify alignment, might have problems with SSE + while(imp_particles_sa.size < nParticles+4) + { + imp_particles_sa.pushAutoSize(PhysParticle()); + } + + /* + size_t nPrevPosStride; + const FourVectors *pPrevPos=pParticles->Get4VAttributePtr( PARTICLE_ATTRIBUTE_PREV_XYZ, &nPrevPosStride ); + */ + + //m_pParticleCache->beginFrame(); + //m_pParticleCache->beginTile(nParticles); + + m_pPhysTiler->beginFrame(Point3D(0.0f, 0.0f, 0.0f)); + + // Unswizzle from the FourVectors format into particles + for(int i=0, p=0;i<nBlocks;i++) + { + FourVectors ofs=*pPos; + + PhysParticle* particle = &(imp_particles_sa[p]); + particle->force.clear(); + if(p < nParticles) + { + particle->center = ofs.Vec(0); + particle->group = 0; + particle->neighbor_count = 0; + m_pPhysTiler->insertParticle(particle); + } + p++; + + particle = &(imp_particles_sa[p]); + particle->force.clear(); + if(p < nParticles) + { + particle->center = ofs.Vec(1); + particle->group = 0; + particle->neighbor_count = 0; + m_pPhysTiler->insertParticle(particle); + } + p++; + + particle = &(imp_particles_sa[p]); + particle->force.clear(); + if(p < nParticles) + { + particle->center = ofs.Vec(2); + particle->group = 0; + particle->neighbor_count = 0; + m_pPhysTiler->insertParticle(particle); + } + p++; + + particle = &(imp_particles_sa[p]); + particle->force.clear(); + if(p < nParticles) + { + particle->center = ofs.Vec(3); + particle->group = 0; + particle->neighbor_count = 0; + m_pPhysTiler->insertParticle(particle); + } + p++; + + pPos += nPosStride; + } + + m_pPhysTiler->processTiles(); + + + float timeStep = 1.0f; // This should be customizable + float nearNeighborInteractionRadius = 2.3f; + float nearNeighborInteractionRadiusSq = nearNeighborInteractionRadius * nearNeighborInteractionRadius; + + PhysParticleCache* pCache = m_pPhysTiler->getParticleCache(); + + // Calculate number of near neighbors for each particle + for(int i = 0; i < nParticles; i++) + { + PhysParticle *b1 = &(imp_particles_sa[i]); + + PhysParticleAndDist* node = pCache->get(b1); + + while(node->particle != NULL) + { + PhysParticle* b2 = node->particle; + + // Compare addresses of the two particles. This makes sure we apply a force only once between a pair of particles. + if(b1 < b2 && node->distSq < nearNeighborInteractionRadiusSq) + { + b1->neighbor_count++; + b2->neighbor_count++; + + } + + node++; + } + } + + // Calculate forces on particles due to other particles + for(int i = 0; i < nParticles; i++) + { + PhysParticle *b1 = &(imp_particles_sa[i]); + + PhysParticleAndDist* node = pCache->get(b1); + + while(node->particle != NULL) + { + PhysParticle* b2 = node->particle; + + // Compare addresses of the two particles. This makes sure we apply a force only once between a pair of particles. + if(b1 < b2) + { + addParticleForce(b1, b2, node->distSq, flStrength, timeStep); + } + + node++; + } + } + + + /* + for(ParticleListNode* bit3 = particles; bit3; bit3 = bit3->next) + { + Particle* b = bit3->particle; + b->prev_group = b->group; // Set prev group + //b1->addDirDragForce(); + b->move(ts); // Move the particle (it should never be used again until next iteration) + } + */ + + m_pPhysTiler->endFrame(); + + // Swizzle forces back into FourVectors format + for(int i=0;i<nBlocks;i++) + { + pAccumulatedForces->X(0) += imp_particles_sa[i*4].force[0]; + pAccumulatedForces->Y(0) += imp_particles_sa[i*4].force[1]; + pAccumulatedForces->Z(0) += imp_particles_sa[i*4].force[2]; + + pAccumulatedForces->X(1) += imp_particles_sa[i*4+1].force[0]; + pAccumulatedForces->Y(1) += imp_particles_sa[i*4+1].force[1]; + pAccumulatedForces->Z(1) += imp_particles_sa[i*4+1].force[2]; + + pAccumulatedForces->X(2) += imp_particles_sa[i*4+2].force[0]; + pAccumulatedForces->Y(2) += imp_particles_sa[i*4+2].force[1]; + pAccumulatedForces->Z(2) += imp_particles_sa[i*4+2].force[2]; + + pAccumulatedForces->X(3) += imp_particles_sa[i*4+3].force[0]; + pAccumulatedForces->Y(3) += imp_particles_sa[i*4+3].force[1]; + pAccumulatedForces->Z(3) += imp_particles_sa[i*4+3].force[2]; + + pAccumulatedForces++; + } + +} + +DEFINE_PARTICLE_OPERATOR( C_OP_LennardJonesForce, "lennard jones force", OPERATOR_GENERIC ); + +BEGIN_PARTICLE_OPERATOR_UNPACK( C_OP_LennardJonesForce ) +DMXELEMENT_UNPACK_FIELD( "interaction radius", "4", float, m_fInteractionRadius ) +DMXELEMENT_UNPACK_FIELD( "surface tension", "1", float, m_fSurfaceTension ) +DMXELEMENT_UNPACK_FIELD( "lennard jones attractive force", "1", float, m_fLennardJonesAttraction ) +DMXELEMENT_UNPACK_FIELD( "lennard jones repulsive force", "1", float, m_fLennardJonesRepulsion ) +DMXELEMENT_UNPACK_FIELD( "max repulsion", "100", float, m_fMaxRepulsion ) +DMXELEMENT_UNPACK_FIELD( "max attraction", "100", float, m_fMaxAttraction ) +END_PARTICLE_OPERATOR_UNPACK( C_OP_LennardJonesForce ) + +#endif + + +void AddBuiltInParticleForceGenerators( void ) +{ + REGISTER_PARTICLE_OPERATOR( FUNCTION_FORCEGENERATOR, C_OP_RandomForce ); + REGISTER_PARTICLE_OPERATOR( FUNCTION_FORCEGENERATOR, C_OP_TwistAroundAxis ); + REGISTER_PARTICLE_OPERATOR( FUNCTION_FORCEGENERATOR, C_OP_AttractToControlPoint ); + #ifdef USE_BLOBULATOR + REGISTER_PARTICLE_OPERATOR( FUNCTION_FORCEGENERATOR, C_OP_LennardJonesForce ); + #endif +} + |