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

vector<float3> addMassMatrix(double rho, double A, double l)
{
	vector < float3 > massNew;
	float3 massEntry;
	massEntry.x = 0;
	massEntry.y = 0;
	massEntry.z = 0.13e2 / 0.35e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 0;
	massEntry.y = 3;
	massEntry.z = 0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 0;
	massEntry.y = 6;
	massEntry.z = 0.9e1 / 0.70e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 0;
	massEntry.y = 9;
	massEntry.z = -0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 1;
	massEntry.y = 1;
	massEntry.z = 0.13e2 / 0.35e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 1;
	massEntry.y = 4;
	massEntry.z = 0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 1;
	massEntry.y = 7;
	massEntry.z = 0.9e1 / 0.70e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 1;
	massEntry.y = 10;
	massEntry.z = -0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 2;
	massEntry.y = 2;
	massEntry.z = 0.13e2 / 0.35e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 2;
	massEntry.y = 5;
	massEntry.z = 0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 2;
	massEntry.y = 8;
	massEntry.z = 0.9e1 / 0.70e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 2;
	massEntry.y = 11;
	massEntry.z = -0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 3;
	massEntry.y = 0;
	massEntry.z = 0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 3;
	massEntry.y = 3;
	massEntry.z = rho * A * l * l * l / 0.105e3;
	massNew.push_back(massEntry);
	massEntry.x = 3;
	massEntry.y = 6;
	massEntry.z = 0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 3;
	massEntry.y = 9;
	massEntry.z = -rho * A * l * l * l / 0.140e3;
	massNew.push_back(massEntry);
	massEntry.x = 4;
	massEntry.y = 1;
	massEntry.z = 0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 4;
	massEntry.y = 4;
	massEntry.z = rho * A * l * l * l / 0.105e3;
	massNew.push_back(massEntry);
	massEntry.x = 4;
	massEntry.y = 7;
	massEntry.z = 0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 4;
	massEntry.y = 10;
	massEntry.z = -rho * A * l * l * l / 0.140e3;
	massNew.push_back(massEntry);
	massEntry.x = 5;
	massEntry.y = 2;
	massEntry.z = 0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 5;
	massEntry.y = 5;
	massEntry.z = rho * A * l * l * l / 0.105e3;
	massNew.push_back(massEntry);
	massEntry.x = 5;
	massEntry.y = 8;
	massEntry.z = 0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 5;
	massEntry.y = 11;
	massEntry.z = -rho * A * l * l * l / 0.140e3;
	massNew.push_back(massEntry);
	massEntry.x = 6;
	massEntry.y = 0;
	massEntry.z = 0.9e1 / 0.70e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 6;
	massEntry.y = 3;
	massEntry.z = 0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 6;
	massEntry.y = 6;
	massEntry.z = 0.13e2 / 0.35e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 6;
	massEntry.y = 9;
	massEntry.z = -0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 7;
	massEntry.y = 1;
	massEntry.z = 0.9e1 / 0.70e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 7;
	massEntry.y = 4;
	massEntry.z = 0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 7;
	massEntry.y = 7;
	massEntry.z = 0.13e2 / 0.35e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 7;
	massEntry.y = 10;
	massEntry.z = -0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 8;
	massEntry.y = 2;
	massEntry.z = 0.9e1 / 0.70e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 8;
	massEntry.y = 5;
	massEntry.z = 0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 8;
	massEntry.y = 8;
	massEntry.z = 0.13e2 / 0.35e2 * rho * A * l;
	massNew.push_back(massEntry);
	massEntry.x = 8;
	massEntry.y = 11;
	massEntry.z = -0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 9;
	massEntry.y = 0;
	massEntry.z = -0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 9;
	massEntry.y = 3;
	massEntry.z = -rho * A * l * l * l / 0.140e3;
	massNew.push_back(massEntry);
	massEntry.x = 9;
	massEntry.y = 6;
	massEntry.z = -0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 9;
	massEntry.y = 9;
	massEntry.z = rho * A * l * l * l / 0.105e3;
	massNew.push_back(massEntry);
	massEntry.x = 10;
	massEntry.y = 1;
	massEntry.z = -0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 10;
	massEntry.y = 4;
	massEntry.z = -rho * A * l * l * l / 0.140e3;
	massNew.push_back(massEntry);
	massEntry.x = 10;
	massEntry.y = 7;
	massEntry.z = -0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 10;
	massEntry.y = 10;
	massEntry.z = rho * A * l * l * l / 0.105e3;
	massNew.push_back(massEntry);
	massEntry.x = 11;
	massEntry.y = 2;
	massEntry.z = -0.13e2 / 0.420e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 11;
	massEntry.y = 5;
	massEntry.z = -rho * A * l * l * l / 0.140e3;
	massNew.push_back(massEntry);
	massEntry.x = 11;
	massEntry.y = 8;
	massEntry.z = -0.11e2 / 0.210e3 * rho * A * l * l;
	massNew.push_back(massEntry);
	massEntry.x = 11;
	massEntry.y = 11;
	massEntry.z = rho * A * l * l * l / 0.105e3;
	massNew.push_back(massEntry);

	return massNew;
}

int ANCFSystem::addParticle(Particle* particle)

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

	for(int i=0;i<12;i++)
	{
