 	fext_h.push_back(-rho * A * a * a * GRAVITYx / 0.12e2);
 	fext_h.push_back(-rho * A * a * a * GRAVITYy / 0.12e2);
 	fext_h.push_back(-rho * A * a * a * GRAVITYz / 0.12e2);
-	
-//	//add mass matrix
-//	vector<float3> mass = addMassMatrix(rho,A,a);
-//
-//	for(int i=0;i<mass.size();i++)
-//	{
-//		massI_h.push_back(mass[i].x+12*(elements.size()-1));
-//		massJ_h.push_back(mass[i].y+12*(elements.size()-1));
-//		mass_h.push_back(static_cast<double>(mass[i].z));
-//	}
 
 	for(int i=0;i<12;i++)
 	{

	return 0;
}

//int ANCFSystem::createMass() // used in Newton iteration, nice to keep it separate (but not memory efficient)
//{
//	massI_d = massI_h;
//	massJ_d = massJ_h;
//	mass_d = mass_h;
//
//	// create mass matrix using cusp library (shouldn't change)
//	thrust::device_ptr<int> wrapped_device_I(CASTI1(massI_d));
//	DeviceIndexArrayView row_indices = DeviceIndexArrayView(wrapped_device_I, wrapped_device_I + massI_d.size());
//
//	thrust::device_ptr<int> wrapped_device_J(CASTI1(massJ_d));
//	DeviceIndexArrayView column_indices = DeviceIndexArrayView(wrapped_device_J, wrapped_device_J + massJ_d.size());
//
//	thrust::device_ptr<double> wrapped_device_V(CASTD1(mass_d));
//	DeviceValueArrayView values = DeviceValueArrayView(wrapped_device_V, wrapped_device_V + mass_d.size());
//
//	mass = DeviceView(p_h.size(), p_h.size(), mass_d.size(), row_indices, column_indices, values);
//	mass.sort_by_row();
//
//	return 0;
//}

int ANCFSystem::initializeSystem()
{
	ANCFSystem::updatePhiq();
	ANCFSystem::calculateInitialPhi();
	//ANCFSystem::createMass();

	for(int i=0;i<constraints.size();i++)
	{

{
	this->partitions = partitions;
	return 0;
}

int ANCFSystem::addParticle(Particle* particle)
