deargui-vpl/ref/virtools/Samples/Behaviors/ParticlesSystem/behaviors src/ClothSystem.cpp

714 lines
20 KiB
C++

#include "CKAll.h"
#include "ClothSystem.h"
#include "ParticleManager.h"
#define for if (true) for
#if defined(WIN32) || (defined(_XBOX) && (_XBOX_VER<200))
using namespace std;
#include <fvec.h>
#endif
ClothSystem::ClothSystem(CKContext* iCtx)
{
m_ParticleManager = (ParticleManager*)iCtx->GetManagerByGuid(PARTICLE_MANAGER_GUID);
m_Mesh = (CKMesh*)iCtx->CreateObject(CKCID_MESH, "ClothMesh");
m_Mesh->ModifyObjectFlags(CK_OBJECT_NOTTOBESAVED,0);
iCtx->GetCurrentLevel()->AddObject(m_Mesh);
m_Referential = NULL;
m_DampingFactor = 0.95f;
m_PointInverseMass = 1.0f;
m_Nodes = NULL;
m_Iterations = 1;
m_Gravity.Set(0.0f, -9.8f, 0.0f);
}
ClothSystem::~ClothSystem()
{
CKContext* ctx = m_Mesh->GetCKContext();
ctx->DestroyObject(m_Mesh);
// Destroy the nodes
VxDeleteAligned(m_Nodes);
}
void
ClothSystem::Update(const float iDeltaTime)
{
//--- cast deltat time
float castedDeltaTime = (iDeltaTime>60.0f)? 60.0f :iDeltaTime;
//--- smooth time
castedDeltaTime = (3*m_oldDeltaTime+1*castedDeltaTime)/4;
m_oldDeltaTime = castedDeltaTime;
while (castedDeltaTime > 0.0f) {
float remainingTime = castedDeltaTime > 30.0f ? 30.0f : castedDeltaTime;
castedDeltaTime -= remainingTime;
//--- Physics integration
_Integrate(remainingTime*0.001f);
// Constraint solving
#ifdef WIN32
CKBOOL simd = GetProcessorFeatures() & PROC_SIMD;
if (simd) {
_SatisfyConstraintsSIMD(remainingTime*0.001f);
}
else
#endif
{
_SatisfyConstraints(remainingTime*0.001f);
}
}
///
// Mesh updating
// position updating
CKDWORD stride;
void* ptBuf = m_Mesh->GetPositionsPtr(&stride);
int vCount = m_Width * m_Height;
VxCopyStructure(vCount, ptBuf, stride, sizeof(VxVector), m_Nodes, sizeof(Node));
if (m_BackMaterial) {
VxCopyStructure(vCount, (BYTE*)ptBuf+vCount*stride, stride, sizeof(VxVector), m_Nodes, sizeof(Node));
}
// Normal updating
m_Mesh->ModifierVertexMove(TRUE,TRUE);
}
void
ClothSystem::Tesselate(const Vx2DVector& iResolution, const Vx2DVector& iSize, const Vx2DVector& iMapping, CKMaterial* iMaterial, CKMaterial* iBackMaterial)
{
m_oldDeltaTime = 15.0f;
m_Width = (int)iResolution.x;
m_Height = (int)iResolution.y;
int sX = (int)iResolution.x;
int sY = (int)iResolution.y;
int vCount = sX * sY;
m_BackMaterial = iBackMaterial;
// Initialize nodes array
VxDeleteAligned(m_Nodes);
m_Nodes = (Node*)VxNewAligned(sizeof(Node)*vCount,16);
{ // Vertices
if (m_BackMaterial)
m_Mesh->SetVertexCount(vCount*2);
else
m_Mesh->SetVertexCount(vCount);
float deltaX = 1.0f / (iResolution.x-1);
float deltaY = 1.0f / (iResolution.y-1);
float startY = -0.5f;
int vind = 0;
for (int i=0; i<sY; ++i, startY += deltaY) {
float startX = -0.5f;
for (int j=0; j<sX; ++j, startX += deltaX, ++vind) {
m_Nodes[vind].position = VxVector(startX*iSize.x,-startY*iSize.y,0.0f);
m_Nodes[vind].oldPosition = m_Nodes[vind].position;
VxVector v(0,0,-1.0f);
m_Mesh->SetVertexNormal(vind, &v);
m_Mesh->SetVertexPosition(vind, &m_Nodes[vind].position);
m_Mesh->SetVertexTextureCoordinates(vind, iMapping.x*(0.5f + startX), iMapping.y*(0.5f + startY));
}
}
if (m_BackMaterial) { // Duplicates the back vertices
float startY = -0.5f;
int vind = 0;
for (int i=0; i<sY; ++i, startY += deltaY) {
float startX = -0.5f;
for (int j=0; j<sX; ++j, startX += deltaX, ++vind) {
VxVector v(0,0,1.0f);
m_Mesh->SetVertexNormal(vind+vCount, &v);
m_Mesh->SetVertexPosition(vind+vCount, &m_Nodes[vind].position);
m_Mesh->SetVertexTextureCoordinates(vind+vCount, iMapping.x*(0.5f + startX), iMapping.y*(0.5f + startY));
}
}
}
}
{ // Face indices
int fCount = (sX-1) * (sY-1) * 2;
if (m_BackMaterial)
m_Mesh->SetFaceCount(fCount*2);
else
m_Mesh->SetFaceCount(fCount);
int find = 0;
for (int i=0; i<sY-1; ++i) {
for (int j=0; j<sX-1; ++j) {
int a = j + i*sX;
int b = j + (i+1)*sX;
m_Mesh->SetFaceVertexIndex(find,a,a+1,b+1);
m_Mesh->SetFaceMaterial(find++,iMaterial);
m_Mesh->SetFaceVertexIndex(find,a,b+1,b);
m_Mesh->SetFaceMaterial(find++,iMaterial);
}
}
if (m_BackMaterial) {
int find = fCount;
for (int i=0; i<sY-1; ++i) {
for (int j=0; j<sX-1; ++j) {
int a = vCount + j + i*sX;
int b = vCount + j + (i+1)*sX;
m_Mesh->SetFaceVertexIndex(find,a,b+1,a+1);
m_Mesh->SetFaceMaterial(find++,m_BackMaterial);
m_Mesh->SetFaceVertexIndex(find,a,b ,b+1);
m_Mesh->SetFaceMaterial(find++,m_BackMaterial);
}
}
}
}
m_Mesh->BuildFaceNormals();
// calculate constraint strut lengths
// (strut : A structural element used to brace or strengthen a framework by resisting longitudinal compression.)
m_StrutLength.x = iSize.x / (iResolution.x-1.0f);
m_StrutLength.y = iSize.y / (iResolution.y-1.0f);
m_StrutLength.z = sqrtf(m_StrutLength.x*m_StrutLength.x +
m_StrutLength.y*m_StrutLength.y);
ClearHooks();
}
void
ClothSystem::ClearHooks()
{
int vCount = m_Width * m_Height;
CK3dEntity* nullHook = NULL;
VxFillStructure(vCount,&(m_Nodes->hook),sizeof(Node),4,&nullHook);
}
void
ClothSystem::SetMaterial(CKMaterial* iFrontMaterial)
{
CKMaterial* oldMaterial = m_Mesh->GetMaterial(0);
if (oldMaterial != iFrontMaterial)
m_Mesh->ReplaceMaterial(oldMaterial, iFrontMaterial);
}
///////////////////////////////////
// PRIVATE PART
///////////////////////////////////
void
ClothSystem::_Integrate(const float iDeltaTime)
{
float timeStep2 = iDeltaTime*iDeltaTime;
// force vector
VxVector force;
///
// Gravity
if (m_Referential)
m_Referential->InverseTransformVector(&force,&m_Gravity);
else
force = m_Gravity;
///
// Winds Integration
{
CKAttributeManager* am = m_Mesh->GetCKContext()->GetAttributeManager();
const XObjectPointerArray& winds = am->GetAttributeListPtr(m_ParticleManager->m_GlobalWindAttribute);
for (CKObject** it = winds.Begin(); it != winds.End(); ++it) {
CK3dEntity* ent = (CK3dEntity*)*it;
if (!ent->IsActiveInCurrentScene()) continue;
// we get the attribute value
CKParameterOut* pout = ent->GetAttributeParameter(m_ParticleManager->m_GlobalWindAttribute);
float globalWindforce;
pout->GetValue(&globalWindforce);
globalWindforce *= 1000.0f;
VxVector temp;
m_Referential->InverseTransformVector(&temp,&ent->GetWorldMatrix()[2]);
force += temp*globalWindforce;
}
}
///
// Local Winds Integration
{
CKAttributeManager* am = m_Mesh->GetCKContext()->GetAttributeManager();
const XObjectPointerArray& winds = am->GetAttributeListPtr(m_ParticleManager->m_LocalWindAttribute);
for (CKObject** it = winds.Begin(); it != winds.End(); ++it) {
CK3dEntity* ent = (CK3dEntity*)*it;
if (!ent->IsActiveInCurrentScene()) continue;
// we get the attribute value
CKParameterOut* pout = ent->GetAttributeParameter(m_ParticleManager->m_LocalWindAttribute);
Vx2DVector v2(1.0f, 1.0f);
if(pout) pout->GetValue(&v2);
float localWindforce=v2.x;
float decayForce=v2.y;
VxVector wind;
m_Referential->InverseTransformVector(&wind,&ent->GetWorldMatrix()[2]);
wind.Normalize();
VxVector windOrigin;
m_Referential->InverseTransform(&windOrigin,&ent->GetWorldMatrix()[3]);
float radius2 = ent->GetRadius();
radius2 *= radius2;
wind *= localWindforce * timeStep2;
VxVector localPos;
VxMatrix mat;
Vx3DMultiplyMatrix(mat,ent->GetInverseWorldMatrix(),m_Referential->GetWorldMatrix());
int vCount = m_Width * m_Height;
for (int i=0; i < vCount; ++i) {
Vx3DMultiplyMatrixVector(&localPos, mat, &m_Nodes[i].position);
if( localPos.z > 0) { // we check if we are in front of the wind emitter
if( (localPos.x*localPos.x + localPos.y*localPos.y) < radius2) {
float att = 1.0f/(1.0f+localPos.z*decayForce);
m_Nodes[i].oldPosition = m_Nodes[i].position;
m_Nodes[i].position += wind * att;
}
}
}
}
}
force *= timeStep2;
// pointers for traversal
XPtrStrided<VxVector> pt(m_Nodes, sizeof(Node));
// temporary vectors
VxVector temp, velocity;
const int vCount = m_Width * m_Height;
// iterate through all control points
for (int i=0; i < vCount; ++i) {
if (m_Nodes[i].hook) { // Hooked point
m_Referential->InverseTransform(&m_Nodes[i].position, &m_Nodes[i].hook->GetWorldMatrix()[3]);
m_Nodes[i].oldPosition = m_Nodes[i].position;
} else { // Free point
// save current control point location
temp = pt[i];
// update control point by the formula
// x += (x - old_x)*dampingFactor + force*timeStep^2
velocity = pt[i] - m_Nodes[i].oldPosition;
velocity *= 1.0f+(1.0f-m_DampingFactor)*iDeltaTime;
//velocity *= m_DampingFactor;
pt[i] += velocity;
pt[i] += force;
// store old control point location
m_Nodes[i].oldPosition = temp;
}
}
}
inline void mulCSDelta( VxVector& delta, const float strutLength2, const float deltaTimeFactor )
{
delta *= (strutLength2 / (SquareMagnitude(delta)+ strutLength2) - 0.5f)*deltaTimeFactor;
}
void
ClothSystem::_SatisfyConstraints(const float iDeltaTime)
{
VxVector delta;
float deltaTimeFactor = iDeltaTime*50.0f;
// TODO : ajouter de la friction (constante + facteur dependant de la penetration)
for(int iter = 0; iter < m_Iterations; ++iter) {
// iterate through control points
int vcount = m_Width * m_Height;
XPtrStrided<VxVector4> pt(&m_Nodes[0].position,sizeof(Node));
XPtrStrided<VxVector4> nextpt(&m_Nodes[0].position,sizeof(Node));
nextpt += m_Width;
// current isFree flag
XPtrStrided<CK3dEntity*> hook(&m_Nodes[0].hook,sizeof(Node));
// isFree flag one row down
XPtrStrided<int> nexthook(&m_Nodes[0].hook,sizeof(Node));
nexthook += m_Width;
int row = 0, col = 0;
CKContext* ctx = m_Mesh->GetCKContext();
CKAttributeManager* am = ctx->GetAttributeManager();
/***************************************************************
Collision Detection
***************************************************************/
if (m_ParticleManager) {
///
// Infinite plane
const XObjectPointerArray& planes = am->GetAttributeListPtr(m_ParticleManager->m_DInfinitePlaneAttribute);
_SatisfyInfinitePlaneDeflectors(planes);
///
// Sphere Deflectors
const XObjectPointerArray& spheres = am->GetAttributeListPtr(m_ParticleManager->m_DSphereAttribute);
_SatisfySphereDeflectors(spheres);
///
// Box Deflectors
const XObjectPointerArray& boxes = am->GetAttributeListPtr(m_ParticleManager->m_DBoxAttribute);
_SatisfyBoxDeflectors(boxes);
///
// Mesh Deflectors
const XObjectPointerArray& meshes = am->GetAttributeListPtr(m_ParticleManager->m_DObjectAttribute);
_SatisfyMeshDeflectors(meshes);
}
VxTimeProfiler tp;
/***************************************************************
Cloth Mesh Constraints with sqrt approximation
The constraint between xi and xj is satisfied by the formula:
delta = [(xj - xi) * strutLength^2 / (|(xj - xi)|^2 + strutLenth^2)] - 1/2
xi -= delta
xj += delta
***************************************************************/
Node* node = m_Nodes;
Node* nextNode = node + m_Width;
for (int y=1; y<m_Height; ++y) {
for (int x=1; x<m_Width; ++x, ++node, ++nextNode) {
delta = *(VxVector*)(node+1) - *(VxVector*)node;
//delta *= (m_StrutLength2.x / (SquareMagnitude(delta)+ m_StrutLength2.x) - 0.5f);
mulCSDelta( delta, m_StrutLength2.x, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!node[1].hook)
*(VxVector*)(node+1) += delta;
delta = *(VxVector*)(nextNode) - *(VxVector*)node;
//delta *= (m_StrutLength2.y / (SquareMagnitude(delta)+ m_StrutLength2.y) - 0.5f);
mulCSDelta( delta, m_StrutLength2.y, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!nextNode[0].hook)
*(VxVector*)(nextNode) += delta;
delta = *(VxVector*)(nextNode+1) - *(VxVector*)node;
//delta *= (m_StrutLength2.z / (SquareMagnitude(delta)+ m_StrutLength2.z) - 0.5f);
mulCSDelta( delta, m_StrutLength2.z, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!nextNode[1].hook)
*(VxVector*)(nextNode+1) += delta;
}
// last column
delta = *(VxVector*)(nextNode) - *(VxVector*)node;
//delta *= (m_StrutLength2.y / (SquareMagnitude(delta)+ m_StrutLength2.y) - 0.5f);
mulCSDelta( delta, m_StrutLength2.y, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!nextNode[0].hook)
*(VxVector*)(nextNode) += delta;
++node;
++nextNode;
}
// last row
for (int x=1; x<m_Width; ++x, ++node, ++nextNode) {
delta = *(VxVector*)(node+1) - *(VxVector*)node;
//delta *= (m_StrutLength2.x / (SquareMagnitude(delta)+ m_StrutLength2.x) - 0.5f);
mulCSDelta( delta, m_StrutLength2.x, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!node[1].hook)
*(VxVector*)(node+1) += delta;
}
/*
float time = tp.Current();
char buffer[32];
sprintf(buffer,"DeltaTime = %f time = %f",m_Mesh->GetCKContext()->GetTimeManager()->GetLastDeltaTime(),time);
CKRenderManager* rm = m_Mesh->GetCKContext()->GetRenderManager();
rm->RegisterText(buffer,0xff00ffaa,40.0f);
*/
}
}
#if defined(WIN32)
inline F32vec4 add_horizontal_ps(F32vec4 &a)
{
return _mm_add_ps(_mm_shuffle_ps(a, a, 0),_mm_add_ps(_mm_shuffle_ps(a, a, 0x55),_mm_shuffle_ps(a, a, 0xaa)));
}
inline void mulCSDeltaSIMD( F32vec4& deltasimd, const F32vec4& strutLength2, const F32vec4& ftemp, const F32vec4& deltaTimeFactorSIMD, const F32vec4& zero5 )
{
//delta *= (strutLength2 / (SquareMagnitude(delta)+ strutLength2) - 0.5f)*deltaTimeFactor;
deltasimd *= ((strutLength2 / (strutLength2 + ftemp)) - zero5);
deltasimd = _mm_mul_ps( deltasimd, deltaTimeFactorSIMD );
}
void
ClothSystem::_SatisfyConstraintsSIMD(const float iDeltaTime)
{
VxVector delta;
float deltaTimeFactor = iDeltaTime*50.0f;
F32vec4 deltaTimeFactorSIMD( deltaTimeFactor );
F32vec4 strutLength2x(m_StrutLength2.x);
F32vec4 strutLength2y(m_StrutLength2.y);
F32vec4 strutLength2z(m_StrutLength2.z);
F32vec4 zero5(0.5f);
F32vec4 deltasimd;
// TODO : ajouter de la friction (constante + facteur dependant de la penetration)
for(int iter = 0; iter < m_Iterations; ++iter) {
// iterate through control points
int vcount = m_Width * m_Height;
XPtrStrided<VxVector4> pt(&m_Nodes[0].position,sizeof(Node));
XPtrStrided<VxVector4> nextpt(&m_Nodes[0].position,sizeof(Node));
nextpt += m_Width;
// current isFree flag
XPtrStrided<CK3dEntity*> hook(&m_Nodes[0].hook,sizeof(Node));
// isFree flag one row down
XPtrStrided<int> nexthook(&m_Nodes[0].hook,sizeof(Node));
nexthook += m_Width;
int row = 0, col = 0;
CKContext* ctx = m_Mesh->GetCKContext();
CKAttributeManager* am = ctx->GetAttributeManager();
/***************************************************************
Collision Detection
***************************************************************/
if (m_ParticleManager) {
///
// Infinite plane
const XObjectPointerArray& planes = am->GetAttributeListPtr(m_ParticleManager->m_DInfinitePlaneAttribute);
_SatisfyInfinitePlaneDeflectors(planes);
///
// Sphere Deflectors
const XObjectPointerArray& spheres = am->GetAttributeListPtr(m_ParticleManager->m_DSphereAttribute);
_SatisfySphereDeflectors(spheres);
///
// Box Deflectors
const XObjectPointerArray& boxes = am->GetAttributeListPtr(m_ParticleManager->m_DBoxAttribute);
_SatisfyBoxDeflectors(boxes);
///
// Mesh Deflectors
const XObjectPointerArray& meshes = am->GetAttributeListPtr(m_ParticleManager->m_DObjectAttribute);
_SatisfyMeshDeflectors(meshes);
}
VxTimeProfiler tp;
/***************************************************************
Cloth Mesh Constraints with sqrt approximation
The constraint between xi and xj is satisfied by the formula:
delta = [(xj - xi) * strutLength^2 / (|(xj - xi)|^2 + strutLenth^2)] - 1/2
xi -= delta
xj += delta
***************************************************************/
Node* node = m_Nodes;
Node* nextNode = node + m_Width;
for (int y=1; y<m_Height; ++y) {
for (int x=1; x<m_Width; ++x, ++node, ++nextNode) {
int allFree = !((int)node[0].hook | (int)node[1].hook | (int)nextNode[0].hook | (int)nextNode[1].hook);
if (allFree) {
// Horizontal
deltasimd = *(F32vec4*)(node+1) - *(F32vec4*)(node);
F32vec4 ftemp = add_horizontal_ps(deltasimd*deltasimd);
//deltasimd *= ((strutLength2x / (strutLength2x + ftemp)) - zero5);
mulCSDeltaSIMD( deltasimd, strutLength2x, ftemp, deltaTimeFactorSIMD, zero5 );
*(F32vec4*)node -= deltasimd;
*(F32vec4*)(node+1) += deltasimd;
// Vertical
deltasimd = *(F32vec4*)(nextNode) - *(F32vec4*)(node);
ftemp = add_horizontal_ps(deltasimd*deltasimd);
//deltasimd *= ((strutLength2y / (strutLength2y + ftemp)) - zero5);
mulCSDeltaSIMD( deltasimd, strutLength2y, ftemp, deltaTimeFactorSIMD, zero5 );
*(F32vec4*)node -= deltasimd;
*(F32vec4*)(nextNode) += deltasimd;
// Diagonal
deltasimd = *(F32vec4*)(nextNode+1) - *(F32vec4*)(node);
ftemp = add_horizontal_ps(deltasimd*deltasimd);
//deltasimd *= ((strutLength2z / (strutLength2z + ftemp)) - zero5);
mulCSDeltaSIMD( deltasimd, strutLength2z, ftemp, deltaTimeFactorSIMD, zero5 );
*(F32vec4*)node -= deltasimd;
*(F32vec4*)(nextNode+1) += deltasimd;
} else {
// Horizontal
deltasimd = *(F32vec4*)(node+1) - *(F32vec4*)(node);
F32vec4 ftemp = add_horizontal_ps(deltasimd*deltasimd);
deltasimd *= ((strutLength2x / (strutLength2x + ftemp)) - zero5);
if (!node[0].hook)
*(F32vec4*)node -= deltasimd;
if (!node[1].hook)
*(F32vec4*)(node+1) += deltasimd;
// Vertical
deltasimd = *(F32vec4*)(nextNode) - *(F32vec4*)(node);
ftemp = add_horizontal_ps(deltasimd*deltasimd);
//deltasimd *= ((strutLength2y / (strutLength2y + ftemp)) - zero5);
mulCSDeltaSIMD( deltasimd, strutLength2y, ftemp, deltaTimeFactorSIMD, zero5 );
if (!node[0].hook)
*(F32vec4*)node -= deltasimd;
if (!nextNode[0].hook)
*(F32vec4*)(nextNode) += deltasimd;
// Diagonal
deltasimd = *(F32vec4*)(nextNode+1) - *(F32vec4*)(node);
ftemp = add_horizontal_ps(deltasimd*deltasimd);
//deltasimd *= ((strutLength2z / (strutLength2z + ftemp)) - zero5);
mulCSDeltaSIMD( deltasimd, strutLength2z, ftemp, deltaTimeFactorSIMD, zero5 );
if (!node[0].hook)
*(F32vec4*)node -= deltasimd;
if (!nextNode[1].hook)
*(F32vec4*)(nextNode+1) += deltasimd;
}
}
// last column
delta = *(VxVector*)(nextNode) - *(VxVector*)node;
//delta *= (m_StrutLength2.y / (SquareMagnitude(delta)+ m_StrutLength2.y) - 0.5f);
mulCSDelta( delta, m_StrutLength2.y, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!nextNode[0].hook)
*(VxVector*)(nextNode) += delta;
++node;
++nextNode;
}
// last row
for (int x=1; x<m_Width; ++x, ++node, ++nextNode) {
delta = *(VxVector*)(node+1) - *(VxVector*)node;
//delta *= (m_StrutLength2.x / (SquareMagnitude(delta)+ m_StrutLength2.x) - 0.5f);
mulCSDelta( delta, m_StrutLength2.x, deltaTimeFactor );
if (!node[0].hook)
*(VxVector*)node -= delta;
if (!node[1].hook)
*(VxVector*)(node+1) += delta;
}
/*
float time = tp.Current();
char buffer[32];
sprintf(buffer,"DeltaTime = %f time = %f",m_Mesh->GetCKContext()->GetTimeManager()->GetLastDeltaTime(),time);
CKRenderManager* rm = m_Mesh->GetCKContext()->GetRenderManager();
rm->RegisterText(buffer,0xff00ffaa,40.0f);
*/
}
}
#endif // WIN32