...

const double MaxSizeWindow =  120 ;
const double HalfMaxSizeWindow =  0 . 5  * MaxSizeWindow;
const int AmountFixedPointsMin =  5 ;
const double g =  9 . 80665 ;
const double dt =  0 . 01 ;

struct PointPositions
{
	bool fixed;

	double XPrimary, YPrimary, ZPrimary,
		XCurrent, YCurrent, ZCurrent;

	double Mass;

	double ForceX, ForceY, ForceZ;

	double AccelerationX, AccelerationY, AccelerationZ,
		SpeedX, SpeedY, SpeedZ;
};

class ModelingMaterial
{
	int LengthMaterial,
		WidthMaterial,
		PointsX,
		PointsY;
	
	double Ks, Kb, Kt;

	double C;

	double Density;

	PointPositions *tmp;

public:
	ModelingMaterial(const int, const int, const int, const int, 
					const double, const double, const double, const double, const double);
	~ModelingMaterial(void);
	void StartDrawing(void);
	void SetFixedPoints(void);
	void CalculateInsideForce(const int, const int, const int, const int, 
							double &, double &, double &);
	void CalculateAllInsideForces(void);
	void CalculatePointMass(void);
	void CalculateNewPositions(void);
	void StartSimulation(void);
	void CheckDistances(void);
	void ReDrawPicture(void);
};

...

void ModelingMaterial::CheckDistances()
{
	for(int x = PointsY -  1 ; x >=  0 ; x--)
		for(int y = PointsX -  1 ; y >=  0 ; y--)
			for(int k = x -  1 ; k <= x +  1 ; k++)
			{
				if(k <  0  || k > PointsY -  1 ) continue;

				for(int m = y -  1 ; m <= y +  1 ; m++)
				{
					if(m <  0  || m > PointsX -  1 ) continue;
					
					if(k == x && m == y) continue;

					if((k == x +  1  || k == x -  1 ) && (m == y +  1  || m == y -  1 )) continue;

					while(tmp[x * PointsX + y].ZCurrent < tmp[k * PointsX + m].ZCurrent -  0 . 02 )
					{	
							tmp[k * PointsX + m].ZCurrent -=  0 . 001 ;
					}
				}
			}
}

void ModelingMaterial::CalculateNewPositions()
{
	double Fx, Fy, Fz, Ax, Ay, Az, Vx, Vy, Vz;

	for(int i =  0 ; i < PointsY; i++)
		for(int j =  0 ; j < PointsX; j++)
		{
			if(tmp[i * PointsX + j].fixed == true) continue;

			Fx = tmp[i * PointsX + j].ForceX + C * tmp[i * PointsX + j].SpeedX;

			Fy = tmp[i * PointsX + j].ForceY + C * tmp[i * PointsX + j].SpeedY;

			Fz = tmp[i * PointsX + j].ForceZ + C * tmp[i * PointsX + j].SpeedZ + 
				tmp[i * PointsX + j].Mass * g;

			Ax = Fx / tmp[i * PointsX + j].Mass;

			Ay = Fy / tmp[i * PointsX + j].Mass;

			Az = Fz / tmp[i * PointsX + j].Mass;

			Vx = tmp[i * PointsX + j].SpeedX + dt /  2  
				* (tmp[i * PointsX + j].AccelerationX + Ax);

			Vy = tmp[i * PointsX + j].SpeedY + dt /  2  
				* (tmp[i * PointsX + j].AccelerationY + Ay);

			Vz = tmp[i * PointsX + j].SpeedZ + dt /  2  
				* (tmp[i * PointsX + j].AccelerationZ + Az);

			tmp[i * PointsX + j].XCurrent += dt /  2  * (tmp[i * PointsX + j].SpeedX + Vx);
			tmp[i * PointsX + j].YCurrent += dt /  2  * (tmp[i * PointsX + j].SpeedY + Vy);
			tmp[i * PointsX + j].ZCurrent += dt /  2  * (tmp[i * PointsX + j].SpeedZ + Vz);


			tmp[i * PointsX + j].AccelerationX = Ax;
			tmp[i * PointsX + j].AccelerationY = Ay;
			tmp[i * PointsX + j].AccelerationZ = Az;

			tmp[i * PointsX + j].SpeedX = Vx;
			tmp[i * PointsX + j].SpeedY = Vy;
			tmp[i * PointsX + j].SpeedZ = Vz;
		}

		CalculateAllInsideForces();
}

void ModelingMaterial::StartSimulation()
{
	for(double t = dt; t <  1 . 1 ; t += dt)
	{
		CalculateNewPositions();
		CheckDistances();
		ReDrawPicture();
	}
}

void reshape(int w, int h)
{
	glViewport( 0 ,  0 , w, h);
  
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();

	if (w >= h)
	{
		double f = static_cast<double>(h) / static_cast<double>(w);
		glOrtho(-HalfMaxSizeWindow, HalfMaxSizeWindow, 
				-HalfMaxSizeWindow * f, HalfMaxSizeWindow * f, - 500 ,  500 );
	}
	else
	{
		double f = static_cast<double>(w) / static_cast<double>(h);
		glOrtho(-HalfMaxSizeWindow * f, HalfMaxSizeWindow * f, 
				-HalfMaxSizeWindow, HalfMaxSizeWindow, - 500 ,  500 );
	}

	glMatrixMode(GL_MODELVIEW);
}

void display(void)
{
	ModelingMaterial ob( 80 ,  50 ,  40 ,  25 ,  50 ,  50 ,  50 ,  0 . 150 , - 0 . 5 ); //density of cotton kg/(m*m)
	ob.StartDrawing();
	ob.SetFixedPoints();
	ob.CalculatePointMass();
	ob.StartSimulation();
}

ModelingMaterial::ModelingMaterial(const int LM, const int WM, const int PX, const int PY, 
									const double ks, const double kb, const double kt, 
									const double D, const double c)
{
	LengthMaterial = LM;
	WidthMaterial = WM;
	PointsX = PX;
	PointsY = PY;
	Ks = ks;
	Kb = kb;
	Kt = kt;
	Density = D;
	C = c;


	tmp = new PointPositions[PointsX * PointsY];

	if(!tmp) exit(- 1 );

	for(int i =  0 ; i < PointsY; i++)
		for(int j =  0 ; j < PointsX; j++)
	{
			tmp[i * PointsX + j].XPrimary =  0 ;
			tmp[i * PointsX + j].YPrimary =  0 ;
			tmp[i * PointsX + j].ZPrimary =  0 ;

			tmp[i * PointsX + j].XCurrent =  0 ;
			tmp[i * PointsX + j].YCurrent =  0 ;
			tmp[i * PointsX + j].ZCurrent =  0 ;

			tmp[i * PointsX + j].fixed = false;

			tmp[i * PointsX + j].Mass =  0 ;

			tmp[i * PointsX + j].ForceX =  0 ; 
			tmp[i * PointsX + j].ForceY =  0 ; 
			tmp[i * PointsX + j].ForceZ =  0 ; 

			tmp[i * PointsX + j].AccelerationX =  0 ; 
			tmp[i * PointsX + j].AccelerationY =  0 ; 
			tmp[i * PointsX + j].AccelerationZ =  0 ;

			tmp[i * PointsX + j].SpeedX =  0 ; 
			tmp[i * PointsX + j].SpeedY =  0 ; 
			tmp[i * PointsX + j].SpeedZ =  0 ;
	}
}