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

MatrixXi duplicatematrix(int n){
	MatrixXi D;
	D.setZero(n*n,n*(n-1)/2 + n);
	int k = n;
	int l = n;
	for(int i=0;i<n;i++){
		D(i*n+i,i) = 1;
		for(int j=i+1;j<n;j++){
			D(i*n + j,k++) = 1;
      D(j*n + i,l++) = 1;
		}
	
	}
	return D;
}

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());
}

int Simulation::initializeSimulation(double deltaT, int iterations, char method, MatrixXi& TT, MatrixXd& TV, MatrixXd& B, vector<int>& moveVertices, vector<int> fixVertices, double youngs, double poissons){
	iters = iterations;

	if (method =='e'){
		integrator = new Verlet();
		cout<<"Initialized Verlet"<<endl;	
	}else if(method == 'i'){
		integrator = new ImplicitEuler();
		cout<<"Initialized Implicit Euler"<<endl;
	}
	else if(method == 'n'){
		integrator = new ImplicitNewmark();
		cout<<"Initialized Implicit Newmark"<<endl;
	}
	else{
		cout<<"Method not supported yet"<<endl;
		exit(0);
	}
	VectorXd force;
	force.resize(3*TV.rows());
	force.setZero();

	setInitPosition(force, fixVertices, moveVertices);

	if(moveVertices.size()>0 or fixVertices.size()>0){
		//cout << "DOING STUFFS" << endl;
		MatrixXd newTV;
		newTV.resize(TV.rows(), TV.cols());
		newTV.setZero();
		MatrixXi newTT;
		newTT.resize(TT.rows(), TT.cols());
		newTT.setZero();
		//cout << "MoveVertsSize :: " << moveVertices.size() << endl;

		//TODO: Make this shit more efficient
		//Hash maps or something
		vector<int> vertexNewIndices;
		for(int i=0; i<TV.rows(); i++){
			bool flag =false;
			for(unsigned int j=0; j<fixVertices.size(); j++){
				if(i==fixVertices[j]){
					flag = true;
				}
			}
			for(unsigned int j=0; j<moveVertices.size(); j++){
				if(i==moveVertices[j]){
					flag = true;
				}
			}
			// if vertex not fixed or moved, re-index to front
			//[v, v, v, v...., f, f, f...., m, m, m...,m]
			if(!flag){
				vertexNewIndices.push_back(i);
			}
		}
		//re-index fixed verts
		for(unsigned int j=0; j<fixVertices.size(); j++){
			vertexNewIndices.push_back(fixVertices[j]);
		}
		//re-index move verts
		for(unsigned int j=0; j<moveVertices.size(); j++){
			vertexNewIndices.push_back(moveVertices[j]);
		}

		//these are the new indices for the fixed verts
		vector<int> newfixIndices;
		for(unsigned int i= vertexNewIndices.size() - (moveVertices.size() + fixVertices.size()); i<(vertexNewIndices.size()-moveVertices.size()); i++){
			newfixIndices.push_back(i);
		}

		//new indices for the moving verts
		vector<int> newMoveIndices;
		for(unsigned int i= vertexNewIndices.size() - moveVertices.size(); i<vertexNewIndices.size(); i++){
			newMoveIndices.push_back(i);
		}

		//cout << "NewMoveIndicesSize :: " << newMoveIndices.size() << endl;

		VectorXd new_force;
		new_force.resize(3*TV.rows());
		reIndexTVandTT(vertexNewIndices, fixVertices.size(), moveVertices.size(), TV, TT, force, newTV, newTT, new_force);

		igl::barycenter(newTV, newTT, B);
		//Initialize Solid Mesh
		M.initializeMesh(newTT, newTV, youngs, poissons);
		if(moveVertices.size() != 0){
			// cout<<"Move vertices "<<moveVertices.size()<<endl;
			// cout<<"fix verts "<<fixVertices.size()<<endl;
			binarySearchYoungs(newMoveIndices, newTV, newTT, fixVertices.size(), B);
			// syntheticTests(newMoveIndices, newTV, newTT, fixVertices.size(), B);
		
		}
		
		integrator->initializeIntegrator(deltaT, M, newTV, newTT);
		this->external_force = new_force;
		integrator->fixVertices(newfixIndices);
		int ignorePastIndex = newTV.rows() - newfixIndices.size();
		staticSolveNewtonsForces(newTV, newTT, B, new_force, ignorePastIndex);


//.........这里部分代码省略.........