Games103的作业4：水波

这个作业做得我头疼，奇奇怪怪的bug一堆，目前暂时没有弄水作用于木块，我猜应该是用浮力对木块作用的力矩来解决

下面记录下代码

using UnityEngine;
using System.Collections;

public class wave_motion : MonoBehaviour 
{
	int size 		= 100;
	float rate 		= 0.005f;
	float gamma		= 0.004f;
	float damping 	= 0.996f;
	float[,] 	old_h;
	float[,]	low_h;
	float[,]	vh;
	float[,]	b;

	bool [,]	cg_mask;
	float[,]	cg_p;
	float[,]	cg_r;
	float[,]	cg_Ap;
	bool 	tag=true;

	Vector3 	cube_v = Vector3.zero;
	Vector3 	cube_w = Vector3.zero;


	// Use this for initialization
	void Start () 
	{
		Mesh mesh = GetComponent<MeshFilter> ().mesh;
		mesh.Clear ();

		Vector3[] X=new Vector3[size*size];

		for (int i=0; i<size; i++)
		for (int j=0; j<size; j++) 
		{
			X[i*size+j].x=i*0.1f-size*0.05f;
			X[i*size+j].y=0;
			X[i*size+j].z=j*0.1f-size*0.05f;
		}

		int[] T = new int[(size - 1) * (size - 1) * 6];
		int index = 0;
		for (int i=0; i<size-1; i++)
		for (int j=0; j<size-1; j++)
		{
			T[index*6+0]=(i+0)*size+(j+0);
			T[index*6+1]=(i+0)*size+(j+1);
			T[index*6+2]=(i+1)*size+(j+1);
			T[index*6+3]=(i+0)*size+(j+0);
			T[index*6+4]=(i+1)*size+(j+1);
			T[index*6+5]=(i+1)*size+(j+0);
			index++;
		}
		mesh.vertices  = X;
		mesh.triangles = T;
		mesh.RecalculateNormals ();

		low_h 	= new float[size,size];
		old_h 	= new float[size,size];
		vh 	  	= new float[size,size];
		b 	  	= new float[size,size];

		cg_mask	= new bool [size,size];
		cg_p 	= new float[size,size];
		cg_r 	= new float[size,size];
		cg_Ap 	= new float[size,size];

		for (int i=0; i<size; i++)
		for (int j=0; j<size; j++) 
		{
			low_h[i,j]=99999;
			old_h[i,j]=0;
			vh[i,j]=0;
		}
	}

	void A_Times(bool[,] mask, float[,] x, float[,] Ax, int li, int ui, int lj, int uj)
	{
		for(int i=li; i<=ui; i++)
		for(int j=lj; j<=uj; j++)
		if(i>=0 && j>=0 && i<size && j<size && mask[i,j])
		{
			Ax[i,j]=0;
			if(i!=0)		Ax[i,j]-=x[i-1,j]-x[i,j];
			if(i!=size-1)	Ax[i,j]-=x[i+1,j]-x[i,j];
			if(j!=0)		Ax[i,j]-=x[i,j-1]-x[i,j];
			if(j!=size-1)	Ax[i,j]-=x[i,j+1]-x[i,j];
		}
	}

	float Dot(bool[,] mask, float[,] x, float[,] y, int li, int ui, int lj, int uj)
	{
		float ret=0;
		for(int i=li; i<=ui; i++)
		for(int j=lj; j<=uj; j++)
		if(i>=0 && j>=0 && i<size && j<size && mask[i,j])
		{
			ret+=x[i,j]*y[i,j];
		}
		return ret;
	}

	void Conjugate_Gradient(bool[,] mask, float[,] b, float[,] x, int li, int ui, int lj, int uj)
	{
		//Solve the Laplacian problem by CG.
		A_Times(mask, x, cg_r, li, ui, lj, uj);

		for(int i=li; i<=ui; i++)
		for(int j=lj; j<=uj; j++)
		if(i>=0 && j>=0 && i<size && j<size && mask[i,j])
		{
			cg_p[i,j]=cg_r[i,j]=b[i,j]-cg_r[i,j];
		}

		float rk_norm=Dot(mask, cg_r, cg_r, li, ui, lj, uj);

		for(int k=0; k<128; k++)
		{
			if(rk_norm<1e-10f)	break;
			A_Times(mask, cg_p, cg_Ap, li, ui, lj, uj);
			float alpha=rk_norm/Dot(mask, cg_p, cg_Ap, li, ui, lj, uj);

			for(int i=li; i<=ui; i++)
			for(int j=lj; j<=uj; j++)
			if(i>=0 && j>=0 && i<size && j<size && mask[i,j])
			{
				x[i,j]   +=alpha*cg_p[i,j];
				cg_r[i,j]-=alpha*cg_Ap[i,j];
			}

			float _rk_norm=Dot(mask, cg_r, cg_r, li, ui, lj, uj);
			float beta=_rk_norm/rk_norm;
			rk_norm=_rk_norm;

			for(int i=li; i<=ui; i++)
			for(int j=lj; j<=uj; j++)
			if(i>=0 && j>=0 && i<size && j<size && mask[i,j])
			{
				cg_p[i,j]=cg_r[i,j]+beta*cg_p[i,j];
			}
		}

	}

	void Shallow_Wave(float[,] old_h, float[,] h, float [,] new_h)
	{
		//Step 1:
		//TODO: Compute new_h based on the shallow wave model.
		for (int i = 0; i < size; i++)
			for (int j = 0; j < size; j++)
			{
				new_h[i, j] = h[i, j] + (h[i, j] - old_h[i, j]) * damping;
				if (i - 1 >= 0)
					new_h[i, j] += (h[i - 1, j] - h[i, j]) * rate;
				if (i + 1 < size)
					new_h[i, j] += (h[i + 1, j] - h[i, j]) * rate;
				if (j - 1 >= 0)
					new_h[i, j] += (h[i, j - 1] - h[i, j]) * rate;
				if (j + 1 < size)
					new_h[i, j] += (h[i, j + 1] - h[i, j]) * rate;
			}

		//Step 2: Block->Water coupling
		GameObject cube = GameObject.Find("Cube");
		//TODO: for block 1, calculate low_h.
		Vector3 cube_p = cube.transform.position;
		Mesh cube_m = cube.GetComponent<MeshFilter>().mesh;
		Bounds bound = cube_m.bounds;

		int li = (int)((cube_p.x) * 10 + 50 - 2);	//世界坐标里长度0.1对应一方格，水面大小为100X100
		int ui = (int)((cube_p.x) * 10 + 50 + 2);
		int lj = (int)((cube_p.z) * 10 + 50 - 2);
		int uj = (int)((cube_p.z) * 10 + 50 + 2);	

		for (int i = li ; i <= ui; i++)      
			for (int j = lj; j <= uj; j++)


				//for (int i = li-3; i <= ui+3; i++)		//我并不怎么理解这么-3 +3的作用，但他的确可以消除残留痕迹，而反之将+3-3这些等价变换到ui，li，uj，lj却不能做到消除残留痕迹的效果
			//for (int j = lj-3; j <= uj+3; j++)
				if (i >= 0 && j >= 0 && i < size && j < size)
				{
				Vector3 p = new Vector3(i * 0.1f - 5f, -11f, j * 0.1f - 5f);
				Vector3 q = new Vector3(i * 0.1f - 5f, -10f, j * 0.1f - 5f);
				
				p = cube.transform.InverseTransformPoint(p);
				q = cube.transform.InverseTransformPoint(q);

				Ray ray = new Ray(p, q - p);
				float dis= 99999;
				bound.IntersectRay(ray,out dis);
				low_h[i, j] = -11f+dis;
			}


		//TODO: then set up b and cg_mask for conjugate gradient.
		for(int i=0;i<size;i++)
			for(int j=0;j<size;j++)
            {
				if(low_h[i,j]>h[i,j])
                {
					b[i, j] = 0;
					cg_mask[i, j] = false;
					vh[i, j] = 0;
                }
				else
                {
					b[i, j] =  (new_h[i, j] - low_h[i, j])/rate;
					cg_mask[i, j] = true;
                }

            }


		//TODO: Solve the Poisson equation to obtain vh (virtual height).
		Conjugate_Gradient(cg_mask, b, vh, li - 1, ui + 1, lj - 1, uj + 1);

		//TODO: for block 2, calculate low_h.
		//TODO: then set up b and cg_mask for conjugate gradient.
		//TODO: Solve the Poisson equation to obtain vh (virtual height).

		//TODO: Diminish vh.
		for (int i = 0; i < size; i++)
			for (int j = 0; j < size; j++)	
			{
				if (cg_mask[i, j])
					vh[i, j] *= gamma;
			}
		//TODO: Update new_h by vh.
		for (int i = 0; i < size; i++)
			for (int j = 0; j < size; j++)
			{
				if ((i - 1) >= 0)
					new_h[i, j] += (vh[i - 1, j] - vh[i, j]) * rate ;
				if ((i + 1) < size)
					new_h[i, j] += (vh[i + 1, j] - vh[i, j]) * rate;
				if ((j - 1) >= 0)
					new_h[i, j ] += (vh[i, j - 1] - vh[i, j]) * rate ;
				if ((j + 1) < size)
					new_h[i, j ] += (vh[i, j + 1] - vh[i, j]) * rate;
			}
		//Step 3
		//TODO: old_h <- h; h <- new_h;
		for (int i=0;i<size;i++)
			for(int j=0;j<size;j++)
            {
				old_h[i, j] = h[i, j];
				h[i, j] = new_h[i, j];
            }
		//Step 4: Water->Block coupling.
		//More TODO here.
	}
	

	// Update is called once per frame
	void Update () 
	{
		Mesh mesh = GetComponent<MeshFilter> ().mesh;
		Vector3[] X    = mesh.vertices;
		float[,] new_h = new float[size, size];
		float[,] h     = new float[size, size];

		//TODO: Load X.y into h.
		for (int i = 0; i < size; i++)
			for (int j = 0; j < size; j++)
				h[i, j] = X[i * size + j].y;

		if (Input.GetKeyDown ("r")) 
		{
			//TODO: Add random water.
			h[15, 15] += 0.1f;
		}
	
		for(int l=0; l<8; l++)
		{
			Shallow_Wave(old_h, h, new_h);
		}

		//TODO: Store h back into X.y and recalculate normal.
		for (int i = 0; i < size; i++)
			for (int j = 0; j < size; j++)
				X[i*size+j].y = h[i, j];
		mesh.vertices = X;
		mesh.RecalculateNormals();
	}
}