709 SLIM

最新推荐文章于 2020-08-28 09:34:00 发布

seamanj

最新推荐文章于 2020-08-28 09:34:00 发布

阅读量246

点赞数

CC 4.0 BY-SA版权

分类专栏： paper IGL library

本文链接：https://blog.youkuaiyun.com/seamanj/article/details/84039387

paper 同时被 2 个专栏收录

106 篇文章

订阅专栏

IGL library

21 篇文章

订阅专栏


IGL_INLINE void igl::slim_buildA(const Eigen::SparseMatrix<double> &Dx,
          const Eigen::SparseMatrix<double> &Dy,
          const Eigen::SparseMatrix<double> &Dz,
          const Eigen::MatrixXd &W,
    std::vector<Eigen::Triplet<double> > & 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<double>::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<double>(dx_r, dx_c, val * W(dx_r, 0)));
        IJV.push_back(Eigen::Triplet<double>(dx_r, v_n + dx_c, val * W(dx_r, 1)));

        IJV.push_back(Eigen::Triplet<double>(2 * f_n + dx_r, dx_c, val * W(dx_r, 2)));
        IJV.push_back(Eigen::Triplet<double>(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<double>::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<double>(f_n + dy_r, dy_c, val * W(dy_r, 0)));
        IJV.push_back(Eigen::Triplet<double>(f_n + dy_r, v_n + dy_c, val * W(dy_r, 1)));

        IJV.push_back(Eigen::Triplet<double>(3 * f_n + dy_r, dy_c, val * W(dy_r, 2)));
        IJV.push_back(Eigen::Triplet<double>(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<double>::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<double>(dx_r, dx_c, val * W(dx_r, 0)));
        IJV.push_back(Eigen::Triplet<double>(dx_r, v_n + dx_c, val * W(dx_r, 1)));
        IJV.push_back(Eigen::Triplet<double>(dx_r, 2 * v_n + dx_c, val * W(dx_r, 2)));

        IJV.push_back(Eigen::Triplet<double>(3 * f_n + dx_r, dx_c, val * W(dx_r, 3)));
        IJV.push_back(Eigen::Triplet<double>(3 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 4)));
        IJV.push_back(Eigen::Triplet<double>(3 * f_n + dx_r, 2 * v_n + dx_c, val * W(dx_r, 5)));

        IJV.push_back(Eigen::Triplet<double>(6 * f_n + dx_r, dx_c, val * W(dx_r, 6)));
        IJV.push_back(Eigen::Triplet<double>(6 * f_n + dx_r, v_n + dx_c, val * W(dx_r, 7)));
        IJV.push_back(Eigen::Triplet<double>(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<double>::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<double>(f_n + dy_r, dy_c, val * W(dy_r, 0)));
        IJV.push_back(Eigen::Triplet<double>(f_n + dy_r, v_n + dy_c, val * W(dy_r, 1)));
        IJV.push_back(Eigen::Triplet<double>(f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 2)));

        IJV.push_back(Eigen::Triplet<double>(4 * f_n + dy_r, dy_c, val * W(dy_r, 3)));
        IJV.push_back(Eigen::Triplet<double>(4 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 4)));
        IJV.push_back(Eigen::Triplet<double>(4 * f_n + dy_r, 2 * v_n + dy_c, val * W(dy_r, 5)));

        IJV.push_back(Eigen::Triplet<double>(7 * f_n + dy_r, dy_c, val * W(dy_r, 6)));
        IJV.push_back(Eigen::Triplet<double>(7 * f_n + dy_r, v_n + dy_c, val * W(dy_r, 7)));
        IJV.push_back(Eigen::Triplet<double>(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<double>::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<double>(2 * f_n + dz_r, dz_c, val * W(dz_r, 0)));
        IJV.push_back(Eigen::Triplet<double>(2 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 1)));
        IJV.push_back(Eigen::Triplet<double>(2 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 2)));

        IJV.push_back(Eigen::Triplet<double>(5 * f_n + dz_r, dz_c, val * W(dz_r, 3)));
        IJV.push_back(Eigen::Triplet<double>(5 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 4)));
        IJV.push_back(Eigen::Triplet<double>(5 * f_n + dz_r, 2 * v_n + dz_c, val * W(dz_r, 5)));

        IJV.push_back(Eigen::Triplet<double>(8 * f_n + dz_r, dz_c, val * W(dz_r, 6)));
        IJV.push_back(Eigen::Triplet<double>(8 * f_n + dz_r, v_n + dz_c, val * W(dz_r, 7)));
        IJV.push_back(Eigen::Triplet<double>(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)\|$
注意
$Ji=[Dxiu,Dyiu,Dziu,Dxiv,Dyiv,Dziv,Dxiw,Dyiw,Dziw]\mathbf J_i = [Dx_iu, Dy_iu, Dz_iu, Dx_iv, Dy_iv, Dz_iv, Dx_iw, Dy_iw, Dz_iw]$
其中 $Dx_i$ 为 $D x$ 的第i行
$D x$ 为#F x #V的gradient矩阵
对于某一个tetrahedra来说
公式(35)主变成

$∥[W11W12W13W21W22W23W31W32W33]([DxiuDyiuDziuDxivDyivDzivDxiwDyiwDziw]−[R11R12R13R21R22R23R31R32R33])∥\|\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})\|$
然后将对应的元素相乘, 提出u,v,w 再将矩阵元素拉成竖起条状

$[W11Dxiu+W12Dxiv+W13DxiwW11Dyiu+W12Dyiv+W13DyiwW11Dziu+W12Dziv+W13DziwW21Dxiu+W22Dxiv+W23DxiwW21Dyiu+W22Dyiv+W23DyiwW21Dziu+W22Dziv+W23DziwW31Dxiu+W32Dxiv+W33DxiwW31Dyiu+W32Dyiv+W33DyiwW31Dziu+W32Dziv+W33Dziw]=[W11DxiW12DxiW13DxiW11DyiW12DyiW13DyiW11DziW12DziW13DziW21DxiW22DxiW23DxiW21DyiW22DyiW23DyiW21DziW22DziW23DziW31DxiW32DxiW33DxiW31DyiW32DyiW33DyiW31DziW32DziW33Dzi][uvw]\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}$
将右边的矩阵中i沿#F的维数中竖直展开就变成了矩阵A

同理right hand side 就变成了

IGL_INLINE void buildRhs(igl::SLIMData& s, const Eigen::SparseMatrix<double> &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;
    }
  }
}