class CDestructionBody2: public CGradationBody2{
public:

	// for construction
	CDestructionBody2(){ }
	CDestructionBody2( CVoxel2 &voxel, Double centerX, Double centerY, Double angle, Double scale = 1.0 ){ this->construct( voxel, centerX, centerY, angle, scale ); }
	void construct( CVoxel2 &voxel, Double centerX, Double centerY, Double angle, Double scale = 1.0 );
	void renew( Double gravityConstant, Double youngsModulus, Double dampingCoefficient, Double timeStep, Double threshold = -1.0 );
	void disconnect( int x, int y, Double threshold );
	void disconnect( Double threshold );
	bool ifDisconnected( int x, int y ){ return joint[jsn(x,y)] && joint[jsn(x,y)]->disconnected; }
	bool ifDisconnected( CElement2 *doubleElement, int x, int y, int n );
	void reconstruct( Double youngsModulus, Double dampingCoefficient );
	bool checkReconstruction( CElement2 *doubleElement, int level, int x, int y, int n );
	void reconstructPhase1( CElement2 *doubleElement, int level, int x, int y, int n, Double youngsModulus, Double dampingCoefficient );
	bool reconstructPhase2( CElement2 *doubleElement, int level, int x, int y, int n, Double youngsModulus, Double dampingCoefficient );
};

class CDestructionBody: public CGradationBody{
public:

	// for construction
	CDestructionBody(){ }
	CDestructionBody( CVoxel &voxel, Double centerX, Double centerY, Double centerZ, Double angle, Double scale = 1.0 ){ this->construct( voxel, centerX, centerY, centerZ, angle, scale ); }
	void construct( CVoxel &voxel, Double centerX, Double centerY, Double centerZ, Double angle, Double scale = 1.0 );
	void renew( Double gravityConstant, Double youngsModulus, Double dampingCoefficient, Double timeStep, Double threshold = -1.0 );
	void disconnect( int x, int y, int z, Double threshold );
	void disconnect( Double threshold );
	bool ifDisconnected( int x, int y, int z ){ return joint[jsn(x,y,z)] && joint[jsn(x,y,z)]->disconnected; }
	bool ifDisconnected( CElement *doubleElement, int x, int y, int z, int n );
	void reconstruct( Double youngsModulus, Double dampingCoefficient );
	bool checkReconstruction( CElement *doubleElement, int level, int x, int y, int z, int n );
	void reconstructPhase1( CElement *doubleElement, int level, int x, int y, int z, int n, Double youngsModulus, Double dampingCoefficient );
	bool reconstructPhase2( CElement *doubleElement, int level, int x, int y, int z, int n, Double youngsModulus, Double dampingCoefficient );
};

void CDestructionBody2::construct( CVoxel2 &voxel, Double centerX, Double centerY, Double angle, Double scale ){
	this->constructSizes( voxel.xSize, voxel.ySize );
//	levelSize = 1;//for test
	this->constructElement( elementMarks[levelSize] );


	int elementCount = 0;
	int y;
	for( y=0; y<ySizes[0]; y++ ){
	for( int x=0; x<xSizes[0]; x++ ){
		if( voxel.ifSurface( x, y ) ){
			element[sn(0, x, y)] = new CElement2( x, y, centerX, centerY, angle, scale );
			elementCount++;
		}
	}}

	CVoxel2 voxel2( voxel );
	for( y=0; y<ySizes[0]; y++ ){
	for( int x=0; x<xSizes[0]; x++ ){
		if( element[sn(0,x,y)] ) voxel2.set( x, y, 0 );
	}}
//	voxel.print();

	printf( "elementCount = %d\n", elementCount + this->registerElement( voxel2, centerX, centerY, angle, scale ) );

	this->constructJoint();
	CVoxelBody2::registerJoint();
	CGradationBody2::registerJoint();
}
void CDestructionBody::construct( CVoxel &voxel, Double centerX, Double centerY, Double centerZ, Double angle, Double scale ){
	this->constructSizes( voxel.xSize, voxel.ySize, voxel.zSize );
//	levelSize = 1;//for test
	this->constructElement( elementMarks[levelSize] );


	int elementCount = 0;
	int z;
	for( z=0; z<zSizes[0]; z++ ){
	for( int y=0; y<ySizes[0]; y++ ){
	for( int x=0; x<xSizes[0]; x++ ){
		if( voxel.ifSurface( x, y, z ) ){
			element[sn(0, x, y, z)] = new CElement( x, y, z, centerX, centerY, centerZ, angle, scale );
			elementCount++;
		}
	}}}

	voxel.print();
	for( z=0; z<zSizes[0]; z++ ){
	for( int y=0; y<ySizes[0]; y++ ){
	for( int x=0; x<xSizes[0]; x++ ){
		if( element[sn(0,x,y,z)] ) voxel.set( x, y, z, 0 );
	}}}
	voxel.print();

	printf( "elementCount = %d\n", elementCount + this->registerElement( voxel, centerX, centerY, centerZ, angle, scale ) );

	this->constructJoint();
	CVoxelBody::registerJoint();
	CGradationBody::registerJoint();
}

void CDestructionBody2::renew( Double gravityConstant, Double youngsModulus, Double dampingCoefficient, Double timeStep, Double threshold ){//1.906
	this->setElasticForce( youngsModulus, dampingCoefficient );//0.906
	if( threshold > 0.0 ){
		this->disconnect( threshold );//0.225
		this->reconstruct( youngsModulus, dampingCoefficient );//0.234
	}
	this->adjustAcceleration();//0.18
	this->addGravityAndRenew(gravityConstant, timeStep);//0.344
}
void CDestructionBody::renew( Double gravityConstant, Double youngsModulus, Double dampingCoefficient, Double timeStep, Double threshold ){
	this->setElasticForce( youngsModulus, dampingCoefficient );
	if( threshold > 0.0 ){
		this->disconnect( threshold );
		this->reconstruct( youngsModulus, dampingCoefficient );
	}
	this->adjustAcceleration();
	this->addGravityAndRenew(gravityConstant, timeStep);
}

void CDestructionBody2::disconnect( int x, int y, Double threshold ){
	if( ( x>0 && y>0 && element[sn(0,x-1,y-1)] )
	  || ( x>0 && y<ySize && element[sn(0,x-1,y)] )
	  || ( x<xSize && y>0 && element[sn(0,x,y-1)] )
	  || ( x<xSize && y<ySize && element[sn(0,x,y)] ) ){
		joint[jsn(x,y)]->disconnect( threshold );
	}
}
void CDestructionBody::disconnect( int x, int y, int z, Double threshold ){
	if(  ( x>0 && y>0 && z>0 && element[sn(0,x-1,y-1,z-1)] )
	  || ( x>0 && y<ySize && z>0 && element[sn(0,x-1,y,z-1)] )
	  || ( x>0 && y>0 && z<zSize && element[sn(0,x-1,y-1,z)] )
	  || ( x>0 && y<ySize && z<zSize && element[sn(0,x-1,y,z)] )
	  || ( x<xSize && y>0 && z>0 && element[sn(0,x,y-1,z-1)] )
	  || ( x<xSize && y<ySize && z>0 && element[sn(0,x,y,z-1)] )
	  || ( x<xSize && y>0 && z<zSize && element[sn(0,x,y-1,z)] )
	  || ( x<xSize && y<ySize && z<zSize && element[sn(0,x,y,z)] ) ){
		joint[jsn(x,y,z)]->disconnect( threshold );
	}
}

void CDestructionBody2::disconnect( Double threshold ){//0.225
	for( int y=0; y<ySize+1; y++ ){
	for( int x=0; x<xSize+1; x++ ){
		if( joint[jsn(x,y)] ) this->disconnect( x, y, threshold );
	}}
}
void CDestructionBody::disconnect( Double threshold ){
	for( int z=0; z<zSize+1; z++ ){
	for( int y=0; y<ySize+1; y++ ){
	for( int x=0; x<xSize+1; x++ ){
		if( joint[jsn(x,y,z)] ) this->disconnect( x, y, z, threshold );
	}}}
}

bool CDestructionBody2::ifDisconnected( CElement2 *doubleElement, int x, int y, int n ){
	if( doubleElement== NULL ) return false;
	int xo=x*n;
	int yo=y*n;
	int u = n/2;
	bool disconnected = false;
	for( int j=0; j<3; j++ ){
	for( int i=0; i<3; i++ ){
		if( doubleElement->vertex[sn33(i,j)] && this->ifDisconnected(xo+u*i,yo+u*j) ) disconnected = true;
	}}
	return disconnected;
}
bool CDestructionBody::ifDisconnected( CElement *doubleElement, int x, int y, int z, int n ){
	if( doubleElement== NULL ) return false;
	int xo=x*n;
	int yo=y*n;
	int zo=z*n;
	int u = n/2;
	bool disconnected = false;
	for( int k=0; k<3; k++ ){
	for( int j=0; j<3; j++ ){
	for( int i=0; i<3; i++ ){
		if( doubleElement->vertex[sn333(i,j,k)] && this->ifDisconnected(xo+u*i,yo+u*j,zo+u*k) ) disconnected = true;
	}}}
	return disconnected;
}

void CDestructionBody2::reconstruct( Double youngsModulus, Double dampingCoefficient ){//0.765
	int n = 1;
	for( int level=1; level<levelSize; level++ ){
		n *= 2;

		int y;
		for( y=0; y<ySizes[level]; y++ ){
		for( int x=0; x<xSizes[level]; x++ ){
			if( element[sn(level,x,y)] ){
				if( this->checkReconstruction( element[sn(level,x,y)], level, x, y, n ) ){
					this->reconstructPhase1( element[sn(level,x,y)], level, x, y, n, youngsModulus, dampingCoefficient );
				}
			}
		}}
		for( y=0; y<ySizes[level]; y++ ){
		for( int x=0; x<xSizes[level]; x++ ){
			if( element[sn(level,x,y)] ){
				if( this->reconstructPhase2( element[sn(level,x,y)], level, x, y, n, youngsModulus, dampingCoefficient ) ){
					delete element[sn(level,x,y)];
					element[sn(level,x,y)] = NULL;
				}
			}
		}}
	}
}
void CDestructionBody::reconstruct( Double youngsModulus, Double dampingCoefficient ){
	int n = 1;
	for( int level=1; level<levelSize; level++ ){
		n *= 2;

		int z;
		for( z=0; z<zSizes[level]; z++ ){
		for( int y=0; y<ySizes[level]; y++ ){
		for( int x=0; x<xSizes[level]; x++ ){
			if( element[sn(level,x,y,z)] ){
				if( this->checkReconstruction( element[sn(level,x,y,z)], level, x, y, z, n ) ){
					this->reconstructPhase1( element[sn(level,x,y,z)], level, x, y, z, n, youngsModulus, dampingCoefficient );
				}
			}
		}}}
		for( z=0; z<zSizes[level]; z++ ){
		for( int y=0; y<ySizes[level]; y++ ){
		for( int x=0; x<xSizes[level]; x++ ){
			if( element[sn(level,x,y,z)] ){
				if( this->reconstructPhase2( element[sn(level,x,y,z)], level, x, y, z, n, youngsModulus, dampingCoefficient ) ){
					delete element[sn(level,x,y,z)];
					element[sn(level,x,y,z)] = NULL;
				}
			}
		}}}
	}
}

bool CDestructionBody2::checkReconstruction( CElement2 *doubleElement, int level, int x, int y, int n ){
	int x2=x*2;
	int y2=y*2;
	int x4=x*4;
	int y4=y*4;
	int u = n/2;

	for( int j=0; j<2; j++ ){
	for( int i=0; i<2; i++ ){
		int i2 = i*2;
		int j2 = j*2;
		int i4 = i*4;
		int j4 = j*4;
		if( doubleElement->vertex[sn33(i2,j2)] ){
			if( level == 1 ){
				if( this->ifDisconnected(x2+i2,y2+j2) || this->ifDisconnected(x2+1,y2+j2) || this->ifDisconnected(x2+i2,y2+1) ){
					return true;
				}
			}else{
				bool mx = element[sn(level-2,x4+i2-1,y4+j2)] || element[sn(level-2,x4+i2-1,y4+j2+1)];
				bool mX = element[sn(level-2,x4+i2+2,y4+j2)] || element[sn(level-2,x4+i2+2,y4+j2+1)];
				bool my = element[sn(level-2,x4+i2,y4+j2-1)] || element[sn(level-2,x4+i2+1,y4+j2-1)];
				bool mY = element[sn(level-2,x4+i2,y4+j2+2)] || element[sn(level-2,x4+i2+1,y4+j2+2)];
				if( mx || mX || my || mY ){
					return true;
				}
			}
		}
	}}
	return false;
}
bool CDestructionBody::checkReconstruction( CElement *doubleElement, int level, int x, int y, int z, int n ){
	int x2=x*2;
	int y2=y*2;
	int z2=z*2;
	int x4=x*4;
	int y4=y*4;
	int z4=z*4;
	int u = n/2;

	for( int k=0; k<2; k++ ){
	for( int j=0; j<2; j++ ){
	for( int i=0; i<2; i++ ){
		int i2 = i*2;
		int j2 = j*2;
		int k2 = k*2;
		int i4 = i*4;
		int j4 = j*4;
		int k4 = k*4;
		if( doubleElement->vertex[sn333(i2,j2,k2)] ){
			if( level == 1 ){
				if( this->ifDisconnected(x2+i2,y2+j2,z2+k2) || this->ifDisconnected(x2+1,y2+j2,z2+k2) || this->ifDisconnected(x2+i2,y2+1,z2+k2) || this->ifDisconnected(x2+i2,y2+j2,z2+1)
				|| this->ifDisconnected(x2+i2,y2+1,z2+1) || this->ifDisconnected(x2+1,y2+j2,z2+1) || this->ifDisconnected(x2+1,y2+1,z2+k2) ){
					return true;
				}
			}else{
				bool mx = element[sn(level-2,x4+i2-1,y4+j2,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2+1)];
				bool mX = element[sn(level-2,x4+i2+2,y4+j2,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2+1)];
				bool my = element[sn(level-2,x4+i2,y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2,y4+j2-1,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2+1)];
				bool mY = element[sn(level-2,x4+i2,y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2,y4+j2+2,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2+1)];
				bool mz = element[sn(level-2,x4+i2,y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2-1)];
				bool mZ = element[sn(level-2,x4+i2,y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2+2)];

				bool mxy = element[sn(level-2,x4+i2-1,y4+j2,  z4+k2)] || element[sn(level-2,x4+i2,  y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2,  z4+k2+1)] || element[sn(level-2,x4+i2,  y4+j2-1,z4+k2+1)];
				bool mXy = element[sn(level-2,x4+i2+2,y4+j2,  z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2,  z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2+1)];
				bool mxY = element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2,  y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2+1)] || element[sn(level-2,x4+i2,  y4+j2+2,z4+k2+1)];
				bool mXY = element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2+1)];

				bool myz = element[sn(level-2,x4+i2,y4+j2-1,z4+k2  )] || element[sn(level-2,x4+i2,y4+j2,  z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2  )] || element[sn(level-2,x4+i2+1,y4+j2,  z4+k2-1)];
				bool mYz = element[sn(level-2,x4+i2,y4+j2+2,z4+k2  )] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2  )] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2-1)];
				bool myZ = element[sn(level-2,x4+i2,y4+j2-1,z4+k2+1)] || element[sn(level-2,x4+i2,y4+j2,  z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2,  z4+k2+2)];
				bool mYZ = element[sn(level-2,x4+i2,y4+j2+2,z4+k2+1)] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2+2)];

				bool mzx = element[sn(level-2,x4+i2,  y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2-1,y4+j2,z4+k2  )] || element[sn(level-2,x4+i2  ,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2  )];
				bool mZx = element[sn(level-2,x4+i2,  y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2-1,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2  ,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2+1)];
				bool mzX = element[sn(level-2,x4+i2+1,y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2+2,y4+j2,z4+k2  )] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2  )];
				bool mZX = element[sn(level-2,x4+i2+1,y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2+2,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2+1)];

				if( mx || mX || my || mY || mz || mZ || mxy || mXy || mxY || mXY || myz || mYz || myZ || mYZ || mzx || mZx || mzX || mZX ){
					return true;
				}
			}
		}
	}}}
	return false;
}

void CDestructionBody2::reconstructPhase1( CElement2 *doubleElement, int level, int x, int y, int n, Double youngsModulus, Double dampingCoefficient ){
	doubleElement->constructMiddlePoints();
	int x2=x*2;
	int y2=y*2;
	int x4=x*4;
	int y4=y*4;
	int u = n/2;

	for( int j=0; j<2; j++ ){
	for( int i=0; i<2; i++ ){
		int i2 = i*2;
		int j2 = j*2;
		int i4 = i*4;
		int j4 = j*4;
		if( doubleElement->vertex[sn33(i2,j2)] ){
			if( level == 1 ){
				if( this->ifDisconnected(x2+i2,y2+j2) || this->ifDisconnected(x2+1,y2+j2) || this->ifDisconnected(x2+i2,y2+1) ){
					element[sn(level-1,x2+i,y2+j)] = new CElement2( *doubleElement->vertex[sn33(i,j)], *doubleElement->vertex[sn33(i+1,j)],*doubleElement->vertex[sn33(i,j+1)],*doubleElement->vertex[sn33(i+1,j+1)] );
					CVoxelBody2::registerJoint( element[sn(level-1,x2+i,y2+j)], x2+i, y2+j );
					element[sn(level-1,x2+i,y2+j)]->setElasticForce( youngsModulus, dampingCoefficient );
					doubleElement->deleteVertex( sn33(i2,j2) );
				}
			}else{
				bool mx = element[sn(level-2,x4+i2-1,y4+j2)] || element[sn(level-2,x4+i2-1,y4+j2+1)];
				bool mX = element[sn(level-2,x4+i2+2,y4+j2)] || element[sn(level-2,x4+i2+2,y4+j2+1)];
				bool my = element[sn(level-2,x4+i2,y4+j2-1)] || element[sn(level-2,x4+i2+1,y4+j2-1)];
				bool mY = element[sn(level-2,x4+i2,y4+j2+2)] || element[sn(level-2,x4+i2+1,y4+j2+2)];
				if( mx || mX || my || mY ){
					element[sn(level-1,x2+i,y2+j)] = new CElement2( *doubleElement->vertex[sn33(i,j)], *doubleElement->vertex[sn33(i+1,j)],*doubleElement->vertex[sn33(i,j+1)],*doubleElement->vertex[sn33(i+1,j+1)], mx, mX, my, mY );
					CGradationBody2::registerJoint( element[sn(level-1,x2+i,y2+j)], x2+i, y2+j, u );
					element[sn(level-1,x2+i,y2+j)]->setElasticForce( youngsModulus, dampingCoefficient );
					doubleElement->deleteVertex( sn33(i2,j2) );
				}
			}
		}
	}}
}
void CDestructionBody::reconstructPhase1( CElement *doubleElement, int level, int x, int y, int z, int n, Double youngsModulus, Double dampingCoefficient ){
	doubleElement->constructMiddlePoints();
	int x2=x*2;
	int y2=y*2;
	int z2=z*2;
	int x4=x*4;
	int y4=y*4;
	int z4=z*4;
	int u = n/2;

	for( int k=0; k<2; k++ ){
	for( int j=0; j<2; j++ ){
	for( int i=0; i<2; i++ ){
		int i2 = i*2;
		int j2 = j*2;
		int k2 = k*2;
		int i4 = i*4;
		int j4 = j*4;
		int k4 = k*4;
		if( doubleElement->vertex[sn333(i2,j2,k2)] ){
			if( level == 1 ){
				if( this->ifDisconnected(x2+i2,y2+j2,z2+k2) || this->ifDisconnected(x2+1,y2+j2,z2+k2) || this->ifDisconnected(x2+i2,y2+1,z2+k2) || this->ifDisconnected(x2+i2,y2+j2,z2+1)
				|| this->ifDisconnected(x2+1,y2+1,z2+k2) || this->ifDisconnected(x2+i2,y2+1,z2+1) || this->ifDisconnected(x2+1,y2+j2,z2+1) ){
					element[sn(level-1,x2+i,y2+j,z2+k)] = new CElement( *doubleElement->vertex[sn333(i,j,k)], *doubleElement->vertex[sn333(i+1,j,k)],*doubleElement->vertex[sn333(i,j+1,k)],*doubleElement->vertex[sn333(i+1,j+1,k)],
						                                                *doubleElement->vertex[sn333(i,j,k+1)], *doubleElement->vertex[sn333(i+1,j,k+1)],*doubleElement->vertex[sn333(i,j+1,k+1)],*doubleElement->vertex[sn333(i+1,j+1,k+1)]);
					CVoxelBody::registerJoint( element[sn(level-1,x2+i,y2+j,z2+k)], x2+i, y2+j, z2+k );
					element[sn(level-1,x2+i,y2+j,z2+k)]->setElasticForce( youngsModulus, dampingCoefficient );
					doubleElement->deleteVertex( sn333(i2,j2,k2) );
				}
			}else{
				bool mx = element[sn(level-2,x4+i2-1,y4+j2,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2+1)];
				bool mX = element[sn(level-2,x4+i2+2,y4+j2,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2+1)];
				bool my = element[sn(level-2,x4+i2,y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2,y4+j2-1,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2+1)];
				bool mY = element[sn(level-2,x4+i2,y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2,y4+j2+2,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2+1)];
				bool mz = element[sn(level-2,x4+i2,y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2-1)];
				bool mZ = element[sn(level-2,x4+i2,y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2+2)];

				bool mxy = element[sn(level-2,x4+i2-1,y4+j2,  z4+k2)] || element[sn(level-2,x4+i2,  y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2,  z4+k2+1)] || element[sn(level-2,x4+i2,  y4+j2-1,z4+k2+1)];
				bool mXy = element[sn(level-2,x4+i2+2,y4+j2,  z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2,  z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2+1)];
				bool mxY = element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2,  y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2+1)] || element[sn(level-2,x4+i2,  y4+j2+2,z4+k2+1)];
				bool mXY = element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2+1)];

				bool myz = element[sn(level-2,x4+i2,y4+j2-1,z4+k2  )] || element[sn(level-2,x4+i2,y4+j2,  z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2  )] || element[sn(level-2,x4+i2+1,y4+j2,  z4+k2-1)];
				bool mYz = element[sn(level-2,x4+i2,y4+j2+2,z4+k2  )] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2  )] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2-1)];
				bool myZ = element[sn(level-2,x4+i2,y4+j2-1,z4+k2+1)] || element[sn(level-2,x4+i2,y4+j2,  z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2-1,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2,  z4+k2+2)];
				bool mYZ = element[sn(level-2,x4+i2,y4+j2+2,z4+k2+1)] || element[sn(level-2,x4+i2,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2+1,y4+j2+2,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2+2)];

				bool mzx = element[sn(level-2,x4+i2,  y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2-1,y4+j2,z4+k2  )] || element[sn(level-2,x4+i2  ,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2  )];
				bool mZx = element[sn(level-2,x4+i2,  y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2-1,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2  ,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2-1,y4+j2+1,z4+k2+1)];
				bool mzX = element[sn(level-2,x4+i2+1,y4+j2,z4+k2-1)] || element[sn(level-2,x4+i2+2,y4+j2,z4+k2  )] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2-1)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2  )];
				bool mZX = element[sn(level-2,x4+i2+1,y4+j2,z4+k2+2)] || element[sn(level-2,x4+i2+2,y4+j2,z4+k2+1)] || element[sn(level-2,x4+i2+1,y4+j2+1,z4+k2+2)] || element[sn(level-2,x4+i2+2,y4+j2+1,z4+k2+1)];

				if( mx || mX || my || mY || mz || mZ || mxy || mXy || mxY || mXY || myz || mYz || myZ || mYZ || mzx || mZx || mzX || mZX ){
					element[sn(level-1,x2+i,y2+j,z2+k)] = new CElement( *doubleElement->vertex[sn333(i,j,k)], *doubleElement->vertex[sn333(i+1,j,k)],*doubleElement->vertex[sn333(i,j+1,k)],*doubleElement->vertex[sn333(i+1,j+1,k)],
						                                                *doubleElement->vertex[sn333(i,j,k+1)], *doubleElement->vertex[sn333(i+1,j,k+1)],*doubleElement->vertex[sn333(i,j+1,k+1)],*doubleElement->vertex[sn333(i+1,j+1,k+1)],
																		mx, mX, my, mY, mz, mZ, mxy, mXy, mxY, mXY, myz, mYz, myZ, mYZ, mzx, mZx, mzX, mZX );
					CGradationBody::registerJoint( element[sn(level-1,x2+i,y2+j,z2+k)], x2+i, y2+j, z2+k, u );
					element[sn(level-1,x2+i,y2+j,z2+k)]->setElasticForce( youngsModulus, dampingCoefficient );
					doubleElement->deleteVertex( sn333(i2,j2,k2) );
				}
			}
		}
	}}}
}

bool CDestructionBody2::reconstructPhase2( CElement2 *doubleElement, int level, int x, int y, int n, Double youngsModulus, Double dampingCoefficient ){
	int x2=x*2;
	int y2=y*2;

	if( doubleElement->vertex[sn33(0,0)]==NULL && doubleElement->vertex[sn33(2,0)]==NULL && doubleElement->vertex[sn33(0,2)]==NULL && doubleElement->vertex[sn33(2,2)]==NULL ) return true;

	for(int i=0; i<2; i++){
		if( doubleElement->vertex[sn33(1,i*2)] ){
			if( doubleElement->vertex[sn33(0,i*2)]==NULL && doubleElement->vertex[sn33(2,i*2)]==NULL
			|| doubleElement->vertex[sn33(0,i*2)] && doubleElement->vertex[sn33(2,i*2)] && element[sn(level-1,x2,y2-1+i*3)]==NULL && element[sn(level-1,x2+1,y2-1+i*3)]==NULL ){
				doubleElement->deleteVertex( sn33(1,i*2) );
			}
		}else{
			if( doubleElement->vertex[sn33(0,i*2)] && !doubleElement->vertex[sn33(2,i*2)] || !doubleElement->vertex[sn33(0,i*2)] && doubleElement->vertex[sn33(2,i*2)]
			|| doubleElement->vertex[sn33(0,i*2)] && doubleElement->vertex[sn33(2,i*2)] && (element[sn(level-1,x2,y2-1+i*3)] || element[sn(level-1,x2+1,y2-1+i*3)]) ){
				doubleElement->vertex[sn33(1,i*2)] = new CVertex2( *doubleElement->vertex[sn33(0,i*2)], *doubleElement->vertex[sn33(2,i*2)] );
			}
		}
		if( doubleElement->vertex[sn33(i*2,1)] ){
			if( doubleElement->vertex[sn33(i*2,0)]==NULL && doubleElement->vertex[sn33(i*2,2)]==NULL
			|| doubleElement->vertex[sn33(i*2,0)] && doubleElement->vertex[sn33(i*2,2)] && element[sn(level-1,x2-1+i*3,y2)]==NULL && element[sn(level-1,x2-1+i*3,y2+1)]==NULL ){
				doubleElement->deleteVertex( sn33(i*2,1) );
			}
		}else{
			if( doubleElement->vertex[sn33(i*2,0)] && !doubleElement->vertex[sn33(i*2,2)] || !doubleElement->vertex[sn33(i*2,0)] && doubleElement->vertex[sn33(i*2,2)]
			|| doubleElement->vertex[sn33(i*2,0)] && doubleElement->vertex[sn33(i*2,2)] && (element[sn(level-1,x2-1+i*3,y2)] || element[sn(level-1,x2-1+i*3,y2+1)]) ){
				doubleElement->vertex[sn33(i*2,1)] = new CVertex2( *doubleElement->vertex[sn33(i*2,0)], *doubleElement->vertex[sn33(i*2,2)] );
			}
		}
	}
	if( doubleElement->vertex[sn33(0,0)] && doubleElement->vertex[sn33(2,0)] && doubleElement->vertex[sn33(0,2)] && doubleElement->vertex[sn33(2,2)] && doubleElement->vertex[sn33(1,1)] ){
		doubleElement->deleteVertex( sn33(1,1) );
	}

	doubleElement->reconstruct();
	CGradationBody2::registerJoint( doubleElement, x, y, n );
	return false;
}
bool CDestructionBody::reconstructPhase2( CElement *doubleElement, int level, int x, int y, int z, int n, Double youngsModulus, Double dampingCoefficient ){
	int x2=x*2;
	int y2=y*2;
	int z2=z*2;

	if( doubleElement->vertex[sn333(0,0,0)]==NULL && doubleElement->vertex[sn333(2,0,0)]==NULL && doubleElement->vertex[sn333(0,2,0)]==NULL && doubleElement->vertex[sn333(0,0,2)]==NULL
	&& doubleElement->vertex[sn333(2,2,0)]==NULL && doubleElement->vertex[sn333(2,0,2)]==NULL && doubleElement->vertex[sn333(0,2,2)]==NULL && doubleElement->vertex[sn333(2,2,2)]==NULL ) return true;

	for(int i=0; i<4; i++){
		int a=i%2*2;
		int b=i/2*2;
		int c=i%2*3-1;
		int d=i/2*3-1;
		if( doubleElement->vertex[sn333(1,a,b)] ){
			if( doubleElement->vertex[sn333(0,a,b)]==NULL && doubleElement->vertex[sn333(2,a,b)]==NULL
			|| doubleElement->vertex[sn333(0,a,b)] && doubleElement->vertex[sn333(2,a,b)] && element[sn(level-1,x2,y2+c,z2+i/2)]==NULL && element[sn(level-1,x2+1,y2+c,z2+i/2)]==NULL && element[sn(level-1,x2,y2+i%2,z2+d)]==NULL && element[sn(level-1,x2+1,y2+i%2,z2+d)]==NULL ){
				doubleElement->deleteVertex( sn333(1,a,b) );
			}
		}
		if( doubleElement->vertex[sn333(b,1,a)] ){
			if( doubleElement->vertex[sn333(b,0,a)]==NULL && doubleElement->vertex[sn333(b,2,a)]==NULL
			|| doubleElement->vertex[sn333(b,0,a)] && doubleElement->vertex[sn333(b,2,a)] && element[sn(level-1,x2+i/2,y2,z2+c)]==NULL && element[sn(level-1,x2+i/2,y2+1,z2+c)]==NULL && element[sn(level-1,x2+d,y2,z2+i%2)]==NULL && element[sn(level-1,x2+d,y2+1,z2+i%2)]==NULL ){
				doubleElement->deleteVertex( sn333(b,1,a) );
			}
		}
		if( doubleElement->vertex[sn333(a,b,1)] ){
			if( doubleElement->vertex[sn333(a,b,0)]==NULL && doubleElement->vertex[sn333(a,b,2)]==NULL
			|| doubleElement->vertex[sn333(a,b,0)] && doubleElement->vertex[sn333(a,b,2)] && element[sn(level-1,x2+c,y2+i/2,z2)]==NULL && element[sn(level-1,x2+c,y2+i/2,z2+1)]==NULL && element[sn(level-1,x2+i%2,y2+d,z2)]==NULL && element[sn(level-1,x2+i%2,y2+d,z2+1)]==NULL ){
				doubleElement->deleteVertex( sn333(a,b,1) );
			}
		}
	}

	for(int j=0; j<3; j+=2){
		int a = j/2*3-1;
		if( doubleElement->vertex[sn333(j,1,1)] ){
			if( doubleElement->vertex[sn333(j,0,0)]==NULL && doubleElement->vertex[sn333(j,2,0)]==NULL && doubleElement->vertex[sn333(j,0,2)]==NULL && doubleElement->vertex[sn333(j,2,2)]==NULL
			|| doubleElement->vertex[sn333(j,0,0)] && doubleElement->vertex[sn333(j,2,0)] && doubleElement->vertex[sn333(j,0,2)] && doubleElement->vertex[sn333(j,2,2)] && element[sn(level-1,x2+a,y2,z2)]==NULL && element[sn(level-1,x2+a,y2+1,z2)]==NULL && element[sn(level-1,x2+a,y2,z2+1)]==NULL && element[sn(level-1,x2+a,y2+1,z2+1)]==NULL ){
				doubleElement->deleteVertex( sn333(j,1,1) );
			}
		}
		if( doubleElement->vertex[sn333(1,j,1)] ){
			if( doubleElement->vertex[sn333(0,j,0)]==NULL && doubleElement->vertex[sn333(0,j,2)]==NULL && doubleElement->vertex[sn333(2,j,0)]==NULL && doubleElement->vertex[sn333(2,j,2)]==NULL
			|| doubleElement->vertex[sn333(0,j,0)] && doubleElement->vertex[sn333(0,j,2)] && doubleElement->vertex[sn333(2,j,0)] && doubleElement->vertex[sn333(2,j,2)] && element[sn(level-1,x2,y2+a,z2)]==NULL && element[sn(level-1,x2,y2+a,z2+1)]==NULL && element[sn(level-1,x2+1,y2+a,z2)]==NULL && element[sn(level-1,x2+1,y2+a,z2+1)]==NULL ){
				doubleElement->deleteVertex( sn333(1,j,1) );
			}
		}
		if( doubleElement->vertex[sn333(1,1,j)] ){
			if( doubleElement->vertex[sn333(0,0,j)]==NULL && doubleElement->vertex[sn333(2,0,j)]==NULL && doubleElement->vertex[sn333(0,2,j)]==NULL && doubleElement->vertex[sn333(2,2,j)]==NULL
			|| doubleElement->vertex[sn333(0,0,j)] && doubleElement->vertex[sn333(2,0,j)] && doubleElement->vertex[sn333(0,2,j)] && doubleElement->vertex[sn333(2,2,j)] && element[sn(level-1,x2,y2,z2+a)]==NULL && element[sn(level-1,x2+1,y2,z2+a)]==NULL && element[sn(level-1,x2,y2+1,z2+a)]==NULL && element[sn(level-1,x2+1,y2+1,z2+a)]==NULL ){
				doubleElement->deleteVertex( sn333(1,1,j) );
			}
		}
	}
	if( doubleElement->vertex[sn333(0,0,0)] && doubleElement->vertex[sn333(2,0,0)] && doubleElement->vertex[sn333(0,2,0)] && doubleElement->vertex[sn333(0,0,2)] 
	&& doubleElement->vertex[sn333(2,2,0)] && doubleElement->vertex[sn333(2,0,2)] && doubleElement->vertex[sn333(0,2,2)] && doubleElement->vertex[sn333(2,2,2)] ){
		doubleElement->deleteVertex( sn333(1,1,1) );
	}

	doubleElement->reconstruct();
	CGradationBody::registerJoint( doubleElement, x, y, z, n );
	return false;
}
//-------------------------------------------------------------------------------------------------------------------
class CDestructionBodyF2: public CDestructionBody2{
public:

	CDestructionBodyF2( CVoxel2 &voxel, Double centerX, Double centerY, Double angle, Double scale = 1.0 ){ this->construct( voxel, centerX, centerY, angle, scale ); }
	void renew( Double gravityConstant, Double youngsModulus, Double dampingCoefficient, Double timeStep, Double threshold = -1.0 );
	void disconnect( Double threshold );
	void reconstruct( Double youngsModulus, Double dampingCoefficient );
};

void CDestructionBodyF2::renew( Double gravityConstant, Double youngsModulus, Double dampingCoefficient, Double timeStep, Double threshold ){//2.375->0.117
	this->setElasticForce( youngsModulus, dampingCoefficient );//0.906->0.055
	if( threshold > 0.0 ){
		this->disconnect( threshold );//0.225->0.016
		this->reconstruct( youngsModulus, dampingCoefficient );//0.234->0.02
	}
	this->adjustAcceleration();//0.18->0.016
	this->addGravityAndRenew( gravityConstant, timeStep );//0.344->0.015
}


void CDestructionBodyF2::disconnect( Double threshold ){
	Double thresholdXthreshold = threshold*threshold;
	for( int i=0; i<(xSize+1)*(ySize+1); i++ ){
		if( joint[i] ){
			int x= i%(xSize+1);
			int y= i/(xSize+1);
			if( ( x>0 && y>0 && element[elementMarks[0] + x-1 + xSizes[0]*(y-1)] )//0.20->0.19
			|| ( x>0 && y<ySize && element[elementMarks[0] + x-1 + xSizes[0]*(y)] )
			|| ( x<xSize && y>0 && element[elementMarks[0] + x + xSizes[0]*(y-1)] )
			|| ( x<xSize && y<ySize && element[elementMarks[0] + x + xSizes[0]*(y)] ) ){
				if( joint[i]->connect[0] && joint[i]->vertex[0] && joint[i]->vertex[1]){
					CVertex2 *v1 = joint[i]->vertex[0];
					CVertex2 *v2 = joint[i]->vertex[1];
					Double x = v2->acceleration.x - v1->acceleration.x;
					Double y = v2->acceleration.y - v1->acceleration.y;
					if( x*x+y*y > thresholdXthreshold ){
						joint[i]->connect[0] = false;
						joint[i]->disconnected = true;
					}
				}
				if( joint[i]->connect[1] && joint[i]->vertex[1] && joint[i]->vertex[3]){
					CVertex2 *v1 = joint[i]->vertex[1];
					CVertex2 *v2 = joint[i]->vertex[3];
					Double x = v2->acceleration.x - v1->acceleration.x;
					Double y = v2->acceleration.y - v1->acceleration.y;
					if( x*x+y*y > thresholdXthreshold ){
						joint[i]->connect[1] = false;
						joint[i]->disconnected = true;
					}
				}
				if( joint[i]->connect[2] && joint[i]->vertex[3] && joint[i]->vertex[2]){
					CVertex2 *v1 = joint[i]->vertex[3];
					CVertex2 *v2 = joint[i]->vertex[2];
					Double x = v2->acceleration.x - v1->acceleration.x;
					Double y = v2->acceleration.y - v1->acceleration.y;
					if( x*x+y*y > thresholdXthreshold ){
						joint[i]->connect[2] = false;
						joint[i]->disconnected = true;
					}
				}
				if( joint[i]->connect[3] && joint[i]->vertex[2] && joint[i]->vertex[0]){
					CVertex2 *v1 = joint[i]->vertex[2];
					CVertex2 *v2 = joint[i]->vertex[0];
					Double x = v2->acceleration.x - v1->acceleration.x;
					Double y = v2->acceleration.y - v1->acceleration.y;
					if( x*x+y*y > thresholdXthreshold ){
						joint[i]->connect[3] = false;
						joint[i]->disconnected = true;
					}
				}
			}
		}
	}
}

void CDestructionBodyF2::reconstruct( Double youngsModulus, Double dampingCoefficient ){
	int n = 1;
	for( int level=1; level<levelSize; level++ ){
		n *= 2;
		int u = n/2;

		int y;
		for( y=0; y<ySizes[level]; y++ ){
		for( int x=0; x<xSizes[level]; x++ ){
			if( element[elementMarks[level] + x+xSizes[level]*y] ){
				int x2=x*2;
				int y2=y*2;
				int x4=x*4;
				int y4=y*4;

				CElement2 *e = element[elementMarks[level] + x+xSizes[level]*y];
				bool reconstructed = false;
				int j;
				for( j=0; j<2; j++ ){
				for( int i=0; i<2; i++ ){
					int i2 = i==0 ? 0 : 2;
					int j2 = j==0 ? 0 : 2;
					int i4 = i==0 ? 0 : 4;
					int j4 = j==0 ? 0 : 4;
					if( e->vertex[i2+j2*3] ){
						CVertex2 *v=e->vertex[i2+j2*3];
						if( level == 1 ){
							if( joint[x2+i2+(xSize+1)*(y2+j2)] && joint[x2+i2+(xSize+1)*(y2+j2)]->disconnected// inside of if statement was not accelerated because it occurs with low frequency.
							||  joint[x2+ 1+(xSize+1)*(y2+j2)] && joint[x2+ 1+(xSize+1)*(y2+j2)]->disconnected
							||  joint[x2+i2+(xSize+1)*(y2+ 1)] && joint[x2+i2+(xSize+1)*(y2+ 1)]->disconnected ){
								reconstructed = true;
								//break;
							}
						}else{
							bool mx = element[elementMarks[level-2] + x4+i2-1+xSizes[level-2]*(y4+j2)] || element[elementMarks[level-2] + x4+i2-1+xSizes[level-2]*(y4+j2+1)];
							bool mX = element[elementMarks[level-2] + x4+i2+2+xSizes[level-2]*(y4+j2)] || element[elementMarks[level-2] + x4+i2+2+xSizes[level-2]*(y4+j2+1)];
							bool my = element[elementMarks[level-2] + x4+i2  +xSizes[level-2]*(y4+j2-1)] || element[elementMarks[level-2] + x4+i2+1+xSizes[level-2]*(y4+j2-1)];
							bool mY = element[elementMarks[level-2] + x4+i2  +xSizes[level-2]*(y4+j2+2)] || element[elementMarks[level-2] + x4+i2+1+xSizes[level-2]*(y4+j2+2)];
							if( mx || mX || my || mY ){// inside of if statement was not accelerated because it occurs with low frequency.
								reconstructed = true;
								//break;
							}
						}
					}
				}}

				if( reconstructed ){
					if( e->vertex[0] && e->vertex[2] && !e->vertex[1] ) e->vertex[1] = new CVertex2( *e->vertex[0], *e->vertex[2] );
					if( e->vertex[6] && e->vertex[8] && !e->vertex[7] ) e->vertex[7] = new CVertex2( *e->vertex[6], *e->vertex[8] );
					if( e->vertex[0] && e->vertex[6] && !e->vertex[3] ) e->vertex[3] = new CVertex2( *e->vertex[0], *e->vertex[6] );
					if( e->vertex[2] && e->vertex[8] && !e->vertex[5] ) e->vertex[5] = new CVertex2( *e->vertex[2], *e->vertex[8] );
					if( !e->vertex[4] ) e->vertex[4] = new CVertex2( *e->vertex[1], *e->vertex[3], *e->vertex[5], *e->vertex[7] );

					for( int j=0; j<2; j++ ){
					for( int i=0; i<2; i++ ){
						int i2 = i==0 ? 0 : 2;
						int j2 = j==0 ? 0 : 2;
						int i4 = i==0 ? 0 : 4;
						int j4 = j==0 ? 0 : 4;
						if( e->vertex[i2+j2*3] ){
							CVertex2 *v=e->vertex[i2+j2*3];
							if( level == 1 ){
								if( joint[x2+i2+(xSize+1)*(y2+j2)] && joint[x2+i2+(xSize+1)*(y2+j2)]->disconnected// inside of if statement was not accelerated because it occurs with low frequency.
								||  joint[x2+ 1+(xSize+1)*(y2+j2)] && joint[x2+ 1+(xSize+1)*(y2+j2)]->disconnected
								||  joint[x2+i2+(xSize+1)*(y2+ 1)] && joint[x2+i2+(xSize+1)*(y2+ 1)]->disconnected ){
									CElement2 *e1=element[sn(level-1,x2+i,y2+j)] = new CElement2( *e->vertex[sn33(i,j)], *e->vertex[sn33(i+1,j)],*e->vertex[sn33(i,j+1)],*e->vertex[sn33(i+1,j+1)] );
									CVoxelBody2::registerJoint( e1, x2+i, y2+j );
									e1->setElasticForce( youngsModulus, dampingCoefficient );
									delete v; 
									e->vertex[i2+j2*3] = NULL;
								}
							}else{
								bool mx = element[elementMarks[level-2] + x4+i2-1+xSizes[level-2]*(y4+j2)] || element[elementMarks[level-2] + x4+i2-1+xSizes[level-2]*(y4+j2+1)];
								bool mX = element[elementMarks[level-2] + x4+i2+2+xSizes[level-2]*(y4+j2)] || element[elementMarks[level-2] + x4+i2+2+xSizes[level-2]*(y4+j2+1)];
								bool my = element[elementMarks[level-2] + x4+i2  +xSizes[level-2]*(y4+j2-1)] || element[elementMarks[level-2] + x4+i2+1+xSizes[level-2]*(y4+j2-1)];
								bool mY = element[elementMarks[level-2] + x4+i2  +xSizes[level-2]*(y4+j2+2)] || element[elementMarks[level-2] + x4+i2+1+xSizes[level-2]*(y4+j2+2)];
								if( mx || mX || my || mY ){// inside of if statement was not accelerated because it occurs with low frequency.
									element[sn(level-1,x2+i,y2+j)] = new CElement2( *e->vertex[sn33(i,j)], *e->vertex[sn33(i+1,j)],*e->vertex[sn33(i,j+1)],*e->vertex[sn33(i+1,j+1)], mx, mX, my, mY );
									element[sn(level-1,x2+i,y2+j)]->setElasticForce( youngsModulus, dampingCoefficient );
									CGradationBody2::registerJoint( element[sn(level-1,x2+i,y2+j)], x2+i, y2+j, u );
									e->deleteVertex( sn33(i2,j2) );
								}
							}
						}
					}}
				}
			}
		}}
		for( y=0; y<ySizes[level]; y++ ){
		for( int x=0; x<xSizes[level]; x++ ){
			if( element[elementMarks[level] + x+xSizes[level]*y] ){
				CElement2 *e = element[elementMarks[level] + x+xSizes[level]*y];
				if( !e->vertex[0] && !e->vertex[2] && !e->vertex[6] && !e->vertex[8] ){
					delete e;
					element[elementMarks[level] + x+xSizes[level]*y] = NULL;
				}else{
					int x2=x*2;
					int y2=y*2;

					if( e->vertex[1] ){
						if( !e->vertex[0] && !e->vertex[2]
						|| e->vertex[0] && e->vertex[2] && !element[elementMarks[level-1] + x2+xSizes[level-1]*(y2-1)] && !element[elementMarks[level-1] + x2+1+xSizes[level-1]*(y2-1)] ){
							delete e->vertex[1];
							e->vertex[1] = NULL;
						}
					}else{
						if( e->vertex[0] && !e->vertex[2] || !e->vertex[0] && e->vertex[2]
						|| e->vertex[0] && e->vertex[2] && (element[elementMarks[level-1] + x2+xSizes[level-1]*(y2-1)] || element[elementMarks[level-1] + x2+1+xSizes[level-1]*(y2-1)]) ){
							e->vertex[1] = new CVertex2( *e->vertex[0], *e->vertex[2] );
						}
					}
					if( e->vertex[7] ){
						if( !e->vertex[6] && !e->vertex[8]
						|| e->vertex[6] && e->vertex[8] && !element[elementMarks[level-1] + x2+xSizes[level-1]*(y2+2)] && !element[elementMarks[level-1] + x2+1+xSizes[level-1]*(y2+2)] ){
							delete e->vertex[7];
							e->vertex[7] = NULL;
						}
					}else{
						if( e->vertex[6] && !e->vertex[8] || !e->vertex[6] && e->vertex[8]
						|| e->vertex[6] && e->vertex[8] && (element[elementMarks[level-1] + x2+xSizes[level-1]*(y2+2)] || element[elementMarks[level-1] + x2+1+xSizes[level-1]*(y2+2)]) ){
							e->vertex[7] = new CVertex2( *e->vertex[6], *e->vertex[8] );
						}
					}
					if( e->vertex[3] ){
						if( !e->vertex[0] && !e->vertex[6]
						|| e->vertex[0] && e->vertex[6] && !element[elementMarks[level-1] + x2-1+xSizes[level-1]*y2] && !element[elementMarks[level-1] + x2-1+xSizes[level-1]*(y2+1)] ){
								delete e->vertex[3];
							e->vertex[3] = NULL;
						}
					}else{
						if( e->vertex[0] && !e->vertex[6] || !e->vertex[0] && e->vertex[6]
						|| e->vertex[0] && e->vertex[6] && (element[elementMarks[level-1] + x2-1+xSizes[level-1]*y2] || element[elementMarks[level-1] + x2-1+xSizes[level-1]*(y2+1)]) ){
							e->vertex[3] = new CVertex2( *e->vertex[0], *e->vertex[6] );
						}
					}
					if( e->vertex[5] ){
						if( !e->vertex[2] && !e->vertex[8]
						|| e->vertex[2] && e->vertex[8] && !element[elementMarks[level-1] + x2+2+xSizes[level-1]*y2] && !element[elementMarks[level-1] + x2+2+xSizes[level-1]*(y2+1)] ){
							delete e->vertex[5];
							e->vertex[5] = NULL;
						}
					}else{
						if( e->vertex[2] && !e->vertex[8] || !e->vertex[2] && e->vertex[8]
						|| e->vertex[2] && e->vertex[8] && (element[elementMarks[level-1] + x2+2+xSizes[level-1]*y2] || element[elementMarks[level-1] + x2+2+xSizes[level-1]*(y2+1)]) ){
							e->vertex[5] = new CVertex2( *e->vertex[2], *e->vertex[8] );
						}
					}
					if( e->vertex[0] && e->vertex[2] && e->vertex[6] && e->vertex[8] && e->vertex[4] ){
							delete e->vertex[4];
							e->vertex[4] = NULL;
					}
					int vertexSize = 0;
					if( e->vertex[0] ) e->vertexList[vertexSize++] = e->vertex[0];
					else if(e->vertex[6]) e->vertexList[vertexSize++] = e->vertex[4];
					if( e->vertex[1] ) e->vertexList[vertexSize++] = e->vertex[1];
					if( e->vertex[2] ) e->vertexList[vertexSize++] = e->vertex[2];
					else if(e->vertex[0]) e->vertexList[vertexSize++] = e->vertex[4];
					if( e->vertex[5] ) e->vertexList[vertexSize++] = e->vertex[5];
					if( e->vertex[8] ) e->vertexList[vertexSize++] = e->vertex[8];
					else if(e->vertex[2]) e->vertexList[vertexSize++] = e->vertex[4];
					if( e->vertex[7] ) e->vertexList[vertexSize++] = e->vertex[7];
					if( e->vertex[6] ) e->vertexList[vertexSize++] = e->vertex[6];
					else if(e->vertex[8]) e->vertexList[vertexSize++] = e->vertex[4];
					if( e->vertex[3] ) e->vertexList[vertexSize++] = e->vertex[3];
					e->vertexSize = vertexSize;

					Double cx = e->vertexList[0]->positionOrigin.x;
					Double cy = e->vertexList[0]->positionOrigin.y;
					int i, j;
					for( i=1; i<vertexSize; i++ ){
						cx += e->vertexList[i]->positionOrigin.x;
						cy += e->vertexList[i]->positionOrigin.y;
					}
					Double rVertexSize = vertexSize < 9 ?
						vertexSize < 5 ? vertexSize<3?(vertexSize==1?1.0:0.5):(vertexSize==3?0.333333333:0.25) : vertexSize<7?(vertexSize==5?0.2:0.1666666666):(vertexSize==7?0.142857142857:0.125)
						: vertexSize < 13 ? vertexSize<11?(vertexSize==9?0.111111111:0.1):(vertexSize==11?0.090909090909:0.0833333333333) : vertexSize<15?(vertexSize==13?0.076923076923:0.0714285714285):(vertexSize==15?0.066666666666:0.0625);
					cx *= rVertexSize;
					cy *= rVertexSize;
					//cx /= vertexSize;
					//cy /= vertexSize;
					Double averageRxR = 0.0;
					for( i=0; i<vertexSize; i++ ){
						e->vertexList[i]->R.x = e->vertexList[i]->positionOrigin.x - cx;
						e->vertexList[i]->R.y = e->vertexList[i]->positionOrigin.y - cy;
						averageRxR += e->vertexList[i]->R.x*e->vertexList[i]->R.x+e->vertexList[i]->R.y*e->vertexList[i]->R.y;
					}
					averageRxR *= rVertexSize;
					e->density = 1.0 / averageRxR;

					//e->setElasticForce( youngsModulus, dampingCoefficient );
					//CGradationBody2::registerJoint( e, x, y, n );
					int xo=x*n;
					int yo=y*n;
					for( int j=0; j<2; j++ ){
					for( int i=0; i<2; i++ ){
						int i2 = i==0 ? 0 : 2;
						int j2 = j==0 ? 0 : 2;
						if( e->vertex[i2+j2*3] ){
							int k=xo+u*i2+(xSize+1)*(yo+u*j2);
							if( !joint[k] ) joint[k] = new CJoint2();
							joint[k]->vertex[3-i-j*2] = e->vertex[i2+j2*3];
							if(e->vertex[1+j2*3] ){
								k=xo+u+(xSize+1)*(yo+u*j2);
								if( !joint[k] ) joint[k] = new CJoint2();
								joint[k]->vertex[i+(1-j)*2] = e->vertex[1+j2*3];
							}
							if(e->vertex[i2+3] ){
								k=xo+u*i2+(xSize+1)*(yo+u);
								if( !joint[k] ) joint[k] = new CJoint2();
								joint[k]->vertex[1-i+j*2] = e->vertex[i2+3];
							}
							if( e->vertex[4] ){
								k=xo+u+(xSize+1)*(yo+u);
								if( !joint[k] ) joint[k] = new CJoint2();
								joint[k]->vertex[i+j*2] = e->vertex[4];
							}
						}
					}}

				}
			}
		}}
	}
}