代码大体与http://bbezxcy.iteye.com/blog/2204322相同,修改了部分不稳定的bug,增加了椭球碰撞以及旋转观察的实现
#ifndef GLUT_DISABLE_ATEXIT_HACK
#define GLUT_DISABLE_ATEXIT_HACK
#endif
#define GLEW_STATIC
#include <GL/glew.h>
#include <GL/wglew.h>
#include <GL/freeglut.h>
#include <vector>
#include<cstring>
#include <glm/glm.hpp>
#include <glm/gtc/matrix_transform.hpp> //for matrices
#include <glm/gtc/type_ptr.hpp>
//undefine if u want to use the default bending constraint of pbd
#include<iostream>
using namespace std;
//using namespace glm;
GLfloat rtx = 0.0f, rty = 0.0f, rtz = 0.0;
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 radius = 1; //object space radius of ellipsoid
//float exp1 = 1e-3;
const int numX = 20, numY=20; //一行有numx+1个点
const int total_points = (numX+1)*(numY+1); //总点数
//布料顶点位置 速度
glm::vec3 Pos[total_points];
glm::vec3 Veloc[total_points];
glm::vec3 force[total_points];
glm::mat4 ellipsoid,invreseEllipsoid;
int size = 4;
float hsize = size/2.0f;
const float frameTime = 1.0f/60.0f;
const float mass = 1.0/total_points;
const float globalDamp = 0.98; //速度衰减参数
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; //弹性限度
float fuck;//球体运动
//视角问题
int oldX=0, oldY=0;
const int width = 1024, height = 1024;
GLdouble P[16];
int selected_index = -1;
int state =1 ;
float rX=15, rY=0;
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))