C++ VectorXd::setZero方法代码示例

    /* Solve min||Ax-b|| for a matrix A whose rank is given. */
    VectorXd solve( const VectorXd& b , bool inY=false)
    {
      if( rank==0 ) return VectorXd::Zero(NC);

      /* Approximate solution using no basis transfo (result is meanigless
       * appart from the computation time pov. */
      /*
	VectorXd sol = b.head(rank);
	matrixL().solveInPlace( sol );
	VectorXd res; res.setZero(NC);
	res.head(rank)=sol; return res;
      */

      /* With plain matrices. */
      /*
	VectorXd sol = matrixUr().transpose()*b;
	matrixL().solveInPlace( sol );
	return matrixVr()*sol;
      */

      /* Using the HH representation of V. */
      assert( m_computeThinU || m_computeFullU );
      VectorXd sol; 
      if( inY ) sol.setZero(rank); else sol.setZero(NC); 
      sol.head(rank) = matrixUr().transpose()*b; 
      matrixL().solveInPlace( sol.head(rank) );
      if( ! inY ) sol.applyOnTheLeft(qrv.householderQ());
      return sol;
   }

 ConjugateGradientType(Functional f, const GridDescription &gdin, double kT, VectorXd *data, LineMinimizer lm,
                       double stepsize = 0.1)
   : MinimizerInterface(f, gdin, kT, data), step(stepsize), orig_step(step), direction(*data), oldgrad(*data), linmin(lm) {
   direction.setZero();
   oldgrad.setZero();
   oldgradsqr = 0;
 }

 void minimize(Functional newf, const GridDescription &gdnew, VectorXd *newx = 0) {
   step = orig_step;
   MinimizerInterface::minimize(newf, gdnew, newx);
   if (newx) {
     direction = *newx;
     oldgrad = *newx;
   }
   direction.setZero();
   oldgrad.setZero();
   oldgradsqr = 0;
 }

/** 
 * reconstruct the displacements u in Euler space with RS coordinates y
 * provided.
 * 
 * @param y the RS coordinates.
 * @param u the constructed Euler coordinates represents the displacements.
 * 
 * @return true if construction is success.
 */
bool RS2Euler::reconstruct(const VectorXd &y, VectorXd &u){

  assert (tetmesh != NULL);
  const int node_numx3 = tetmesh->nodes().size()*3;
  bool succ = true;

  // assemble the right_side
  VectorXd b;
  assemble_b(y,b);
  assert_eq(VG_t.cols(),b.size());
  assert_eq(VG_t.rows(),node_numx3);

  VectorXd right_side(node_numx3 + numFixedDofs());
  right_side.segment(0, node_numx3) = VG_t*b;
  right_side.segment(node_numx3,numFixedDofs()).setZero();
  right_side.segment(node_numx3,barycenter_uc.size()) = barycenter_uc;

  // solve A*u = right_side
  assert_eq (right_side.size(), A_solver.rows());
  u = A_solver.solve(right_side);

  // get the result 
  if(A_solver.info()!=Eigen::Success) {
	
  	succ = false;
  	u.resize(node_numx3);
  	u.setZero();
  	ERROR_LOG("failed to solve for A X = P.");
  }else{
	assert_gt(u.size(), node_numx3);
	const VectorXd x = u.head(node_numx3);
	u = x;
  }
  return succ;
}

void ImplicitEuler::ImplicitCalculateForces( MatrixXd& TVk, SparseMatrix<double>& forceGradient, VectorXd& x_k, VectorXd& f){
	// //gravity
	f.setZero();
	for(unsigned int i=0; i<M.tets.size(); i++){
		double vertex_mass = M.tets[i].undeformedVol/4;//assume const density 1
		Vector4i indices = M.tets[i].verticesIndex;
		f(3*indices(0)+1) += vertex_mass*gravity;
		f(3*indices(1)+1) += vertex_mass*gravity; 
		f(3*indices(2)+1) += vertex_mass*gravity;
		f(3*indices(3)+1) += vertex_mass*gravity;
	}

	//elastic
	for(unsigned int i=0; i<M.tets.size(); i++){
		Vector4i indices = M.tets[i].verticesIndex;
		MatrixXd F_tet = M.tets[i].computeElasticForces(TVk, simTime%2);
		f.segment<3>(3*indices(0)) += F_tet.col(0);
		f.segment<3>(3*indices(1)) += F_tet.col(1);
		f.segment<3>(3*indices(2)) += F_tet.col(2);
		f.segment<3>(3*indices(3)) += F_tet.col(3);
	}
	// cout<<f<<endl<<endl;
	//damping
	f += rayleighCoeff*forceGradient*(x_k - x_old)/h;
	// cout<<f<<endl<<endl;
	return;
}

void BuildCausalVector(VectorXd& vec2Build , VectorXd& index){
    vec2Build.setZero();
    for(int i = 0; i <index.size(); i++){
        if(index[i] >0){
            vec2Build[index[i]-1] = 1;
        }
    }
}

VectorXd vec(MatrixXd& M){
	VectorXd V;
	V.setZero(M.rows()*M.cols());
	for(int i=0;i<M.cols();i++){
		V.segment(i*M.rows(), M.rows()) = M.col(i);
	}
	return V;
}

VectorXd Tetra10::Give_weights()
{
    VectorXd w;
    w.setZero(ngauss);

    w(0) = 0.0416666666666667;
    w(1) = 0.0416666666666667;
    w(2) = 0.0416666666666667;
    w(3) = 0.0416666666666667;
    return w;
}

VectorXd vech(MatrixXd& M){
	int n = M.rows();
	VectorXd V;
	V.setZero(n*(n-1)/2+n);
	V.segment(0,n) = M.diagonal();
	int k = n;
	for(int i=0;i<n-1;i++){
		V.segment(k,n-i-1) = M.block(i,i+1,1,n-i-1).transpose();
		k += n-i-1;
	}
	return V;
}

void Simulation::calculateElasticForces(VectorXd &f, MatrixXd& TV){
	f.setZero();
	//elastic
	for(unsigned int i=0; i< M.tets.size(); i++){
		Vector4i indices = M.tets[i].verticesIndex;
		MatrixXd F_tet = M.tets[i].computeElasticForces(TV, 1);
		f.segment<3>(3*indices(0)) += F_tet.col(0);
		f.segment<3>(3*indices(1)) += F_tet.col(1);
		f.segment<3>(3*indices(2)) += F_tet.col(2);
		f.segment<3>(3*indices(3)) += F_tet.col(3);
	}
	return;
}

VectorXd Triangle<ConcreteShape>::QuadratureIntegrationWeight(const int order) {

  VectorXd wn;
  if(order == 3) {
    int num_pts = 12;
    wn.setZero(num_pts);
    quadrature_weights_p3_triangle(wn.data());
  } else {
    std::cerr << "ERROR: Order NOT implemented!\n";
    MPI::COMM_WORLD.Abort(-1);
  }
  return wn;
}

const VectorXd &AniEditDM::getUforConstraint(const int frame)const{
  
  static VectorXd tempt_u0;
  if (interpolator != NULL && totalFrameNum() > 0){
	assert_in(frame,0, totalFrameNum()-1);
	return interpolator->getUforConstraint(frame);
  }else{
	if (tempt_u0.size() != fullDim()){
	  tempt_u0.resize(fullDim());
	  tempt_u0.setZero();
	}
	return tempt_u0;
  }
}

VectorXd multivariateNormal_density(const MatrixXd &x,
                                    const MatrixXd &Mu, const MatrixXd &Sigma) {
    if ( x.rows() != Mu.rows()) {
        std::cerr << "FROM multivariateNormal_density." << std::endl;
        std::cerr << "ERROR: X_ROWS != MU_ROWS." << std::endl;
    }
    VectorXd res;
    int T = x.rows();
    res.setZero(T);
    for(int i = 0; i<x.rows(); i++) {
        res(i) = double_multivariateNormal_density(x.row(i),
                 Mu.row(i),
                 Sigma);
    }
    return res;
}

void ImplicitEuler::ImplicitTVtoX(VectorXd& x_tv, MatrixXd& TVk){
	x_tv.setZero();
	for(unsigned int i = 0; i < M.tets.size(); i++){
		Vector4i indices = M.tets[i].verticesIndex;

		x_tv(3*indices(0)) = TVk.row(indices(0))[0];
		x_tv(3*indices(0)+1) = TVk.row(indices(0))[1];
		x_tv(3*indices(0)+2) = TVk.row(indices(0))[2];

		x_tv(3*indices(1)) = TVk.row(indices(1))[0];
		x_tv(3*indices(1)+1) = TVk.row(indices(1))[1];
		x_tv(3*indices(1)+2) = TVk.row(indices(1))[2];

		x_tv(3*indices(2)) = TVk.row(indices(2))[0];
		x_tv(3*indices(2)+1) = TVk.row(indices(2))[1];
		x_tv(3*indices(2)+2) = TVk.row(indices(2))[2];

		x_tv(3*indices(3)) = TVk.row(indices(3))[0];
		x_tv(3*indices(3)+1) = TVk.row(indices(3))[1];
		x_tv(3*indices(3)+2) = TVk.row(indices(3))[2];
	}
}

void read_SparseMatrix(SparseMatrix<double,0,int>& M, string path,string name){

	FILE * pFile;
	pFile = fopen((path + name+"i.bin").c_str(), "rb");
	if (pFile==NULL) {fputs ("File error",stderr); exit (1);}

	int n,m,ni;
	if(1 != fread(&n,sizeof(int),1,pFile)){
		fputs ("Read error\n",stderr); exit (1);
	}
	if(1 != fread(&m,sizeof(int),1,pFile)){
		fputs ("Read error\n",stderr); exit (1);
	}

	if(1 != fread(&ni,sizeof(int),1,pFile)){
		fputs ("Read error\n",stderr); exit (1);
	}
	M.resize(n,m);
	MatrixXi ij;
	ij.setZero(ni,2);
	if(ni*2 != fread(ij.data(),sizeof(int),ni*2,pFile)){
		fputs ("Read error\n",stderr); exit (1);
	}
	fclose(pFile);
	VectorXd v;
	v.setZero(ni);
	pFile = fopen((path + name+"v.bin").c_str(), "rb");
	if (pFile==NULL) {fputs ("File error",stderr); exit (1);}
	if(ni != fread(v.data(),sizeof(double),ni,pFile)){
		fputs ("Read error\n",stderr); exit (1);
	}
	fclose(pFile);
	vector<T> coef;
	for(int i=0;i<ni;i++){
		coef.push_back(T(ij(i,0),ij(i,1),v(i)));
	}
	M.setFromTriplets(coef.begin(), coef.end());
}