流体模拟四：

各向异性(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;}}// reorderint 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];}}}}

上篇提到我们为了解决不规则粒子放置的问题，我们将扩散平滑步骤应用于核中心的位置。需要获取一个新的表面重建的粒子位置 $\overline{x}$ ：

因此我们把计算完的 $\overline{x}$ 放入m_hPosW数组。然后用m_hOldPos数组保存原粒子位置。这是因为 $\overline{x}$ 仅用于表面重建。之后因为我们需要计算表面重建，所以我们需要这行代码m_gridContainer.insertParticles(&m_pointBuffer);把 $\overline{x}$ 的信息加入到网格。这里要注意的是，为了获取足够多的样本信息，我们搜索的领域距离是两倍的光滑核半径，即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:CalColorFieldfloat 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;}

到这里我们所有的代码就结束了。

我们运行观察结果：

各向同性：

各项异性：

利用简单冯氏照明可以更直接的看到各向同性和各项异性之间的差别：

我们发现表面确实有质的变化。但是，速度上真的是慢了一大截- -。到这里我们的流体模拟就结束啦。

SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(2)相关推荐

SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(1)
流体模拟四: 各向异性(Anisotropic)表面光滑(1) 表面的表示与定义我们是以隐式表面来表示流体表面.定义一个足够大的参考域,它完全包含粒子表示的流体所在的空间D,即.流体表面隐含的定义为 ...
SPH（光滑粒子流体动力学）流体模拟实现六：Position Based Fluid（PBF）
SPH(光滑粒子流体动力学)流体模拟实现六:Position Based Fluid(PBF) PBF方法和前篇提到的PCISPH方法类似,都属于迭代矫正法.PCISPH是通过迭代预测压力,通过压力变 ...
SPH（光滑粒子流体动力学）流体模拟实现七：屏幕空间流体渲染（SSF）
SPH(光滑粒子流体动力学)流体模拟实现七:屏幕空间流体渲染(SSF) 之前都是用的Marching Cube重建流体表面的方法.近来为了做对比实验,尝试了屏幕空间流体渲染的方法.发现屏幕空间的方法不 ...
SPH（光滑粒子流体动力学）流体模拟实现：算法总览
流体模拟(一) 流体模拟算法总体流程: 流体现象广泛存在于自然界.日常生活以及工业生产中,对流体的模拟即流体动画, 一直是基于物理的动画以及计算机图形学的重要研究内容.目前, 基于物理模拟的流体动画 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现1
流体模拟(二) SPH算法实现1: 由于我们计算每个粒子的状态时,都需要获得在它光滑核半径内(邻域内)的所有粒子信息.我们如果遍历每个粒子,计算欧式距离的话,那开销就过于庞大了.因此我们可以将我们的空 ...
SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（1）
流体模拟(三) Marching Cube算法(1) 我们在实现流体表面重建时,需要事先在空间中划分网格,我们的流体系统正好已经完成了此项工作.其次利用Marching Cube算法计算出构成表面的 ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（4）-算法实现2
流体模拟(二) SPH算法实现2: 在前面一节我们完成了粒子缓存类,网格类和邻接表类.我们现在可以正式的整合在我们的流体系统类中了. 流体系统类 class FluidSystem{public:Fl ...
SPH（光滑粒子流体动力学）流体模拟实现二：SPH算法（3）-光滑核函数
流体模拟(二) 光滑核函数: sph中涉及的光滑核可以理解为:在一定的光滑核半径内,所受的力受距离权重的影响,距离越近所受影响越大.其表现形式如图所示. 这里我们便可以将流体看成一个个粒子的集合,每一 ...
SPH（光滑粒子流体动力学）流体模拟实现三：Marching Cube算法（2）
流体模拟(三) Marching Cube算法(2) 我们在之前的流体系统类里新加入一些函数和成员,用来引用我们的MC类,便可以获得生成有表面的流体模型了,效果如图: 添加简单的天空盒以及Phong氏 ...

SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(2)

流体模拟四：

各向异性(Anisotropic)表面光滑(2)

SPH（光滑粒子流体动力学）流体模拟实现四：各向异性(Anisotropic)表面光滑(2)相关推荐

最新文章

热门文章