DirectX11:Position Based Fluid
阿新 • • 發佈：2022-05-21
目錄
DirectX11:Position Based Fluid
前言
這是我本科畢業設計專案，使用DirectX11實現一個基於PBD的流體模擬模擬，閱讀本文假設你對以下內容比較熟悉：（由於內容較多，我也分為了幾篇部落格）
| DirectX11 With Windows SDK--00 目錄 |
|---|
| 流體模擬：Smoothed Particle Hydrodynamics |
| 流體模擬：NeighborHood Search |
| DirectX11：GPU基數排序 |
| 流體模擬：Position Based Fluid |
演算法過程
具體過程
領域搜尋（Neighbor Search）
我們採用空間雜湊的方法對粒子所處的空間網格進行劃分，通過計算其空間雜湊值並將其進行排序，得到當前網格的起始與結束地址。（具體實現可參考：流體模擬：NeighborHood Search）
我們可以通過遍歷當前粒子附近的27個網格得到出其的鄰居粒子（即粒子之間距離少於一定距離，本專案設為粒子半徑），最大鄰居粒子數量設為96。
HLSL核心程式碼
#include "PBFSolverCommon.hlsli" //Find Neighbor Paticle [numthreads(THREAD_NUM_X, 1, 1)] void CS( uint3 DTid : SV_DispatchThreadID ) { if (DTid.x >= g_ParticleNums) { return; } //curr particle pos float3 currPos = g_sortedNewPosition[DTid.x]; float3 sceneMin = g_Bounds[0]; int3 currCellPos = floor((currPos - sceneMin) / g_CellSize); //curr particle index //uint currParitcleIndex = g_ParticleIndex[DTid.x]; int neighborCount = 0; int x = 0, y = 0, z = 0; [unroll(3)] for (z = -1; z <= 1; ++z) { [unroll(3)] for (y = -1; y <= 1; ++y) { [unroll(3)] for (x = -1; x <= 1; ++x) { //find 27 cell neighbor particle int3 neighCellPos = currCellPos + int3(x, y, z); if (neighCellPos.x < 0.0f || neighCellPos.y < 0.0f || neighCellPos.z < 0.0f) { continue; } uint cellHash = GetCellHash(neighCellPos); uint neighborCellStart = g_CellStart[cellHash]; uint neighborCellEnd = g_CellEnd[cellHash]; if (neighborCellStart >= neighborCellEnd) { continue; } for (uint index = neighborCellStart; index < neighborCellEnd; ++index) { //get the cell particle pos float3 neighborPartclePos = g_sortedNewPosition[index]; float3 distance = currPos - neighborPartclePos; float distancesq = dot(distance, distance); if (distancesq < g_ParticleRadiusSq) { //contact int contactsIndex = DTid.x * g_MaxNeighborPerParticle + neighborCount; g_Contacts[contactsIndex] = index; neighborCount++; } if (neighborCount == g_MaxNeighborPerParticle) { g_ContactCounts[DTid.x] = neighborCount; return; } } } } } g_ContactCounts[DTid.x] = neighborCount; }
不可壓縮約束和拉格朗日乘子
流體粒子\(i\)的密度根據SPH方法計算可得：\(\rho_{i}=\sum_{j} m_{j} W\left(p_{i}-p_{j}, h\right)\)
拉格朗日乘子的計算公式為：
\[\nabla_{p k} C_{i}=\frac{1}{\rho_{0}} \begin{cases}\Sigma_{j} \nabla_{p k} W\left(p_{i}-p_{j}, h\right) & \text { if } k=i \\ -\nabla_{p k} W\left(p_{i}-p_{j}, h\right) & \text { if } k=j\end{cases} \] \[\lambda_{i}=-\frac{C_{i}\left(p_{1}, \ldots, p_{n}\right)}{\Sigma_{k}\left|\nabla_{p k} C_{i}\right|^{2}+\epsilon} \]其中密度約束\(C_i\)
HLSL核心程式碼：
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
//curr particle pos
float3 currPos = g_sortedNewPosition[DTid.x];
//curr neighbor count
uint neightborCount = g_ContactCounts[DTid.x];
//clac density
float density = 0;
//clac Lagrange multiplier
float3 gradSum_i = float3(0.0f, 0.0f, 0.0f);
float gradSum_j = 0;
uint i = 0;
for (i = 0; i < neightborCount; ++i)
{
//get the cell particle pos
uint neightborParticleIndex = g_Contacts[DTid.x * g_MaxNeighborPerParticle + i];
float3 neighborPartclePos = g_sortedNewPosition[neightborParticleIndex];
//r=p_i-p_j
float3 r = currPos - neighborPartclePos;
density += WPoly6(r, g_sphSmoothLength);
float3 currGrad = WSpikyGrad(r, g_sphSmoothLength);
currGrad *= g_InverseDensity_0;
gradSum_i += currGrad;
if (neightborParticleIndex != DTid.x)
{
gradSum_j += dot(currGrad, currGrad);
}
}
//debug show
g_Density[DTid.x] = density;
float gradSumTotal = gradSum_j + dot(gradSum_i, gradSum_i);
// evaluate density constraint
float constraint = max(density * g_InverseDensity_0 - 1.0f, 0.0f);
float lambda = -constraint / (gradSumTotal + g_LambdaEps);
g_LambdaMultiplier[DTid.x] = lambda;
}
約束投影與拉伸不穩定性
粒子\(i\)經過上述的約束投影后可以計算出對應的位移向量（包括自身的密度約束以及鄰居粒子的密度約束的共同作用）,可由下面公式計算得出：
\[\begin{aligned} \Delta p_{i} &=\lambda_{i} \nabla_{p_{i}} C_{i}+\Sigma_{j} \lambda_{j} \nabla_{p_{j}} C_{i} \\ &=\frac{1}{\rho}_{0} \Sigma_{j} \lambda_{i} \nabla_{p_{i}} W(r, h)+\left(-\frac{1}{\Sigma_{j}} \lambda_{j} \nabla_{p_{j}} W(r, h)\right) \\ &=\frac{1}{\rho_{0}} \Sigma_{j} \lambda_{i} \nabla_{p_{i}} W(r, h)+\frac{1}{\rho_{0}} \Sigma_{j} \lambda_{j} \nabla_{p_{i}} W(r, h) \\ &=\frac{1}{\rho_{0}} \Sigma_{j}\left(\lambda_{i}+\lambda_{j}\right) \nabla_{p_{i}} W(r, h) \end{aligned} \]以及通過新增排斥項\(S_{corr}\)來避免產生粒子聚集的現象：
\[s_{c o r r}=-k\left(\frac{W\left(p_{i}-p_{j}, h\right)}{W(\Delta q, h)}\right)^{n} \]其中，\(\Delta q\)表示到粒子\(i\)的一個固定距離，通常取\(|\Delta q|=0.1h,...,0.3h\),\(h\)即前面提到的光滑核半徑。此外，公式中的\(k\)可以看作表面張力引數，取值\(k=0.1\),而\(n=4\)。最後公式變成了：
\[\Delta p_{i}=\frac{1}{\rho_{0}} \Sigma_{j}\left(\lambda_{i}+\lambda_{j}+s_{\text {corr }}\right) \nabla_{p_{j}} W\left(p_{i}-p_{j}, h\right) \]HLSL核心程式碼：
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
//curr particle pos
float3 currPos = g_sortedNewPosition[DTid.x];
//curr neighbor count
uint neightborCount = g_ContactCounts[DTid.x];
float currLambda = g_LambdaMultiplier[DTid.x];
float poly6Q = WSpiky(g_DeltaQ, g_sphSmoothLength);
float3 deltaPos = float3(0.0f, 0.0f, 0.0f);
uint i = 0;
for (i = 0; i < neightborCount; ++i)
{
//get the cell particle pos
uint neightborParticleIndex = g_Contacts[DTid.x * g_MaxNeighborPerParticle + i];
float neighborLambda = g_LambdaMultiplier[neightborParticleIndex];
//get the cell particle pos
float3 neighborParticlePos = g_sortedNewPosition[neightborParticleIndex];
//r=p_i-p_j
float3 r = currPos - neighborParticlePos;
float poly6 = WSpiky(r, g_sphSmoothLength);
float diffPoly = poly6 / poly6Q;
float scorr = -g_ScorrK * pow(abs(diffPoly), g_ScorrN);
float coff_j = currLambda + neighborLambda + scorr;
float3 currGrad = WSpikyGrad(r, g_sphSmoothLength);
deltaPos += coff_j * currGrad;
}
deltaPos = deltaPos * g_InverseDensity_0;
g_DeltaPosition[DTid.x] = deltaPos;
}
上述描述的平均約束力確保了收斂性，但是在某些情況下，這種平均會過於激進，並且達到解所需的迭代次數增加。因此，我們需要一個全域性使用者引數\(\omega\)來控制逐次超鬆馳法(SOR)的速度。
\[\boldsymbol{\Delta} \mathbf{x}_{i}=\frac{\omega}{n} \sum_{n} \nabla C_{j} \lambda_{j} \]HLSL核心程式碼：
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID)
{
if (DTid.x >= g_ParticleNums)
{
return;
}
uint totalNeighborNum = g_ContactCounts[DTid.x];
float factor = max(g_SOR * totalNeighborNum, 1.0f);
float3 resPos = g_sortedNewPosition[DTid.x] + g_DeltaPosition[DTid.x] * (1 / factor);
g_sortedNewPosition[DTid.x] = resPos;
}
處理碰撞
粒子與邊界碰撞的處理也一直是一個非常關鍵的問題，一般來說，在SPH方法中採用將邊界粒子也視為粒子（即邊界粒子），將邊界
另一種處理方法是採用SDF(Signed Distance Field)的方式變大空間中容器的位置，然後判斷粒子是否出去了這個空間並將其推回去，從而保證粒子在容器之中。因為時間有限以及場景都是簡單的幾何體，所以本專案直接使用瞭解析SDF函式來解決。若場景複雜，存在多面體模型的話，可以考慮先bake出當前場景的中3D SDF Texture，在Shader中進行取樣進行計算。
因為本場景採用平面對粒子的容器範圍進行限制，所以採用了平面的SDF距離場公式：
我們由高中數學知識可得平面的一般式為：\(A(x-x_0)+B(y-y_0)+C(z-z_0)+h=0\)，其中法線為\((A,B,C)\),\((x_0,y_0,z_0)\)為平面上一點，可得到如下的平面的SDF函式。（更多公式可參考大神inigo quilez的部落格）
//sdf plane function
float sdfPlane(float3 p, float3 n, float h)
{
// n must be normalized
return dot(p, n) + h;
}
所以我們處理粒子的碰撞時，首先在約束求解前得到每個粒子接觸的平面（最大的接觸平面數為6）。在每次迭代求解後，對其進行判斷是否出去了這個空間並將其推回去。
處理碰撞過程中，我還增加了一個friction model，具體公式如下：
\[\boldsymbol{\Delta} \mathbf{x}_{\perp}=\left[\left(\mathbf{x}_{i}^{*}-\mathbf{x}_{i}\right)-\left(\mathbf{x}_{j}^{*}-\mathbf{x}_{j}\right)\right] \perp \mathbf{n} \] \[\boldsymbol{\Delta} \mathbf{x}_{i}=\frac{w_{i}}{w_{i}+w_{j}} \begin{cases}\boldsymbol{\Delta} \mathbf{x}_{\perp}, & \left|\boldsymbol{\Delta} \mathbf{x}_{\perp}\right|<\mu_{s} d \\ \boldsymbol{\Delta} \mathbf{x}_{\perp} \cdot \min \left(\frac{\mu_{k} d}{\left|\Delta \mathbf{x}_{\perp}\right|}, 1\right), & \text { otherwise }\end{cases} \]HLSL核心程式碼：
//CollisionPlane_CS.hlsl
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
float3 currPos = g_sortedNewPosition[DTid.x];
int count = 0;
int i = 0;
[unroll]
for (i = 0; i < g_PlaneNums; ++i)
{
float distance = sdfPlane(currPos, g_Plane[i].xyz, g_Plane[i].w) - g_CollisionDistance;
if (distance < g_CollisionThreshold && count < g_MaxCollisionPlanes)
{
int index = DTid.x * g_MaxCollisionPlanes + count;
g_CollisionPlanes[index] = g_Plane[i];
count++;
}
}
g_CollisionCounts[DTid.x] = count;
}
//SolveContact_CS.hlsl
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
float3 currPos = g_sortedNewPosition[DTid.x];
float3 oldPos = g_sortedOldPosition[DTid.x];
int collisionCount = g_CollisionCounts[DTid.x];
int i = 0;
for (i = 0; i < collisionCount; ++i)
{
int index = DTid.x * g_MaxCollisionPlanes + i;
float4 currPlane = g_CollisionPlanes[index];
float distance = sdfPlane(currPos, currPlane.xyz, currPlane.w) - g_CollisionDistance; //d
if (distance < 0.0f)
{
float3 sdfPos = (-distance) * currPlane.xyz;
//friction model
float3 deltaPos = currPos - oldPos;
float deltaX = dot(deltaPos, currPlane.xyz);
float3 deltaDistane = (-deltaX) * currPlane.xyz + deltaPos; //DeltaX
float deltaLength = dot(deltaDistane, deltaDistane);
[flatten]
if (deltaLength <= (g_StaticFriction * distance)) //|deltaX|< u_s*disctance
{
sdfPos -= deltaDistane;
}
else
{
float dynamicFriction = min((-distance) * 0.01f * rsqrt(deltaLength), 1.0f); //
sdfPos -= dynamicFriction * (deltaDistane);
}
currPos += sdfPos;
}
}
g_UpdatedPosition[DTid.x] = currPos;
}
更新速度
這裡的公式非常簡單：
\[\mathbf{v}^{\mathbf{t}+\mathbf{1}}=\frac{\mathbf{p}^{*}-\mathbf{p}^{\mathbf{t}}}{\Delta t} \]HLSL核心程式碼：
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
float3 oldPos = g_sortedOldPosition[DTid.x];
float3 updatePos = g_sortedNewPosition[DTid.x];
float3 newVec = g_InverseDeltaTime * (updatePos - oldPos);
g_UpdatedVelocity[DTid.x] = newVec;
}
渦輪控制和人工粘性
由於數值耗散，PBD方法通常會引入額外的阻尼，導致整個系統的能來損耗，由此會導致本來該有的一些渦流快速消失。PBF通過vorticity confinement由系統重新注入能量：
\[f_{i}^{\text {vorticity }}=\epsilon\left(N \times \omega_{i}\right) \]上述公式中，\(N=\frac{\eta}{|\eta|}, \eta=\nabla|\omega|_{i}\),而流體粒子的旋度\(w_i\)計算公式如下：
\[\omega_{i}=\nabla \times v=\Sigma_{j}\left(v_{j}-v_{i}\right) \times \nabla_{p_{j}} W\left(p_{i}-p_{j}, h\right) \] \[\eta=\sum_{j} \frac{m_{j}}{\rho_{j}}\left\|\omega_{j}\right\| \nabla_{i} W_{i j} \]Vorticity Confinement的基本思路是：通過新增體積力（body force）\(f_{i}^{\text {vorticity }}\)的方式，在旋度粒子（可理解為比周圍粒子旋轉快的粒子，\(\omega_i\)指向粒子\(i\)的旋轉軸）處加速粒子的旋轉運動，通過這種方式來保持系統的旋度。\(\epsilon\)用來控制Vorticity Confinement的強度。
PBF方法採用XSPH的粘度方法直接更新速度，從而產生粘性阻尼。新增人工粘性（Artificial Viscosity）除了可以增加模擬的數值穩定性，還可以消除非物理振盪（nonphysical oscillations）
\[v_{i}^{\text {new }}=v_{i}+c \Sigma_{j}\left(v_{i}-v_{j}\right) \cdot W\left(p_{i}-p_{j}, h\right) \]HLSL核心程式碼：
#include "PBFSolverCommon.hlsli"
//Clac curl
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
//curr neighbor count
uint neightborCount = g_ContactCounts[DTid.x];
//curr particle pos
float3 currPos = g_sortedNewPosition[DTid.x];
float3 currVec = g_UpdatedVelocity[DTid.x];
float3 currOmega = float3(0.0f, 0.0f, 0.0f);
uint i = 0;
for (i = 0; i < neightborCount; ++i)
{
//get the cell particle pos
uint neightborParticleIndex = g_Contacts[DTid.x * g_MaxNeighborPerParticle + i];
//get the cell particle pos
float3 neighborParticlePos = g_sortedNewPosition[neightborParticleIndex];
//r=p_i-p_j
float3 r = currPos - neighborParticlePos;
//v_j-v_i
float3 deltaVelocity = g_UpdatedVelocity[neightborParticleIndex] - currVec;
float3 currGrad = WSpikyGrad(r, g_sphSmoothLength);
//calc omega
float3 omega_j = cross(deltaVelocity, currGrad);
currOmega += omega_j;
}
float curlLength = length(currOmega);
g_Curl[DTid.x] = float4(currOmega.xyz, curlLength);
}
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
float3 currPos = g_sortedNewPosition[DTid.x];
float3 currVec = g_UpdatedVelocity[DTid.x];
float3 oldVec = g_sortedVelocity[DTid.x];
uint counter = g_ContactCounts[DTid.x];
float3 deltaTotalVec = float3(0.0f, 0.0f, 0.0f);
float3 etaTotal = float3(0.0f, 0.0f, 0.0f);
float density;
for (uint i = 0; i < counter; ++i)
{
//get the cell particle pos
uint neightborParticleIndex = g_Contacts[DTid.x * g_MaxNeighborPerParticle + i];
//get the cell particle pos
float3 neighborParticlePos = g_sortedNewPosition[neightborParticleIndex];
//r=p_i-p_j
float3 r = currPos - neighborParticlePos;
//v_j-v_i
float3 deltaVelocity = g_UpdatedVelocity[neightborParticleIndex] - currVec;
float3 currGrad = WSpikyGrad(r, g_sphSmoothLength);
//vorsitory confinement
float neighCurlLength = g_Curl[neightborParticleIndex].w;
etaTotal += currGrad * neighCurlLength;
//XSPH
float3 deltaVec_j = deltaVelocity * WSpiky(r, g_sphSmoothLength);
deltaTotalVec += deltaVec_j;
//density debug
density += WPoly6(r, g_sphSmoothLength);
}
float3 impulse = float3(0.0f, 0.0f, 0.0f);
//vorticity Confinement
if (length(etaTotal) > 0.0f && g_VorticityConfinement > 0.0f && density>0.0f)
{
float epsilon = g_DeltaTime * g_DeltaTime * g_InverseDensity_0 * g_VorticityConfinement;
float3 currCurl = g_Curl[DTid.x].xyz; //r2
float3 N = normalize(etaTotal);
float3 force = cross(N, currCurl);
impulse += epsilon * force;
}
//XSPH
impulse += g_VorticityC * deltaTotalVec;
// solve plane
uint planeCounts = g_CollisionCounts[DTid.x];
uint resCounts = 0;
float3 restitutionVec = float3(0.0f, 0.0f, 0.0f);
for (uint j = 0; j < planeCounts; ++j)
{
float4 plane = g_CollisionPlanes[DTid.x * g_MaxCollisionPlanes + j];
float distance = sdfPlane(currPos, plane.xyz, plane.w) - 1.001f * g_CollisionDistance;
float oldVecD = dot(oldVec, plane.xyz);
if (distance < 0.0f && oldVecD < 0.0f)
{
float currVecD = dot(currVec, plane.xyz);
float restitutionD = oldVecD * g_Restituion + currVecD;
restitutionVec += plane.xyz * (-restitutionD);
resCounts++;
}
}
resCounts = max(resCounts, 1);
restitutionVec /= resCounts;
impulse += restitutionVec;
g_Impulses[DTid.x] = impulse;
}
最終處理
最終我們只需將求解後的最終結果(位置與速度資訊)輸出即可，本專案還對粒子的最大速度進行一定的限制。
HLSL核心程式碼：
#include "PBFSolverCommon.hlsli"
[numthreads(THREAD_NUM_X, 1, 1)]
void CS( uint3 DTid : SV_DispatchThreadID )
{
if (DTid.x >= g_ParticleNums)
{
return;
}
uint prevIndex = g_Particleindex[DTid.x];
g_SolveredPosition[prevIndex] = g_sortedNewPosition[DTid.x];
float3 currVec = g_UpdatedVelocity[DTid.x];
float3 impulse = g_Impulses[DTid.x];
float3 oldVec = g_oldVelocity[prevIndex];
float3 deltaVec = currVec + impulse - oldVec;
float deltaVecLengthsq = dot(deltaVec, deltaVec);
if (deltaVecLengthsq > (g_MaxVeclocityDelta * g_MaxVeclocityDelta))
{
deltaVec = deltaVec * rsqrt(deltaVecLengthsq) * g_MaxVeclocityDelta;
}
float3 finVec = oldVec + deltaVec;
g_SolveredVelocity[prevIndex] = finVec;
}
C++核心部分
因為程式碼太多，這裡這粗略展示演算法核心過程程式碼，具體程式碼可去下面的倉庫地址檢視：
void FluidSystem::TickLogic(ID3D11DeviceContext* deviceContext, PBFSolver::PBFParams params)
{
m_pPBFSolver->SetPBFParams(params);
for (int i = 0; i < params.subStep; ++i)
{
m_pPBFSolver->PredParticlePosition(deviceContext, *m_pPBFSolverEffect,
m_pParticlePosBuffer->GetBuffer(), m_pParticleVecBuffer->GetBuffer());
//NeighborSearch
m_GpuTimer_NeighBorSearch.Start();
m_pNeighborSearch->BeginNeighborSearch(deviceContext, m_pPBFSolver->GetPredPosition(), m_pParticleIndexBuffer->GetBuffer(), params.cellSize);
m_pNeighborSearch->CalcBounds(deviceContext, *m_pNeighborSearchEffect, m_pPBFSolver->GetPredPosition(), m_pParticleIndexBuffer->GetBuffer(), params.cellSize);
m_pNeighborSearch->RadixSort(deviceContext, *m_pNeighborSearchEffect);
m_pNeighborSearch->FindCellStartAndEnd(deviceContext, *m_pNeighborSearchEffect);
m_pNeighborSearch->EndNeighborSearch();
// Constraint iter solver
m_pPBFSolver->BeginConstraint(deviceContext, *m_pPBFSolverEffect, m_pNeighborSearch->GetSortedParticleIndex(),
m_pNeighborSearch->GetSortedCellStart(), m_pNeighborSearch->GetSortedCellEnd(), m_pNeighborSearch->GetBounds());
m_pPBFSolver->SolverConstraint(deviceContext, *m_pPBFSolverEffect);
m_pPBFSolver->EndConstraint(deviceContext, *m_pPBFSolverEffect);
//update data
m_pParticlePosBuffer->UpdataBufferGPU(deviceContext, m_pPBFSolver->GetSolveredPosition());
m_pParticleVecBuffer->UpdataBufferGPU(deviceContext, m_pPBFSolver->GetSolveredVelocity());
}
m_pPBFSolver->CalcAnisotropy(deviceContext, *m_pPBFSolverEffect);
}
void PBFSolver::SolverConstraint(ID3D11DeviceContext* deviceContext, PBFSolverEffect& effect)
{
effect.SetOutPutUAVByName("g_LambdaMultiplier", m_pLambdaMultiplierBuffer->GetUnorderedAccess());
effect.SetOutPutUAVByName("g_Density", m_pDensityBuffer->GetUnorderedAccess());
effect.SetOutPutUAVByName("g_DeltaPosition", m_pDeltaPositionBuffer->GetUnorderedAccess());
effect.SetOutPutUAVByName("g_UpdatedPosition", m_pUpdatedPositionBuffer->GetUnorderedAccess());
for (int i = 0; i < m_PBFParams.maxSolverIterations; ++i)
{
//calc lagrange multiplier
effect.SetCalcLagrangeMultiplierState();
effect.Apply(deviceContext);
deviceContext->Dispatch(m_BlockNums, 1, 1);
//calc Displacement
effect.SetCalcDisplacementState();
effect.Apply(deviceContext);
deviceContext->Dispatch(m_BlockNums, 1, 1);
//add deltapos
effect.SetADDDeltaPositionState();
effect.Apply(deviceContext);
deviceContext->Dispatch(m_BlockNums, 1, 1);
//solver contacts
effect.SetSolverContactState();
effect.Apply(deviceContext);
deviceContext->Dispatch(m_BlockNums, 1, 1);
m_pSortedNewPostionBuffer->UpdataBufferGPU(deviceContext, m_pUpdatedPositionBuffer->GetBuffer(),m_ParticleNums);
}
}