clothsim.java

实现Java 3D环境下基于质点－弹簧模型的柔体模拟

      {
     normals[((xTotal-1)*2*(yTotal-2)+2*(i-1))*3*3+12]= normalsAverage.x;
      normals[((xTotal-1)*2*(yTotal-2)+2*(i-1))*3*3+13]=normalsAverage.y;
      normals[((xTotal-1)*2*(yTotal-2)+2*(i-1))*3*3+14]=normalsAverage.z; 	
      }
       }
       
       
       }//优化三角面顶点法向量,以使得三角面平滑过渡
       //顶点法向量坐标求平均值
         public float NormalsAverage(float nromalsa,float normalsCout)
        {
         float normalsb=0.0f;
        	normalsb=nromalsa/normalsCout;
        	return normalsb;
        } 
        //顶点法向量坐标相加 
public float NormalsAdd(float nromals1,float normals2)
        {
         float normalsb=0.0f;
         normalsb=	nromals1+normals2;
        	return normalsb;
        }    
       }
		

       
///////////////////////////////////////////////////////////
//              物理模型Physics内部类                    //
//                                                       //
//     创建一个质点位置更新模型,根据当前质点的位置,      //
//     计算质点的内力,外力,然后利用简单的欧拉公式        //
//     计算质点的加速度,速度,质点的位移,计算下个         //
//     时间点的位置                                      //
//                                                       //
//                                                       //
///////////////////////////////////////////////////////////	
		
	//物理模型内部类	
		
		class Physics
{
 
 



//根据内部弹簧位移计算内力，并结合风力，空气阻力计算质点的总受力
 public void calculateF(Point3d pointPosition1[][],Vector3d vectorF1[][],Vector3d vector3dV1[][])
 {
 
 


for(int i=0;i<xTotal;i++)
for(int j=0;j<yTotal;j++)
{
	vectorF1[i][j]=add(new Vector3d(0.0,0.0,0.0),new Vector3d(0.0,0.0,0.0));
//初始化质点力矢量为零矢量
if(j<(yTotal-1))vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i][j+1],lstretch));		
if(j>0)vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i][j-1],lstretch));
if(i!=(xTotal-1))vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i+1][j],lstretch));
if(i!=0)vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i-1][j],lstretch));
//计算质点由于结构弹簧所受的力
if(i!=(xTotal-1)&j!=0)vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i+1][j-1],lshear));
if(i!=(xTotal-1)&j!=(yTotal-1))vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i+1][j+1],lshear));
if(i!=0&j!=0)vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i-1][j-1],lshear));
if(i!=0&j!=(yTotal-1))vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i-1][j+1],lshear));
//计算质点由于剪切弹簧所受的力
if(i>1)vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i-2][j],lbend));
if(i<(xTotal-2))vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i+2][j],lbend));
if(j>1)vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i][j-2],lbend));
if(j<(yTotal-2))vectorF1[i][j]=add(vectorF1[i][j],neighborF(pointPosition1[i][j],pointPosition1[i][j+2],lbend));
//计算质点由于弯曲弹簧所受的力

 vectorF1[i][j]=add(vectorF1[i][j],new Vector3d(0.0,-mg,0.0));
 
//计算重力


Vector3d vectorF1wind=new Vector3d(Fwindx*(1.0-Math.random()*0.8),-0.1*Math.random(),Fwindz*(1.0-Math.random()*0.8));
vectorF1[i][j]=add(vectorF1[i][j],vectorF1wind);


Vector3d vectorFwind=new Vector3d(500.0,0.0,-250.0);

Vector3d a=new Vector3d(vectorFwind.x-vector3dV1[i][j].x,vectorFwind.y-vector3dV1[i][j].y,vectorFwind.z-vector3dV1[i][j].z);
Vector3d b=new Vector3d(0.0,0.0,0.0);
b=pointNormal[i][j];
if(b.length()!=0)

b.normalize() ;
b.dot(a);

vectorFwind.x=-kvi*b.length()*pointNormal[i][j].x;
vectorFwind.y=-kvi*b.length()*pointNormal[i][j].y;
vectorFwind.z=-kvi*b.length()*pointNormal[i][j].z;

//vectorF1[i][j]=add(vectorF1[i][j],vectorFwind);

//计算风力
Vector3d vectorF1air=new Vector3d(vector3dV1[i][j]);
if(vector3dV1[i][j].length()!=0.0)
vectorF1air.normalize();
vectorF1air.x=-kair*vector3dV1[i][j].length()*vector3dV1[i][j].length()*vectorF1air.x;
vectorF1air.y=-kair*vector3dV1[i][j].length()*vector3dV1[i][j].length()*vectorF1air.y;
vectorF1air.z=-kair*vector3dV1[i][j].length()*vector3dV1[i][j].length()*vectorF1air.z;
vectorF1[i][j]=add(vectorF1[i][j],vectorF1air);
//计算空气阻力
vectorF1[i][j].x=vectorF1[i][j].x/m;
vectorF1[i][j].y=vectorF1[i][j].y/m;
vectorF1[i][j].z=vectorF1[i][j].z/m;
//计算质点加速度

	}

	 

	}
	
	
	
	//矢量相加函数
	public Vector3d add(Vector3d vector1,Vector3d vector2)
	{
		Vector3d vector3d0=new Vector3d();
	vector3d0.x=(vector1.x+vector2.x);
vector3d0.y=(vector1.y+vector2.y);
vector3d0.z=(vector1.z+vector2.z);
return vector3d0;
	}
	


 
//	根据质点加速度，计算质点速度
public void calculateV(Vector3d vector3dV1[][])
{
for(int i=0;i<xTotal;i++)
for(int j=0;j<yTotal;j++)	
{
	if(i==0&downFlage==false)
	{
		}
else 
	{
			
			
if(((pointPosition[i][j].x<50.0&pointPosition[i][j].x>-50))&(pointPosition[i][j].z<50&pointPosition[i][j].z>-50)&pointPosition[i][j].y<-0.990)	
	{vector3dV1[i][j].x=vector3dV[i][j].x+vectorF[i][j].x*t/2;
		vector3dV1[i][j].y=0.0;
			vector3dV1[i][j].z=vector3dV[i][j].z+vectorF[i][j].z*t/2;
		}
		else if(pointPosition[i][j].x>-0.015&pointPosition[i][j].x<0.015&pointPosition[i][j].y>-0.990)
		{
		vector3dV1[i][j].x=0.1;
		vector3dV1[i][j].y=vector3dV[i][j].y+vectorF[i][j].y*t/2;
			vector3dV1[i][j].z=0.0;	
		}
			else
		vector3dV1[i][j]=add(vector3dV[i][j],new Vector3d(vectorF[i][j].x*t/2,vectorF[i][j].y*t/2,vectorF[i][j].z*t/2));
}
}


}




//根据质点速度计算质点的位移
public void calculatePosition(Point3d pointPosition1[][],Vector3d vector3dV1[][])
{
for(int i=0;i<xTotal;i++)
for(int j=0;j<yTotal;j++)	
{
pointPosition1[i][j].x=pointPosition[i][j].x+vector3dV1[i][j].x*t;
pointPosition1[i][j].y=pointPosition[i][j].y+vector3dV1[i][j].y*t;
pointPosition1[i][j].z=pointPosition[i][j].z+vector3dV1[i][j].z*t;
}
}



//根据两质点间的距离，由虎克定律计算质点间的受力
	public Vector3d neighborF(Point3d point2,Point3d point1,double constanlength)
	{
		Vector3d vector3d=new Vector3d(point1.x-point2.x,point1.y-point2.y,point1.z-point2.z);
		Vector3d vector3dBack=new Vector3d(vector3d);
		vector3d.normalize();
		
		
		Vector3d vector3dF=new Vector3d(k*(vector3dBack.length()-constanlength)*vector3d.x,k*(vector3dBack.length()-constanlength)*vector3d.y,k*(vector3dBack.length()-constanlength)*vector3d.z);
		
		return vector3dF;
	
		}
		//计算两个向量的平均值
		public void VecotorAverage(Vector3d v1[][],Vector3d v2[][])
        {
         
         for(int i=0;i<xTotal;i++)
for(int j=0;j<yTotal;j++)
{
v1[i][j].x=(v1[i][j].x+v2[i][j].x)/2;	
v1[i][j].y=(v1[i][j].y+v2[i][j].y)/2;
v1[i][j].z=(v1[i][j].z+v2[i][j].z)/2;	
}
        
        } 
        
        //利用简单欧拉公式和改进欧拉公式混合计算刷新质点位置
        public void oulaCalculate()
{if(ouLaFlage%3!=2)
{	calculateF(pointPosition,vectorF,vector3dV);
	 calculateV(vector3dV);
	 calculatePosition(pointPosition,vector3dV);//利用简单欧拉公式刷新计算质点位置
	 ouLaFlage++;
	 if(ouLaFlage==3)ouLaFlage=0;
	 }
	
	else
	{ calculateF(pointPosition,vectorF,vector3dV);
	 calculateV(vectorVBack);
	 calculatePosition(pointPositionBack,vectorVBack);
	 calculateF(pointPositionBack,vectorABack,vectorVBack);
	 VecotorAverage(vectorF,vectorABack);
	  calculateV(vector3dV);
	  
	 calculatePosition(pointPosition,vector3dV);//利用改进欧拉公式刷新计算质点位置
	 ouLaFlage++;
	  if(ouLaFlage==3)ouLaFlage=0;}
}		
		//物理类构造函数
	public Physics()
{
}	
		
}
      
/////////////////////////////////////////////////////////////
//             clothSim键盘响应内部类                      //
//                                                         //
//     响应键盘的输入,调整风速,改变面模型形式              //
//                                                         //
//                                                         //
/////////////////////////////////////////////////////////////	

// 键盘响应行为类
public class ClothSimBehavior extends Behavior{
		          	
        // 构造函数
        ClothSimBehavior()
        {
        				
        }

        
		// initialize 方法
        public void initialize()
        {
			this.wakeupOn(new WakeupOnAWTEvent(KeyEvent.KEY_PRESSED));
					//设置唤醒条件为有键按下
        }

        // processStimulus方法
        public void processStimulus(Enumeration criteria)
        {
           	WakeupCriterion wakeup = null;
			AWTEvent[] event = null;
			// 对获得的事件进行解析，并作出响应的反应
			wakeup = (WakeupCriterion)criteria.nextElement( );
					// 获得当前的唤醒条件
			if ( wakeup instanceof WakeupOnAWTEvent )
			{
				event = ((WakeupOnAWTEvent)wakeup).getAWTEvent( );
						// 获得当前的唤醒事件
				KeyEvent keyevent = (KeyEvent)event[0];
						// 转化为键盘事件
				// 判断是否为0键，并做相应参数变化
				
				if(keyevent.getKeyCode()==KeyEvent.VK_0)
					{
					 k=500.0;//弹性模量常数
 kair=0.5;//空气阻力常数
 kvi=0.25;//空气流参数
Fwindx=0.0;
Fwindz=0.0;
tWait=2;}
if(keyevent.getKeyCode()==KeyEvent.VK_1)
					{
					 k=500.0;//弹性模量常数
 kair=0.018;//空气阻力常数
 kvi=0.25;//空气流参数
Fwindx=0.65;
Fwindz=-0.45;
tWait=2;}
					if(keyevent.getKeyCode()==KeyEvent.VK_2)
					{
					 k=500.0;//弹性模量常数
 kair=0.018;//空气阻力常数
 kvi=0.25;//空气流参数
Fwindx=2.5;
Fwindz=-1.5;
tWait=2;}

	if(keyevent.getKeyCode()==KeyEvent.VK_3)
	{
					 k=1500.0;//弹性模量常数
 kair=0.0025;//空气阻力常数
 kvi=0.25;//空气流参数

Fwindx=10.0;
Fwindz=-7.5;
tWait=3;				}
			
				
if(keyevent.getKeyCode()==KeyEvent.VK_ENTER)
{	
if(clothSimFaceType==1)
{s.setWhichChild(1);//显示线模型
sPoint.setWhichChild(0);//显示质点模型
tWait=0;//设置刷新线程刷新等待时间
clothSimFaceType--;

}

else 
{s.setWhichChild(0);//显示面模型 
sPoint.setWhichChild(Switch.CHILD_NONE);//不显示质点球模型
tWait=0;//设置刷新线程刷新等待时间
clothSimFaceType++;	
}
}


if(keyevent.getKeyCode()==KeyEvent.VK_S)
{if(timerFlagevcout==0)
	{timerFlage=false;
	timerFlagevcout++;}
	else
	{
	timerFlage=true;
	timerFlagevcout--;	}//刷新线程暂停	

	}
if(keyevent.getKeyCode()==KeyEvent.VK_D)
{

	downFlage=true;	
	}						
}			
            this.wakeupOn(new WakeupOnAWTEvent(KeyEvent.KEY_PRESSED));// 设置下一次的唤醒条件
        }
    }
    
}