📄 main.cpp
字号:
/* Copyright (C) Miguel Gomez, 2000. * All rights reserved worldwide. * * This software is provided "as is" without express or implied * warranties. You may freely copy and compile this source into * applications you distribute provided that the copyright text * below is included in the resulting source code, for example: * "Portions Copyright (C) Miguel Gomez, 2000" */#include <stdio.h>#include <stdlib.h>#include <time.h>#include <cmath>#include <GL/glut.h>#include "body.h"#include "mathphys.h"// TEMP HACK: x is right, y is forward, z is up//// 3D Camera//struct Camera : public COORD_FRAME{ float D, F; //near and far plane VECTOR3 V; //velocity VECTOR3 W; //angular velocity Camera() : D(0.001f), //HACK: very near D helps prevent shadow problems F(100) {} Camera( float d, float f, float a ) : D(d), F(f) {}};//// 3 vertex indices//struct Triangle{ unsigned short V[3];//vertex indices //NOTE: No more than 64K vertices //default Triangle() { //more efficient to not initialize } //initialize Triangle( unsigned short v0, unsigned short v1, unsigned short v2 ) { V[0] = v0; V[1] = v1; V[2] = v2; }};//// Triangles, Vertices, Normals, Topology////NOTE: No more than 64K triangles or verticesstruct BoatGeometry{ Triangle *F;//faces unsigned short FCount;//number of faces, face normals VECTOR3 *V;//vertices VECTOR3 *N;//vertex normals float *dA; unsigned short VCount; BoatGeometry() : F(NULL), FCount(0), V(NULL), N(NULL), dA(NULL), VCount(0) {} ~BoatGeometry() { if( F ) delete F; if( V ) delete V; if( N ) delete N; if( dA ) delete dA; }};//// 2D Array of floats//struct FloatArray2D{ float* pData;//pointer to externally allocated data const long Width;//array width FloatArray2D( float* pd, const long w ) : pData(pd), Width(w) {} //indexing operator //read-only const float* operator [] ( const long i ) const { return pData + i*Width; } //write-to float* operator [] ( const long i ) { return pData + i*Width; }};//globalsCamera gCamera;RIGID_BODY gBoat;BoatGeometry gBoatGeometry;//water stuffconst long N = 128 + 1;//number of gridpoints along a sideconst float c = 2;//wave space, meters/secconst float h = 2*float(0.1*c*sqrt(2));//grid cell width (min for stability, assuming dt<=0.1)const float L = h*(N-1);//width of entire gridfloat z[N][N];//z[n] valuesfloat z1[N][N];//z[n-1] and z[n+1] valuesfloat d[N][N];//damping coefficientsFloatArray2D Z( (float*)z, N );//z[n] valuesFloatArray2D Z1( (float*)z1, N );//z[n-1] and z[n+1] valuesVECTOR3 V[N][N];//triangle verticesVECTOR3 VN[N][N];//vertex normalsGLint Tri[N-1][N][2];//triangle strip elements/////////////////////////////////////////////////////////////////////////////////////////void animateWater( FloatArray2D& z, FloatArray2D& z1, const float dt){ //Integrate the solution of the 2D wave equation. //OPTIMIZATION: This code has been written for //clarity, not speed! It would be better to use //pointer arithmetic instead of indexing. //precalculate coefficients const float A = (c*dt/h)*(c*dt/h); const float B = 2 - 4*A; long i, j; //edges are unchanged for( i=1 ; i<N-1 ; i++ ) { for( j=1 ; j<N-1 ; j++ ) { //integrate, replacing z[n-1] with z[n+1] in place z1[i][j] = A*( z[i-1][j] + z[i+1][j] + z[i][j-1] + z[i][j+1] ) + B*z[i][j] - z1[i][j]; //apply damping coefficients z1[i][j] *= d[i][j]; } } //swap pointers Swap( z.pData, z1.pData );}/////////////////////////////////////////////////////////////////////////////////////////void interpolateWater( const float x, const float y, float& zw, VECTOR3& nw ){ const float rh = 1 / h; //fractional index float u = (x + 0.5*L) * rh; float v = (y + 0.5*L) * rh; //lower-left vertex of the enclosing grid cell const long i = long( u ); const long j = long( v ); //interpolation coefficients const float a = u - i; const float b = v - j; const float ab = a*b; //if the position is outside of the grid, give a fake value if( i<0 || N<=i || j<0 || N<=j ) { zw = 0; nw.x = nw.y = 0; nw.z = 1; } else { //bilinearly interpolate zw and nw zw = (1-a-b+ab) * z[i][j] + (b-ab) * z[i][j+1] + (a-ab) * z[i+1][j] + ab * z[i+1][j+1]; nw = (1-a-b+ab) * VN[i][j] + (b-ab) * VN[i][j+1] + (a-ab) * VN[i+1][j] + ab * VN[i+1][j+1]; nw.normalize();//if the interpolation was spherical, this wouldn't be necessary /* //DEBUG: interpolated height should be between z_min and z_max float z_min = Min( Min(z[i][j],z[i][j+1]), Min(z[i+1][j],z[i+1][j+1]) ); float z_max = Max( Max(z[i][j],z[i][j+1]), Max(z[i+1][j],z[i+1][j+1]) ); if( !(z_min<=zw && zw<=z_max) ) { int kd=0; } //*/ }}/////////////////////////////////////////////////////////////////////////////////////////void animateBoat( const float dt ){ const float k = 1000 * 9.8;//water density * gravity VECTOR3 F_net(0,0,0);//net force VECTOR3 T_net(0,0,0);//net torque //gravity F_net.z -= 9.8 * gBoat.M; //buoyancy for( long i=0 ; i<gBoatGeometry.VCount ; i++ ) { //vector from CM to hull vertex, in world space const VECTOR3 r = gBoat.transformVectorToParent( gBoatGeometry.V[i] ); const VECTOR3 p = gBoat.position() + r;//hull vertex VECTOR3 nw;//water normal float zw;//water height //bilinearly interpolate the height //and normal of the water surface //at the (x,y) position interpolateWater( p.x, p.y, zw, nw ); //only calculate buoyancy //for submerged portins if( zw > p.z ) { //velocity of hull vertex with respect to water, assuming water is still const VECTOR3 v = gBoat.V + gBoat.W.cross(r); //drag force is proportional to relative velocity const VECTOR3 Fd = - 2500 * v;//HACK: magic number //buoyant force is proportional to weight of water displaced const VECTOR3 Fb = - k * (zw - p.z) * gBoatGeometry.dA[i] * gBoatGeometry.N[i].z * nw; //total force exerted on this hull vertex const VECTOR3 F = Fd + Fb; F_net += F;//add to force T_net += r.cross(F);//add to torque } } const VECTOR3 dv = (1 / gBoat.M) * F_net * dt; const VECTOR3 dw = gBoat.dwdt( T_net ) * dt; gBoat.translate( gBoat.V * dt ); gBoat.rotate( gBoat.W * dt ); gBoat.V += dv; gBoat.W += dw; gBoat.W *= 0.999f;//HACK: fake energy loss prevents instability}/////////////////////////////////////////////////////////////////////////////////////void idle(){ const float dt = 0.025f;//HACK: constant time step //animate the camera gCamera.translate( gCamera.V * dt ); //animate the water animateWater( Z, Z1, dt ); //animate the boat animateBoat( dt ); //redraw glutPostRedisplay();}/////////////////////////////////////////////////////////////////////////////////////void drawScene(){ long i, j; //draw the water //update the vertices and the normals for( i=0 ; i<N; i++ ) { for( j=0 ; j<N; j++ ) { //only the z needs to be updated V[i][j].z = z[i][j]; //only the x and y need to be updated VN[i][j].x = (z[i][j-1] - z[i][j+1]); VN[i][j].y = (z[i-1][j] - z[i+1][j]); //OPTIMIZATION: For clarity, the //normals are updated here, but it //might be faster to update them //inside the integration loop since //the neighboring points have already //been loaded into the cache. } } //draw triangles glEnableClientState( GL_VERTEX_ARRAY ); glEnableClientState( GL_NORMAL_ARRAY ); glVertexPointer( 3, GL_FLOAT, 0, V ); glNormalPointer( GL_FLOAT, 0, VN ); //allow the driver to normalize the normals. //If the video card has T&L, it might be //faster than doing it in the CPU. glEnable( GL_NORMALIZE ); glColor3f( 0, 0, 1 );//blue glPolygonMode( GL_FRONT, GL_LINE ); glPolygonMode( GL_BACK, GL_LINE ); for( i=0 ; i<N-1; i++ ) { glBegin( GL_TRIANGLE_STRIP ); { for( j=0 ; j<N; j++ ) { glArrayElement( Tri[i][j][0] ); glArrayElement( Tri[i][j][1] ); } } glEnd(); } glDisable( GL_NORMALIZE ); glDisableClientState( GL_NORMAL_ARRAY ); glDisableClientState( GL_VERTEX_ARRAY ); //draw the boat glColor3f( 1, 1, 1 );//white glPolygonMode( GL_FRONT, GL_FILL );//solid glPolygonMode( GL_BACK, GL_FILL ); glEnableClientState( GL_VERTEX_ARRAY ); glEnableClientState( GL_NORMAL_ARRAY ); glVertexPointer( 3, GL_FLOAT, 0, gBoatGeometry.V ); glNormalPointer( GL_FLOAT, 0, gBoatGeometry.N ); glPushMatrix(); { const VECTOR3& R0 = gBoat.R[0]; const VECTOR3& R1 = gBoat.R[1]; const VECTOR3& R2 = gBoat.R[2]; const VECTOR3& O = gBoat.O;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -