流体模拟四:
各向异性(Anisotropic)表面光滑(2)
理论整理清楚了,我们看一下代码实现,我们首先在我们之前的FluidSystem类中添加几个成员函数和对象:
class FluidSystem{
......//之前的代码
void CalAnisotropicKernel();//各项异性核计算
//设置各项异性开关
void setAnisotropic1(){isAnisotropic=1;}
void setAnisotropic0(){isAnisotropic=0;}
//奇异值分解
void RxSVDecomp3(float w[3], float u[9], float v[9], float eps);
//奇异值分解中的相关计算
template<class T>
inline T RX_SIGN2(const T &a, const T &b){ return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a); }
template<class T>
inline T RX_MAX(const T &a, const T &b){ return ((a > b) ? a : b); }
inline float RxPythag(const float a, const float b)
{
float absa = abs(a), absb = abs(b);
return (absa > absb ? absa*(float)sqrt((double)(1.0+(absb/absa)*(absb/absa))) :
(absb == 0.0 ? 0.0 : absb*(float)sqrt((double)(1.0+(absa/absb)*(absa/absb)))));
}
vector<float> m_hEigen; //协方差矩阵奇异值
vector<float> m_hRMatrix; //协方差矩阵的奇异向量矩阵(旋转矩阵)
vector<float> m_hG; //最后的变换矩阵
vector<float> m_hOldPos; //保存粒子的实际位置
bool isAnisotropic; //是否应用各项异性
};
在这里我们只需要添加各项异性计算函数,以及奇异值分解函数。奇异值分解计算方法单独又可以讲好几篇了,不是本文的重点,所以本文就不对其分析了,我们只要知道参数w是计算出来的奇异值(特征值),u为左奇异矩阵(特征向量列向量构成的矩阵),v为右奇异矩阵,eps为允许的最大误差。
我们还定义了5个成员变量,来保存协方差矩阵的奇异值,奇异向量,最后的变换矩阵,实际模拟流体的粒子位置(注意表面重建的粒子位置不同于模拟中的粒子位置)和一个各项异性开关。
我们接下来看一下奇异值分解和各项异性计算的具体实现:
void FluidSystem::RxSVDecomp3(float w[3], float u[9], float v[9], float eps)
{
bool flag;
int i, its, j, jj, k, l, nm;
float anorm, c, f, g, h, s, scale, x, y, z;
float rv1[3];
g = scale = anorm = 0.0;
for(i = 0; i < 3; ++i){
l = i+2;
rv1[i] = scale*g;
g = s = scale = 0.0;
for(k = i; k < 3; ++k) scale += abs(u[k*3+i]);
if(scale != 0.0){
for(k = i; k < 3; ++k){
u[k*3+i] /= scale;
s += u[k*3+i]*u[k*3+i];
}
f = u[i*3+i];
g = -RX_SIGN2(sqrt(s), f);
h = f*g-s;
u[i*3+i] = f-g;
for(j = l-1; j < 3; ++j){
for(s = 0.0, k = i; k < 3; ++k) s += u[k*3+i]*u[k*3+j];
f = s/h;
for(k = i; k < 3; ++k) u[k*3+j] += f*u[k*3+i];
}
for(k = i; k < 3; ++k) u[k*3+i] *= scale;
}
w[i] = scale*g;
g = s = scale = 0.0;
if(i+1 <= 3 && i+1 != 3){
for(k = l-1; k < 3; ++k) scale += abs(u[i*3+k]);
if(scale != 0.0){
for(k = l-1; k < 3; ++k){
u[i*3+k] /= scale;
s += u[i*3+k]*u[i*3+k];
}
f = u[i*3+l-1];
g = -RX_SIGN2(sqrt(s), f);
h = f*g-s;
u[i*3+l-1] = f-g;
for(k = l-1; k < 3; ++k) rv1[k] = u[i*3+k]/h;
for(j = l-1; j < 3; ++j){
for(s = 0.0,k = l-1; k < 3; ++k) s += u[j*3+k]*u[i*3+k];
for(k = l-1; k < 3; ++k) u[j*3+k] += s*rv1[k];
}
for(k = l-1; k < 3; ++k) u[i*3+k] *= scale;
}
}
anorm = RX_MAX(anorm, (abs(w[i])+abs(rv1[i])));
}
for(i = 2; i >= 0; --i){
if(i < 2){
if(g != 0.0){
for(j = l; j < 3; ++j){
v[j*3+i] = (u[i*3+j]/u[i*3+l])/g;
}
for(j = l; j < 3; ++j){
for(s = 0.0, k = l; k < 3; ++k) s += u[i*3+k]*v[k*3+j];
for(k = l; k < 3; ++k) v[k*3+j] += s*v[k*3+i];
}
}
for(j = l; j < 3; ++j) v[i*3+j] = v[j*3+i] = 0.0;
}
v[i*3+i] = 1.0;
g = rv1[i];
l = i;
}
for(i = 2; i >= 0; --i){
l = i+1;
g = w[i];
for(j = l; j < 3; ++j) u[i*3+j] = 0.0;
if(g != 0.0){
g = 1.0/g;
for(j = l; j < 3; ++j){
for(s = 0.0, k = l; k < 3; ++k) s += u[k*3+i]*u[k*3+j];
f = (s/u[i*3+i])*g;
for(k = i; k < 3; ++k) u[k*3+j] += f*u[k*3+i];
}
for(j = i; j < 3; ++j) u[j*3+i] *= g;
}
else{
for(j = i; j < 3; ++j) u[j*3+i] = 0.0;
}
++u[i*3+i];
}
for(k = 2; k >= 0; --k){
for(its = 0; its < 30; ++its){
flag = true;
for(l = k; l >= 0; --l){
nm = l-1;
if(l == 0 || abs(rv1[l]) <= eps*anorm){
flag = false;
break;
}
if(abs(w[nm]) <= eps*anorm) break;
}
if(flag){
c = 0.0;
s = 1.0;
for(i = l; i < k+1; ++i){
f = s*rv1[i];
rv1[i] = c*rv1[i];
if(abs(f) <= eps*anorm) break;
g = w[i];
h = RxPythag(f, g);
w[i] = h;
h = 1.0/h;
c = g*h;
s = -f*h;
for(j = 0; j < 3; ++j){
y = u[j*3+nm];
z = u[j*3+i];
u[j*3+nm] = y*c+z*s;
u[j*3+i] = z*c-y*s;
}
}
}
z = w[k];
if(l == k){
if(z < 0.0){
w[k] = -z;
for(j = 0; j < 3; ++j) v[j*3+k] = -v[j*3+k];
}
break;
}
if(its == 29){
printf("no convergence in 30 svdcmp iterations");
return;
}
x = w[l];
nm = k-1;
y = w[nm];
g = rv1[nm];
h = rv1[k];
f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g = RxPythag(f, 1.0f);
f = ((x-z)*(x+z)+h*((y/(f+RX_SIGN2(g, f)))-h))/x;
c = s = 1.0;
for(j = l; j <= nm; ++j){
i = j+1;
g = rv1[i];
y = w[i];
h = s*g;
g = c*g;
z = RxPythag(f, h);
rv1[j] = z;
c = f/z;
s = h/z;
f = x*c+g*s;
g = g*c-x*s;
h = y*s;
y *= c;
for(jj = 0; jj < 3; ++jj){
x = v[jj*3+j];
z = v[jj*3+i];
v[jj*3+j] = x*c+z*s;
v[jj*3+i] = z*c-x*s;
}
z = RxPythag(f, h);
w[j] = z;
if(z){
z = 1.0/z;
c = f*z;
s = h*z;
}
f = c*g+s*y;
x = c*y-s*g;
for(jj = 0; jj < 3; ++jj){
y = u[jj*3+j];
z = u[jj*3+i];
u[jj*3+j] = y*c+z*s;
u[jj*3+i] = z*c-y*s;
}
}
rv1[l] = 0.0;
rv1[k] = f;
w[k] = x;
}
}
// reorder
int inc = 1;
float sw;
float su[3], sv[3];
do{
inc *= 3;
inc++;
}while(inc <= 3);
do{
inc /= 3;
for(i = inc; i < 3; ++i){
sw = w[i];
for(k = 0; k < 3; ++k) su[k] = u[k*3+i];
for(k = 0; k < 3; ++k) sv[k] = v[k*3+i];
j = i;
while (w[j-inc] < sw){
w[j] = w[j-inc];
for(k = 0; k < 3; ++k) u[k*3+j] = u[k*3+j-inc];
for(k = 0; k < 3; ++k) v[k*3+j] = v[k*3+j-inc];
j -= inc;
if (j < inc) break;
}
w[j] = sw;
for(k = 0; k < 3; ++k) u[k*3+j] = su[k];
for(k = 0; k < 3; ++k) v[k*3+j] = sv[k];
}
}while(inc > 1);
for(k = 0; k < 3; ++k){
s = 0;
for(i = 0; i < 3; ++i) if(u[i*3+k] < 0.) s++;
for(j = 0; j < 3; ++j) if(v[j*3+k] < 0.) s++;
if(s > 3){
for(i = 0; i < 3; ++i) u[i*3+k] = -u[i*3+k];
for(j = 0; j < 3; ++j) v[j*3+k] = -v[j*3+k];
}
}
}
我们主要看一下下面的计算各项异性代码:
void FluidSystem::CalAnisotropicKernel()
{
if(m_pointBuffer.size() == 0) return;
float h0 = m_smoothRadius;
float h = 2.0*h0;
float lambda = 0.9;
// 计算更新位置
m_hPosW.resize(m_pointBuffer.size()*3,0.0);
for(unsigned int i = 0; i < m_pointBuffer.size(); ++i){
glm::vec3 pos0;
pos0[0] = m_pointBuffer.get(i)->pos.x;
pos0[1] = m_pointBuffer.get(i)->pos.y;
pos0[2] = m_pointBuffer.get(i)->pos.z;
// 相邻的粒子
glm::vec3 posw(0.0);
float sumw = 0.0f;
int cell[64];//2*smoothRadius的领域网格
m_gridContainer.findTwoCells(pos0, h/m_unitScale, cell);
for(int j=0;j<64;j++){
if (cell[j]==-1) continue;
int pndx=m_gridContainer.getGridData(cell[j]);
while (pndx!=-1) {
glm::vec3 pos1;
pos1[0] = m_pointBuffer.get(pndx)->pos.x;
pos1[1] = m_pointBuffer.get(pndx)->pos.y;
pos1[2] = m_pointBuffer.get(pndx)->pos.z;
float r = glm::length(pos1-pos0)*m_unitScale;
if(r < h){
float q = 1-r/h;
float wij = q*q*q;
posw += pos1*wij;
sumw += wij;
}
pndx=m_pointBuffer.get(pndx)->next;
}
}
int DIM=3;
m_hPosW[DIM*i+0] = posw[0]/sumw;
m_hPosW[DIM*i+1] = posw[1]/sumw;
m_hPosW[DIM*i+2] = posw[2]/sumw;
m_hPosW[DIM*i+0] = (1-lambda)*m_pointBuffer.get(i)->pos.x+lambda*m_hPosW[DIM*i+0];
m_hPosW[DIM*i+1] = (1-lambda)*m_pointBuffer.get(i)->pos.y+lambda*m_hPosW[DIM*i+1];
m_hPosW[DIM*i+2] = (1-lambda)*m_pointBuffer.get(i)->pos.z+lambda*m_hPosW[DIM*i+2];
}
m_hOldPos.resize(3*m_pointBuffer.size(),0.0);
for(int i = 0; i < m_pointBuffer.size(); ++i){
m_hOldPos[3*i+0]=m_pointBuffer.get(i)->pos.x;
m_hOldPos[3*i+1]=m_pointBuffer.get(i)->pos.y;
m_hOldPos[3*i+2]=m_pointBuffer.get(i)->pos.z;
}
for(int i = 0; i < m_hPosW.size()/3; ++i){
Point *p=m_pointBuffer.get(i);
p->pos=glm::vec3(m_hPosW[3*i],m_hPosW[3*i+1],m_hPosW[3*i+2]);
}
m_gridContainer.insertParticles(&m_pointBuffer);
// 协方差矩阵 计算
for(unsigned int i = 0; i < m_pointBuffer.size(); ++i){
glm::vec3 xi;
xi[0] = m_hPosW[3*i+0];
xi[1] = m_hPosW[3*i+1];
xi[2] = m_hPosW[3*i+2];
glm::vec3 xiw(0.0);
float sumw = 0.0f;
int cell[64];//2*smoothRadius的领域网格
m_gridContainer.findTwoCells(xi, h/m_unitScale, cell);
for(int j=0;j<64;j++){
if (cell[j]==-1) continue;
int pndx=m_gridContainer.getGridData(cell[j]);
while (pndx!=-1) {
glm::vec3 xj;
xj[0] = m_hPosW[3*pndx+0];
xj[1] = m_hPosW[3*pndx+1];
xj[2] = m_hPosW[3*pndx+2];
float r = glm::length(xi-xj)*m_unitScale;
if(r < h){
float q = 1-r/h;
float wij = q*q*q;
xiw += xj*wij;
sumw += wij;
}
pndx=m_pointBuffer.get(pndx)->next;
}
}
xiw /= sumw;
glm::mat3 c(0.0);
int n = 0;
sumw = 0.0f;
for(int j=0;j<64;j++){
if (cell[j]==-1) continue;
int pndx=m_gridContainer.getGridData(cell[j]);
while (pndx!=-1) {
glm::vec3 xj;
xj[0] = m_hPosW[3*pndx+0];
xj[1] = m_hPosW[3*pndx+1];
xj[2] = m_hPosW[3*pndx+2];
float r = glm::length(xi-xj)*m_unitScale;
if(r < h){
float q = 1-r/h;
float wij = q*q*q;
glm::vec3 dxj = xj-xiw;
for(int k = 0; k < 3; ++k){
c[0][k] += wij*dxj[k]*dxj[0];
c[1][k] += wij*dxj[k]*dxj[1];
c[2][k] += wij*dxj[k]*dxj[2];
}
n++;
sumw += wij;
}
pndx=m_pointBuffer.get(pndx)->next;
}
}
c/=sumw;
// 奇异值分解
float w[3], u[9], v[9];
for(int k = 0; k < 3; ++k){
u[k*3+0] = c[0][k];
u[k*3+1] = c[1][k];
u[k*3+2] = c[2][k];
}
RxSVDecomp3(w, u, v, 1.0e-10);
// 特征值Σ
glm::vec3 sigma;
for(int j = 0; j < 3; ++j){
sigma[j] = w[j];
}
// 特征向量(旋转矩阵R)
glm::mat3 R;
for(int j = 0; j < 3; ++j){
for(int k = 0; k < 3; ++k){
R[j][k] = u[k*3+j];
}
}
// 保存的值用于调试
for(int j = 0; j < 3; ++j){
m_hEigen[3*i+j] = sigma[j];
}
for(int j = 0; j < 3; ++j){
for(int k = 0; k < 3; ++k){
m_hRMatrix[9*i+3*j+k] = R[k][j];
}
}
int ne = 25;//m_params.KernelParticles*0.8;
float ks = 1400;
float kn = 0.5;
float kr = 4.0;
if(n > ne){
for(int j = 1; j < 3; ++j){
sigma[j] = RX_MAX(sigma[j], sigma[0]/kr);
}
sigma *= ks;
}
else{
for(int j = 0; j < 3; ++j){
sigma[j] = kn*1.0;
}
}
// 核变形矩阵G.
glm::mat3 G(0.0);
glm::mat3 Sigma(1.0/sigma.x,0,0,0,1.0/sigma.y,0,0,0,1.0/sigma.z);
G=R*Sigma*glm::transpose(R);
double max_diag = -1.0e10;
for(int j = 0; j < 3; ++j){
for(int k = 0; k < 3; ++k){
if(G[k][j] > max_diag) max_diag = G[k][j];
}
}
G /= max_diag;
//G=glm::inverse(G);
for(int j = 0; j < 3; ++j){
for(int k = 0; k < 3; ++k){
m_hG[9*i+3*j+k] = G[k][j];
}
}
}
}
上篇提到我们为了解决不规则粒子放置的问题,我们将扩散平滑步骤应用于核中心的位置。需要获取一个新的表面重建的粒子位置:
因此我们把计算完的放入m_hPosW数组。然后用m_hOldPos数组保存原粒子位置。这是因为仅用于表面重建。之后因为我们需要计算表面重建,所以我们需要这行代码m_gridContainer.insertParticles(&m_pointBuffer);把的信息加入到网格。这里要注意的是,为了获取足够多的样本信息,我们搜索的领域距离是两倍的光滑核半径,即2h。之后便是计算WPCA的协方差矩阵。
我们需要对所有流体粒子计算一个形变矩阵,所以我们会遍历每个粒子。在每次遍历中,我们先根据:
其中wij为:
计算出xiw,然后新建一个对象glm::mat3 c(0.0);利用:
计算其协方差矩阵。
我们然后把c矩阵计算奇异值分解我们便可以获得其中的特征值Σ,以及特征向量构成的矩阵R。由于我们需要的形变矩阵G是:
其中:
根据文章,我们定义了int ne = 25;; float ks = 1400; float kn = 0.5;float kr = 4.0;最后利用m_hG保存我们每个粒子的变形矩阵G。到这里我们我们各项异性变形矩阵G就计算完成了。
我们对之前的tick()函数做如下修改:
void FluidSystem::tick(){
m_gridContainer.insertParticles(&m_pointBuffer);//每帧刷新粒子位置
glm::vec3 tem=m_sphWallBox.min;
_computerPressure();
_computerForce();
_advance();
if(isAnisotropic)
CalAnisotropicKernel();
CalImplicitFieldDevice(m_rexSize, tem, glm::vec3((1.0/8.0)/m_gridContainer.getDelta()), m_field);
clearSuf();//清空表面数据
m_mcMesh.CreateMeshV(m_field, tem, (1.0/8.0)/m_gridContainer.getDelta(), m_rexSize, m_thre, m_vrts, m_nrms, m_face);
if(isAnisotropic){
//返回最初点位置
for(int i = 0; i < m_hOldPos.size()/3; ++i){
Point* p = m_pointBuffer.get(i);
p->pos=glm::vec3(m_hOldPos[3*i],m_hOldPos[3*i+1],m_hOldPos[3*i+2]);
}
}
culAll();
}
我们在MC算法之前判断是否使用各项异性,如果是,则先计算各项异性形变矩阵G,然后在计算完各项异性后,从新将粒子跟新为流体模拟位置。
我们的形变矩阵是在计算颜色场使用的,因此我们还需要修改计算颜色场函数:
double FluidSystem::CalColorField(double x, double y, double z)
{
// MRK:CalColorField
float c = 0.0;
glm::vec3 pos(x, y, z);
if(pos[0] < m_sphWallBox.min[0]) return c;
if(pos[0] > m_sphWallBox.max[0]) return c;
if(pos[1] < m_sphWallBox.min[1]) return c;
if(pos[1] > m_sphWallBox.max[1]) return c;
if(pos[2] < m_sphWallBox.min[2]) return c;
if(pos[2] > m_sphWallBox.max[2]) return c;
float h = m_smoothRadius;
int cell[8];
m_gridContainer.findCells(pos, h/m_unitScale, cell);
if(!isAnisotropic){
// 近傍粒子(各向同性)
for(int i=0;i<8;i++){
if(cell[i] < 0) continue;
int pndx=m_gridContainer.getGridData(cell[i]);
while(pndx!=-1){
Point* p=m_pointBuffer.get(pndx);
float r = glm::distance(pos,p->pos)*m_unitScale;
float q = h*h-r*r;
if(q>0){
c += m_pointMass*m_kernelPoly6*q*q*q;
}
pndx=p->next;
}
}
}else{
// 近傍粒子(各向异性)
for(int i=0;i<8;i++){
if(cell[i] < 0) continue;
int pndx=m_gridContainer.getGridData(cell[i]);
while(pndx!=-1){
glm::mat3 G(0.0);
for(int j=0;j<3;++j){
for(int k=0;k<3;++k){
G[k][j]=m_hG[pndx*9+3*j+k];
}
}
Point* p=m_pointBuffer.get(pndx);
glm::vec3 rij=pos-p->pos;
rij=G*rij;//形变后的椭球体距离
float r = glm::length(rij)*m_unitScale;
float q = h*h-r*r;
if(q>0){
float detG=glm::determinant(G);//计算行列式
c += detG*m_pointMass*m_kernelPoly6*q*q*q;
}
pndx=p->next;
}
}
}
return c;
}
到这里我们所有的代码就结束了。
我们运行观察结果:
各向同性:
各项异性:
利用简单冯氏照明可以更直接的看到各向同性和各项异性之间的差别:
我们发现表面确实有质的变化。但是,速度上真的是慢了一大截- -。到这里我们的流体模拟就结束啦。