本文共 11138 字，大约阅读时间需要 37 分钟。
IGL_INLINE void igl::slim_buildA(const Eigen::SparseMatrix &Dx, const Eigen::SparseMatrix &Dy, const Eigen::SparseMatrix &Dz, const Eigen::MatrixXd &W, std::vector > & IJV){ const int dim = (W.cols() == 4) ? 2 : 3; const int f_n = W.rows(); const int v_n = Dx.cols(); // formula (35) in paper if (dim == 2) { IJV.reserve(4 * (Dx.outerSize() + Dy.outerSize())); /*A = [W11*Dx, W12*Dx; W11*Dy, W12*Dy; W21*Dx, W22*Dx; W21*Dy, W22*Dy];*/ for (int k = 0; k < Dx.outerSize(); ++k) { for (Eigen::SparseMatrix ::InnerIterator it(Dx, k); it; ++it) { int dx_r = it.row(); int dx_c = it.col(); double val = it.value(); IJV.push_back(Eigen::Triplet (dx_r, dx_c, val * W(dx_r, 0))); IJV.push_back(Eigen::Triplet (dx_r, v_n + dx_c, val * W(dx_r, 1))); IJV.push_back(Eigen::Triplet (2 * f_n + dx_r, dx_c, val * W(dx_r, 2))); IJV.push_back(Eigen::Triplet (2 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 3))); } } for (int k = 0; k < Dy.outerSize(); ++k) { for (Eigen::SparseMatrix ::InnerIterator it(Dy, k); it; ++it) { int dy_r = it.row(); int dy_c = it.col(); double val = it.value(); IJV.push_back(Eigen::Triplet (f_n + dy_r, dy_c, val * W(dy_r, 0))); IJV.push_back(Eigen::Triplet (f_n + dy_r, v_n + dy_c, val * W(dy_r, 1))); IJV.push_back(Eigen::Triplet (3 * f_n + dy_r, dy_c, val * W(dy_r, 2))); IJV.push_back(Eigen::Triplet (3 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 3))); } } } else { /*A = [ W11*Dx, W12*Dx, W13*Dx; W11*Dy, W12*Dy, W13*Dy; W11*Dz, W12*Dz, W13*Dz; W21*Dx, W22*Dx, W23*Dx; W21*Dy, W22*Dy, W23*Dy; W21*Dz, W22*Dz, W23*Dz; W31*Dx, W32*Dx, W33*Dx; W31*Dy, W32*Dy, W33*Dy; W31*Dz, W32*Dz, W33*Dz;];*/ // A is a 9#F x 3#V matrix. //DX is a #F x #V matrix, Dx.outerSize is #V, it's a column major matrix. //but reserve size should be 9 * (Dx.outerSize()*4 + Dy.outerSize()*4 + Dz.outerSize()*4) //cause each vertex is shared by 4 tetrahedra. IJV.reserve(9 * (Dx.outerSize() + Dy.outerSize() + Dz.outerSize())); for (int k = 0; k < Dx.outerSize(); k++) { for (Eigen::SparseMatrix ::InnerIterator it(Dx, k); it; ++it) { int dx_r = it.row();//local row index in Dx int dx_c = it.col();//local col index in Dx double val = it.value(); IJV.push_back(Eigen::Triplet (dx_r, dx_c, val * W(dx_r, 0))); IJV.push_back(Eigen::Triplet (dx_r, v_n + dx_c, val * W(dx_r, 1))); IJV.push_back(Eigen::Triplet (dx_r, 2 * v_n + dx_c, val * W(dx_r, 2))); IJV.push_back(Eigen::Triplet (3 * f_n + dx_r, dx_c, val * W(dx_r, 3))); IJV.push_back(Eigen::Triplet (3 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 4))); IJV.push_back(Eigen::Triplet (3 * f_n + dx_r, 2 * v_n + dx_c, val * W(dx_r, 5))); IJV.push_back(Eigen::Triplet (6 * f_n + dx_r, dx_c, val * W(dx_r, 6))); IJV.push_back(Eigen::Triplet (6 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 7))); IJV.push_back(Eigen::Triplet (6 * f_n + dx_r, 2 * v_n + dx_c, val * W(dx_r, 8))); } } for (int k = 0; k < Dy.outerSize(); k++) { for (Eigen::SparseMatrix ::InnerIterator it(Dy, k); it; ++it) { int dy_r = it.row(); int dy_c = it.col(); double val = it.value(); IJV.push_back(Eigen::Triplet (f_n + dy_r, dy_c, val * W(dy_r, 0))); IJV.push_back(Eigen::Triplet (f_n + dy_r, v_n + dy_c, val * W(dy_r, 1))); IJV.push_back(Eigen::Triplet (f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 2))); IJV.push_back(Eigen::Triplet (4 * f_n + dy_r, dy_c, val * W(dy_r, 3))); IJV.push_back(Eigen::Triplet (4 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 4))); IJV.push_back(Eigen::Triplet (4 * f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 5))); IJV.push_back(Eigen::Triplet (7 * f_n + dy_r, dy_c, val * W(dy_r, 6))); IJV.push_back(Eigen::Triplet (7 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 7))); IJV.push_back(Eigen::Triplet (7 * f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 8))); } } for (int k = 0; k < Dz.outerSize(); k++) { for (Eigen::SparseMatrix ::InnerIterator it(Dz, k); it; ++it) { int dz_r = it.row(); int dz_c = it.col(); double val = it.value(); IJV.push_back(Eigen::Triplet (2 * f_n + dz_r, dz_c, val * W(dz_r, 0))); IJV.push_back(Eigen::Triplet (2 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 1))); IJV.push_back(Eigen::Triplet (2 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 2))); IJV.push_back(Eigen::Triplet (5 * f_n + dz_r, dz_c, val * W(dz_r, 3))); IJV.push_back(Eigen::Triplet (5 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 4))); IJV.push_back(Eigen::Triplet (5 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 5))); IJV.push_back(Eigen::Triplet (8 * f_n + dz_r, dz_c, val * W(dz_r, 6))); IJV.push_back(Eigen::Triplet (8 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 7))); IJV.push_back(Eigen::Triplet (8 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 8))); } } }}/// Slim Implementation
为什么A长这样？
看公式（35）
∥ W ( J − Λ ) ∥ \|\mathbf W(\mathbf J - \Lambda)\| ∥W(J−Λ)∥ 注意
J i = [ D x i u , D y i u , D z i u , D x i v , D y i v , D z i v , D x i w , D y i w , D z i w ] \mathbf J_i = [Dx_iu, Dy_iu, Dz_iu, Dx_iv, Dy_iv, Dz_iv, Dx_iw, Dy_iw, Dz_iw] Ji=[Dxiu,Dyiu,Dziu,Dxiv,Dyiv,Dziv,Dxiw,Dyiw,Dziw] 其中
D x i Dx_i Dxi为
D x Dx Dx的第i行
D x Dx Dx为#F x #V的gradient矩阵
对于某一个tetrahedra来说
公式(35)主变成
∥ [ W 11 W 12 W 13 W 21 W 22 W 23 W 31 W 32 W 33 ] ( [ D x i u D y i u D z i u D x i v D y i v D z i v D x i w D y i w D z i w ] − [ R 11 R 12 R 13 R 21 R 22 R 23 R 31 R 32 R 33 ] ) ∥ \|\begin{bmatrix} \mathbf W_{11} & \mathbf W_{12} & \mathbf W_{13}\\ \mathbf W_{21} & \mathbf W_{22} & \mathbf W_{23}\\ \mathbf W_{31} & \mathbf W_{32} & \mathbf W_{33}\\\end{bmatrix} (\begin{bmatrix} Dx_iu & Dy_iu & Dz_iu\\ Dx_iv& Dy_iv& Dz_iv\\ Dx_iw & Dy_iw &Dz_iw\\\end{bmatrix}- \begin{bmatrix} \mathbf R_{11} & \mathbf R_{12} & \mathbf R_{13}\\ \mathbf R_{21} & \mathbf R_{22} & \mathbf R_{23}\\ \mathbf R_{31} & \mathbf R_{32} & \mathbf R_{33}\\\end{bmatrix})\| ∥⎣⎡W11W21W31W12W22W32W13W23W33⎦⎤(⎣⎡DxiuDxivDxiwDyiuDyivDyiwDziuDzivDziw⎦⎤−⎣⎡R11R21R31R12R22R32R13R23R33⎦⎤)∥
然后将对应的元素相乘, 提出u,v,w 再将矩阵元素拉成竖起条状
[ W 11 D x i u + W 12 D x i v + W 13 D x i w W 11 D y i u + W 12 D y i v + W 13 D y i w W 11 D z i u + W 12 D z i v + W 13 D z i w W 21 D x i u + W 22 D x i v + W 23 D x i w W 21 D y i u + W 22 D y i v + W 23 D y i w W 21 D z i u + W 22 D z i v + W 23 D z i w W 31 D x i u + W 32 D x i v + W 33 D x i w W 31 D y i u + W 32 D y i v + W 33 D y i w W 31 D z i u + W 32 D z i v + W 33 D z i w ] = [ W 11 D x i W 12 D x i W 13 D x i W 11 D y i W 12 D y i W 13 D y i W 11 D z i W 12 D z i W 13 D z i W 21 D x i W 22 D x i W 23 D x i W 21 D y i W 22 D y i W 23 D y i W 21 D z i W 22 D z i W 23 D z i W 31 D x i W 32 D x i W 33 D x i W 31 D y i W 32 D y i W 33 D y i W 31 D z i W 32 D z i W 33 D z i ] [ u v w ] \begin{bmatrix} \mathbf W_{11} Dx_iu + \mathbf W_{12}Dx_iv + \mathbf W_{13}Dx_iw\\ \mathbf W_{11} Dy_iu + \mathbf W_{12}Dy_iv + \mathbf W_{13}Dy_iw\\ \mathbf W_{11} Dz_iu + \mathbf W_{12}Dz_iv + \mathbf W_{13}Dz_iw\\ \mathbf W_{21} Dx_iu + \mathbf W_{22}Dx_iv + \mathbf W_{23}Dx_iw\\ \mathbf W_{21} Dy_iu + \mathbf W_{22}Dy_iv + \mathbf W_{23}Dy_iw\\ \mathbf W_{21} Dz_iu + \mathbf W_{22}Dz_iv + \mathbf W_{23}Dz_iw\\ \mathbf W_{31} Dx_iu + \mathbf W_{32}Dx_iv + \mathbf W_{33}Dx_iw\\ \mathbf W_{31} Dy_iu + \mathbf W_{32}Dy_iv + \mathbf W_{33}Dy_iw\\ \mathbf W_{31} Dz_iu + \mathbf W_{32}Dz_iv + \mathbf W_{33}Dz_iw\\ \end{bmatrix}= \begin{bmatrix} \mathbf W_{11} Dx_i & \mathbf W_{12}Dx_i & \mathbf W_{13}Dx_i\\ \mathbf W_{11} Dy_i & \mathbf W_{12}Dy_i & \mathbf W_{13}Dy_i\\ \mathbf W_{11} Dz_i & \mathbf W_{12}Dz_i & \mathbf W_{13}Dz_i\\ \mathbf W_{21} Dx_i & \mathbf W_{22}Dx_i & \mathbf W_{23}Dx_i\\ \mathbf W_{21} Dy_i & \mathbf W_{22}Dy_i & \mathbf W_{23}Dy_i\\ \mathbf W_{21} Dz_i & \mathbf W_{22}Dz_i & \mathbf W_{23}Dz_i\\ \mathbf W_{31} Dx_i & \mathbf W_{32}Dx_i & \mathbf W_{33}Dx_i\\ \mathbf W_{31} Dy_i & \mathbf W_{32}Dy_i & \mathbf W_{33}Dy_i\\ \mathbf W_{31} Dz_i & \mathbf W_{32}Dz_i & \mathbf W_{33}Dz_i\\ \end{bmatrix}\begin{bmatrix}u \\ v \\ w\end{bmatrix} ⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡W11Dxiu+W12Dxiv+W13DxiwW11Dyiu+W12Dyiv+W13DyiwW11Dziu+W12Dziv+W13DziwW21Dxiu+W22Dxiv+W23DxiwW21Dyiu+W22Dyiv+W23DyiwW21Dziu+W22Dziv+W23DziwW31Dxiu+W32Dxiv+W33DxiwW31Dyiu+W32Dyiv+W33DyiwW31Dziu+W32Dziv+W33Dziw⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡W11DxiW11DyiW11DziW21DxiW21DyiW21DziW31DxiW31DyiW31DziW12DxiW12DyiW12DziW22DxiW22DyiW22DziW32DxiW32DyiW32DziW13DxiW13DyiW13DziW23DxiW23DyiW23DziW33DxiW33DyiW33Dzi⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤⎣⎡uvw⎦⎤
将右边的矩阵中i沿#F的维数中竖直展开就变成了矩阵A
同理right hand side 就变成了
IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix &A) { Eigen::VectorXd f_rhs(s.dim * s.dim * s.f_n); f_rhs.setZero(); if (s.dim == 2) { /*b = [W11*R11 + W12*R21; (formula (36)) W11*R12 + W12*R22; W21*R11 + W22*R21; W21*R12 + W22*R22];*/ for (int i = 0; i < s.f_n; i++) { f_rhs(i + 0 * s.f_n) = s.W(i, 0) * s.Ri(i, 0) + s.W(i, 1) * s.Ri(i, 1); f_rhs(i + 1 * s.f_n) = s.W(i, 0) * s.Ri(i, 2) + s.W(i, 1) * s.Ri(i, 3); f_rhs(i + 2 * s.f_n) = s.W(i, 2) * s.Ri(i, 0) + s.W(i, 3) * s.Ri(i, 1); f_rhs(i + 3 * s.f_n) = s.W(i, 2) * s.Ri(i, 2) + s.W(i, 3) * s.Ri(i, 3); } } else { /*b = [W11*R11 + W12*R21 + W13*R31; W11*R12 + W12*R22 + W13*R32; W11*R13 + W12*R23 + W13*R33; W21*R11 + W22*R21 + W23*R31; W21*R12 + W22*R22 + W23*R32; W21*R13 + W22*R23 + W23*R33; W31*R11 + W32*R21 + W33*R31; W31*R12 + W32*R22 + W33*R32; W31*R13 + W32*R23 + W33*R33;];*/ for (int i = 0; i < s.f_n; i++)// i is the offset in face. { f_rhs(i + 0 * s.f_n) = s.W(i, 0) * s.Ri(i, 0) + s.W(i, 1) * s.Ri(i, 1) + s.W(i, 2) * s.Ri(i, 2); f_rhs(i + 1 * s.f_n) = s.W(i, 0) * s.Ri(i, 3) + s.W(i, 1) * s.Ri(i, 4) + s.W(i, 2) * s.Ri(i, 5); f_rhs(i + 2 * s.f_n) = s.W(i, 0) * s.Ri(i, 6) + s.W(i, 1) * s.Ri(i, 7) + s.W(i, 2) * s.Ri(i, 8); f_rhs(i + 3 * s.f_n) = s.W(i, 3) * s.Ri(i, 0) + s.W(i, 4) * s.Ri(i, 1) + s.W(i, 5) * s.Ri(i, 2); f_rhs(i + 4 * s.f_n) = s.W(i, 3) * s.Ri(i, 3) + s.W(i, 4) * s.Ri(i, 4) + s.W(i, 5) * s.Ri(i, 5); f_rhs(i + 5 * s.f_n) = s.W(i, 3) * s.Ri(i, 6) + s.W(i, 4) * s.Ri(i, 7) + s.W(i, 5) * s.Ri(i, 8); f_rhs(i + 6 * s.f_n) = s.W(i, 6) * s.Ri(i, 0) + s.W(i, 7) * s.Ri(i, 1) + s.W(i, 8) * s.Ri(i, 2); f_rhs(i + 7 * s.f_n) = s.W(i, 6) * s.Ri(i, 3) + s.W(i, 7) * s.Ri(i, 4) + s.W(i, 8) * s.Ri(i, 5); f_rhs(i + 8 * s.f_n) = s.W(i, 6) * s.Ri(i, 6) + s.W(i, 7) * s.Ri(i, 7) + s.W(i, 8) * s.Ri(i, 8); } } //f_rhs is #9F x 1 // [ u ] // [ v ] // [ w ] Eigen::VectorXd uv_flat(s.dim *s.v_n); for (int i = 0; i < s.dim; i++) for (int j = 0; j < s.v_n; j++) uv_flat(s.v_n * i + j) = s.V_o(j, i); //first align the vertex dimension, then align the x,y,z dimension. //WGL_M is the area of each element s.rhs = (f_rhs.transpose() * s.WGL_M.asDiagonal() * A).transpose() + s.proximal_p * uv_flat; } }}
转载地址：http://yoxdi.baihongyu.com/