[opengl]弹簧质点模型的opengl实现 opengl

实现了一种基于弹簧质点模型的布料仿真程序
论文见：http://wenku.baidu.com/link?url=8rBDraTsWTtTL8cFU0vjhXMyv4RF0npjeizz2CQQI4DvvTrxsN3bNOK91_1jRw7TuVadHuds5VnWzQ8CxP0QSOcY6sEUMKUib44crkbil0K

对于超弹问题，依然没能很好地解决，希望大神给予指导

#ifndef GLUT_DISABLE_ATEXIT_HACK #define GLUT_DISABLE_ATEXIT_HACK #endif #define GLEW_STATIC #include #include #include #include #include #include #include //for matrices #include //undefine if u want to use the default bending constraint of pbd #include using namespace std; //using namespace glm; GLfloat rtx = 0.0f, rty = 0.0f, rtz = 0.0; GLfloat sim_step = 0.1; glm::vec3 Up=glm::vec3(0,1,0), Right, viewDir; const int GRID_SIZE=10; //地板砖边长 GLdouble MV[16]; GLint viewport[4]; GLdouble PP[16]; bool isfix[1000]; int dist = -23; //椭球参数 int iStacks = 30; int iSlices = 30; float fRadius = 1; //float exp1 = 1e-3; const int numX = 20, numY=20; //一行有numx+1个点 const int total_points = (numX+1)*(numY+1); //总点数 //布料顶点位置速度 const int fck = 10000; glm::vec3 Pos[total_points]; glm::vec3 Veloc[total_points]; glm::vec3 force[total_points]; int size = 4; float hsize = size/2.0f; float currentTime,newTime,startTime; const float mass = 1.0/total_points; const glm::vec3 gvat = glm::vec3(0,-9.8,0); //重力加速度 const float Spring_K = 2.5; //弹性系数 const float len0 = 4.0/numX; //单边长度 const float tolera = 1.08; //弹性限度void initGL() { //初始化顶点位置 memset(Pos,0,sizeof(Pos)); memset(Veloc,0,sizeof(Veloc)); memset(force,0,sizeof(force)); //fill in positions int count1 = 0; int u = numX + 1; int v = numY + 1; for(int j=0; j<=numY; j++) { for(int i=0; i<=numX; i++) { Pos[count1++] = glm::vec3( ((float(i)/(u-1)) *2-1)* hsize, size+1, ((float(j)/(v-1) )* size)); printf("(%.1lf ,%.1lf)",((float(i)/(u-1)) *2-1)* hsize,((float(j)/(v-1) )* size)); }printf("\n"); //悬挂点为X[0] 和 X[numX] } memset(isfix,0,sizeof(isfix)); isfix[0] = isfix[numX] = 1; //处理时间 startTime = (float)glutGet(GLUT_ELAPSED_TIME); cout<1){//上 force[k] += SpringForce(k, k - 2*numX - 2, 2); } if(j1){//左 force[k] += SpringForce(k, k - 2, 2); } k++; } } //剪切弹簧 k = 0; for(i = 0; i <= numY; i ++){ for(j = 0; j <= numX; j ++){ if(i>0&&j>0){//左上 force[k] += SpringForce(k, k - numX - 2, 3); } if(i>0&&j0){//坐下 force[k] += SpringForce(k, k + numX, 3); } k ++; } } }void CalcPos(){//计算新的位置 newTime = (float) glutGet(GLUT_ELAPSED_TIME); float frameTime = newTime-currentTime; frameTime/=1000; currentTime = newTime; glm::vec3 acc = glm::vec3(0); //得到的frame就是间隔时间毫秒 for(size_t i = 0; i <= total_points; i ++){ if(isfix[i])continue; acc = force[i] / mass; //得到加速度向量 Veloc[i] *= 0.98; Veloc[i] = Veloc[i] + acc*frameTime; //得到新的速度值 Pos[i] = Pos[i] + Veloc[i] * frameTime; } }void gao(int a,int b,int is){//直线弹簧要考虑是否是isfix！ glm::vec3 tmp = Pos[b] - Pos[a]; float dis = glm::length(tmp); float inilen; if(is == 1)inilen = len0; if(is == 2) inilen = len0 * 2; if(is == 3) inilen = len0 * 1.414213; inilen *= tolera; if(dis<=inilen)return; glm::vec3 tmp1 = tmp*(inilen/dis); if(isfix[a]&&isfix[b]){ return; }else if(isfix[a]){ Pos[b] = Pos[a] + tmp1; }else if(isfix[b]){ Pos[a] = Pos[b] - tmp1; }else{ glm::vec3 mid = (Pos[a] + Pos[b]) * 0.5f; Pos[a] = mid - tmp1 * 0.5f; Pos[b] = mid + tmp1 * 0.5f; } }//限制每一根弹簧长度 void DynamicConstrain(){ int i,j,k = 0; for(i = 0; i <= numY; i++){ for(j = 0; j <= numX; j ++){ if(j!=numX){//右 gao(k,k+1,1); } if(i!=numY){//下 gao(k,k + numX + 1,1); }if(j0){//左下 gao(k, k + numX, 3); } k++; } } } void StepPhysics(){ ComputeForces(); CalcPos(); glutPostRedisplay(); //DynamicConstrain(); Sleep(5); }void OnReshape(int nw, int nh) { glViewport(0,0,nw, nh); glMatrixMode(GL_PROJECTION); glLoadIdentity(); gluPerspective(60, (GLfloat)nw / (GLfloat)nh, 1.f, 100.0f); glGetIntegerv(GL_VIEWPORT, viewport); glGetDoublev(GL_PROJECTION_MATRIX, PP); glMatrixMode(GL_MODELVIEW); } int main(int argc, char * argv[]) { glutInit(&argc, argv); glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGBA | GLUT_DEPTH); glutInitWindowSize(1024, 1024); glutCreateWindow("Mass-Spring Model"); initGL(); glutDisplayFunc(OnRender); //指定窗口形状变化时的回调函数 glutReshapeFunc(OnReshape); //指定程序空闲时调用函数 glutIdleFunc(StepPhysics); glEnable(GL_DEPTH_TEST); glutMainLoop(); return 0; }

实现效果