 
 		oglcamera.Update();
 
-		for (int i = 0; i < sys.particles.size(); i++) {
+		for (int i = 0; i < sys->particles.size(); i++) {
 			glColor3f(0.0f, 1.0f, 0.0f);
 			glPushMatrix();
-			float3 pos = sys.getXYZPositionParticle(i);
+			float3 pos = sys->getXYZPositionParticle(i);
 			glTranslatef(pos.x, pos.y, pos.z);
-			glutSolidSphere(sys.particles[i].getRadius(), 30, 30);
+			glutSolidSphere(sys->particles[i].getRadius(), 30, 30);
 			glPopMatrix();
 
 			//indicate velocity
-			glLineWidth(sys.elements[i].getRadius()*500);
+			glLineWidth(sys->elements[i].getRadius()*500);
 			glColor3f(1.0f,0.0f,0.0f);
 			glBegin(GL_LINES);
 			glVertex3f(pos.x,pos.y,pos.z);
-			float3 vel = sys.getXYZVelocityParticle(i);
+			float3 vel = sys->getXYZVelocityParticle(i);
 			//cout << "v:" << vel.x << " " << vel.y << " " << vel.z << endl;
-			pos +=2*sys.particles[i].getRadius()*normalize(vel);
+			pos +=2*sys->particles[i].getRadius()*normalize(vel);
 			glVertex3f(pos.x,pos.y,pos.z);
 			glEnd();
 			glFlush();
 		}
 
-		for(int i=0;i<sys.elements.size();i++)
+		for(int i=0;i<sys->elements.size();i++)
 		{
-			int xiDiv = sys.numContactPoints;
+			int xiDiv = sys->numContactPoints;
 
 			double xiInc = 1/(static_cast<double>(xiDiv-1));
 

				for(int j=0;j<xiDiv;j++)
				{
					glPushMatrix();
					float3 position = sys.getXYZPosition(i,xiInc*j);
					glTranslatef(position.x,position.y,position.z);
					glutSolidSphere(sys.elements[i].getRadius(),10,10);
					glPopMatrix();
				}
			}
			else
			{
				int xiDiv = sys.numContactPoints;
				double xiInc = 1/(static_cast<double>(xiDiv-1));
				glLineWidth(sys.elements[i].getRadius()*500);
				glColor3f(0.0f,1.0f,0.0f);
				glBegin(GL_LINE_STRIP);
				for(int j=0;j<sys.numContactPoints;j++)
				{
					float3 position = sys.getXYZPosition(i,xiInc*j);
					glVertex3f(position.x,position.y,position.z);
				}
				glEnd();

	
	// if you don't want to visualize, then output the data
	int fileIndex = 0;
	while(sys->time < t_end)
	{
		if(sys->getTimeIndex()%outputInterval==0)
		{
			stringstream ss;
			//cout << "Frame: " << fileIndex << endl;
			ss << data_folder << "/" << fileIndex << ".txt";
			sys->writeToFile(ss.str());
			fileIndex++;
		}

    // Force a preconditioner update if needed
    if ((sys->preconditionerUpdateModulus > 0) && (sys->timeIndex % sys->preconditionerUpdateModulus == 0)) {
      //mySolver->update(lhs.values);
      delete sys->mySolver;
      sys->mySolver = new SpikeSolver(sys->partitions, sys->solverOptions);
      sys->mySolver->setup(sys->lhs);
      sys->precUpdated = true;
      printf("Preconditioner updated (step condition)!\n");
    }

		sys->DoTimeStep();
		ofile << sys->time                 << ", "
		      << sys->stepTime             << ", "
		      << sys->stepNewtonIterations << ", "
		      << sys->stepKrylovIterations << ", "
		      << sys->precUpdated          << " ,     ";
		for (size_t i = 0; i < sys->stepNewtonIterations; ++i)
			ofile << sys->spikeSolveTime[i] << ", " << sys->spikeNumIter[i] << ",     ";
		ofile << endl;
	}
	printf("Total time to simulate: %f [s]\n",sys->timeToSimulate);
	ofile.close();

	return 0;
