`
暴风雪
  • 浏览: 388077 次
  • 性别: Icon_minigender_2
  • 来自: 杭州
社区版块
存档分类
最新评论

[opengl]弹簧质点法模拟柔性布料以及椭球碰撞的opengl实现

 
阅读更多

代码大体与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)) *2-1)* hsize,((float(j)/(v-1) )* size));
		}printf("\n");
		//悬挂点为X[0] 和 X[numX]
	}
	memset(isfix,0,sizeof(isfix));
	isfix[0] = isfix[numX] = 1;
	fuck = 0;


}
void DrawEllipsoid(){//画绿色椭球
    //设置椭球
    ellipsoid = glm::translate(glm::mat4(1),glm::vec3(0,fuck,0));
    fuck += 0.01;
	if(fuck>=4.0)fuck = -1;
	ellipsoid = glm::rotate(ellipsoid, 45.0f ,glm::vec3(1,0,0));
	ellipsoid = glm::scale(ellipsoid, glm::vec3(fRadius,fRadius,fRadius/2));
	invreseEllipsoid = glm::inverse(ellipsoid);

	fuck +=0.005;

    //画出椭球
    glColor3f(0,1,0);
	glPushMatrix();
		glMultMatrixf(glm::value_ptr(ellipsoid));
        //半径,经线条数,维线条数
        glutWireSphere(fRadius, iSlices, iStacks);
	glPopMatrix();
}
void DrawGrid()  //画地板
{
	glBegin(GL_LINES);
	glColor3f(0.5f, 0.5f, 0.5f);
	for(int i=-GRID_SIZE;i<=GRID_SIZE;i++)
	{
		glVertex3f((float)i,0,(float)-GRID_SIZE);
		glVertex3f((float)i,0,(float)GRID_SIZE);

		glVertex3f((float)-GRID_SIZE,0,(float)i);
		glVertex3f((float)GRID_SIZE,0,(float)i);
	}
	glEnd();
}

void drawTextile(){  //画布料
	//draw polygons
	int k = 0;
    glBegin(GL_LINES);
    glColor3f(1,1,1);
	for(int i = 0;i <= numX;i ++){
        for(int j = 0;j <= numY;j ++){
//                cout<<i<<" "<<j<<";";
            if(j != numX){
                glm::vec3 p1 = Pos[k];
                glm::vec3 p2 = Pos[k+1];
                glVertex3f(p1.x,p1.y,p1.z);
                glVertex3f(p2.x,p2.y,p2.z);
//                printf("add %d %d\n",k,k+1);
            }
            if(i != numY){
                glm::vec3 p1 = Pos[k];
                glm::vec3 p2 = Pos[k+numX+1];
                glVertex3f(p1.x,p1.y,p1.z);
                glVertex3f(p2.x,p2.y,p2.z);
            }
            k++;
        }
//        cout<<endl;
	}
    glEnd();

    glPointSize(3);
    glBegin(GL_POINTS);
    glColor3f(1,0,0);
	for(size_t i = 0;i < total_points;i ++){
        glm::vec3 p = Pos[i];
		glVertex3f(p.x,p.y,p.z);
	}
    glEnd();
}

void OnRender(void)
{
    glClear(GL_COLOR_BUFFER_BIT| GL_DEPTH_BUFFER_BIT);
    glLoadIdentity();
    //设置视角
	glTranslatef(0,0,dist);
//	glRotatef(15,1,0,0);
	glRotatef(rX,1,0,0);
	glRotatef(rY,0,1,0);
    glGetDoublev(GL_MODELVIEW_MATRIX, MV);
	viewDir.x = (float)-MV[2];
	viewDir.y = (float)-MV[6];
	viewDir.z = (float)-MV[10];
	Right = glm::cross(viewDir, Up);
	//做出地面
    DrawGrid();
    drawTextile();
    DrawEllipsoid();
//    drawblock();
    glutSwapBuffers();
}
//求出b点对a点的作用力
glm::vec3 SpringForce(int a, int b, int is){
    float inilen;
    if(is == 1)inilen = len0;
    if(is == 2) inilen = len0 * 2;
    if(is == 3) inilen = len0 * 1.414213;
    glm::vec3 res = glm::vec3(0);
    glm::vec3 tmp = Pos[b] - Pos[a];
    float dis = glm::length(tmp);
    res = tmp/dis;
    res *= (dis - inilen);
    res *= Spring_K;
    return res;
}
void ComputeForces(){
    //重力
	for(size_t i=0;i<total_points;i++) {
		force[i] = glm::vec3(0);
		if(!isfix[i])
            force[i] += mass * gvat;
	}
	int i,j,k = 0;
    //结构弹簧 上下左右
    for(i = 0;i <= numY;i++){
        for(j = 0;j <= numX;j ++){
            if(i!=0){//上
                force[k] += SpringForce(k, k - numX - 1, 1);
            }
            if(j!=numX){//右
                force[k] += SpringForce(k, k+1, 1);
            }
            if(i!=numY){//下
                force[k] += SpringForce(k, k + numX + 1, 1);
            }
            if(j!=0){//左
                force[k] += SpringForce(k, k - 1, 1);
            }
            k++;
        }
    }
    //柔性弹簧
    k = 0;
    for(i = 0;i <= numY;i++){
        for(j = 0;j <= numX;j ++){
            if(i>1){//上
                force[k] += SpringForce(k, k - 2*numX - 2, 2);
            }
            if(j<numX-1){//右
                force[k] += SpringForce(k, k+2, 2);
            }
            if(i<numY-1){//下
                force[k] += SpringForce(k, k + 2*numX + 2, 2);
            }
            if(j>1){//左
                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&&j<numX){  //右上
                force[k] += SpringForce(k, k - numX, 3);
            }
            if(i<numY&&j<numX){  //右下
                force[k] += SpringForce(k, k + numX + 2, 3);
            }
            if(i<numY&&j>0){   //坐下
                force[k] += SpringForce(k, k + numX, 3);
            }
            k ++;
        }
    }
}

void EllipsoidCollision() { //检测与处理布料和椭球之间的碰撞
	for(size_t i=0;i<total_points;i++) {
        glm::vec3 tPos(Pos[i]);
		glm::vec4 X_0 = (invreseEllipsoid*glm::vec4(tPos,1));
		glm::vec3 delta0 = glm::vec3(X_0.x, X_0.y, X_0.z);
		float distance = glm::length(delta0);
		if (distance < 1.0f) {
			delta0 = (radius - distance) * delta0 / distance;
			// Transform the delta back to original space
			glm::vec3 delta;
			glm::vec3 transformInv;
			transformInv = glm::vec3(ellipsoid[0].x, ellipsoid[1].x, ellipsoid[2].x);
			transformInv /= glm::dot(transformInv, transformInv);
			delta.x = glm::dot(delta0, transformInv);
			transformInv = glm::vec3(ellipsoid[0].y, ellipsoid[1].y, ellipsoid[2].y);
			transformInv /= glm::dot(transformInv, transformInv);
			delta.y = glm::dot(delta0, transformInv);
			transformInv = glm::vec3(ellipsoid[0].z, ellipsoid[1].z, ellipsoid[2].z);
			transformInv /= glm::dot(transformInv, transformInv);
			delta.z = glm::dot(delta0, transformInv);
			tPos +=  delta ;
            Veloc[i] += (tPos - Pos[i])/frameTime;
            Veloc[i] *= globalDamp;
            Pos[i] = tPos;
		}
	}
}

void CalcPos(){  //计算新的位置
//	frameTime/=1000;
//	printf("%.3lf\n",frameTime);
	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] = Veloc[i] + acc*frameTime;  //得到新的速度值
        Veloc[i] *= globalDamp;
        Pos[i] = Pos[i] + Veloc[i] * frameTime;
    }
}

void StepPhysics(){
    ComputeForces();
    CalcPos();
    EllipsoidCollision();
    glutPostRedisplay();
    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);
}

void OnMouseDown(int button, int s, int x, int y)
{
	if (s == GLUT_DOWN)
	{
		oldX = x;
		oldY = y;
		int window_y = (height - y);
		float norm_y = float(window_y)/float(height/2.0);
		int window_x = x ;
		float norm_x = float(window_x)/float(width/2.0);

		float winZ=0;
		glReadPixels( x, height-y, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ );
		if(winZ==1)
			winZ=0;
		double objX=0, objY=0, objZ=0;
		gluUnProject(window_x,window_y, winZ,  MV,  P, viewport, &objX, &objY, &objZ);
		glm::vec3 pt(objX,objY, objZ);
		size_t i=0;
		for(i=0;i<total_points;i++) {
			if( glm::distance(Pos[i],pt)<0.1) {
				selected_index = i;
				printf("Intersected at %d\n",i);
				break;
			}
		}
	}

	if(button == GLUT_MIDDLE_BUTTON)
		state = 0;
	else
		state = 1;

	if(s==GLUT_UP) {
		selected_index= -1;
		glutSetCursor(GLUT_CURSOR_INHERIT);
	}
}

void OnMouseMove(int x, int y)
{
	if(selected_index == -1) {
		if (state == 0)
			dist *= (1 + (y - oldY)/60.0f);
		else
		{
			rY += (x - oldX)/5.0f;
			rX += (y - oldY)/5.0f;
		}
	} else {
		float delta = 1500/abs(dist);
		float valX = (x - oldX)/delta;
		float valY = (oldY - y)/delta;
		if(abs(valX)>abs(valY))
			glutSetCursor(GLUT_CURSOR_LEFT_RIGHT);
		else
			glutSetCursor(GLUT_CURSOR_UP_DOWN);

		Veloc[selected_index] = glm::vec3(0);
		Pos[selected_index].x += Right[0]*valX ;
		float newValue = Pos[selected_index].y+Up[1]*valY;
		if(newValue>0)
			Pos[selected_index].y = newValue;
		Pos[selected_index].z += Right[2]*valX + Up[2]*valY;
	}
	oldX = x;
	oldY = y;

	glutPostRedisplay();
}
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);
	glutMouseFunc(OnMouseDown);
	glutMotionFunc(OnMouseMove);
    glutMainLoop();
    return 0;
}

 

2
0
分享到:
评论
1 楼 wsyzyrxp 2016-01-13  
非常感谢 兄弟 帮了我大忙

相关推荐

    椭球碰撞检测 智能检测

    椭球碰撞检测是一种在计算机图形学中常见的物理模拟技术,特别是在游戏开发和三维模拟软件中。这个技术主要用于确定两个或多个椭球是否在空间中相互接触或者交叉,这对于实现复杂物体间的交互至关重要。椭球碰撞检测...

    椭球拟合_基于最小二乘的椭球拟合_椭球拟合

    椭球拟合是一种在三维空间中寻找最佳椭球形状来包容一组数据点的数学方法。在许多领域,如地质学、计算机视觉、图像处理和数据分析中,椭球拟合都有广泛应用。基于最小二乘法的椭球拟合是其中一种经典且有效的技术。...

    椭球法求解凸优化

    用于求解凸优化问题,是一种迭代收敛算法,将问题转化为凸问题后都可以用椭球法

    椭球拟合,椭球拟合 地磁校准,matlab

    这个文件可能包括了数据预处理、椭球参数初始化、最小二乘优化过程以及结果可视化等部分。通过学习和运行这个代码,开发者能够更好地理解椭球拟合的细节,并且能够将这个技术应用到自己的项目中。 在实际应用中,...

    基于椭球膨胀法实现独立坐标系统的建立.pdf

    在实际的工程实践中,文章提到三种常用的椭球膨胀法实现不同投影面的坐标转换: 1. 椭球长半径的改正量为投影面大地高。 2. 投影面的大地高等于卯酉曲率半径的变化量。 3. 以基点上参考椭球的平均曲率半径的变化值求...

    双椭球_双椭球热源模型_

    在有限元分析中,双椭球热源模型是一种高级的焊接热源模型,它被广泛应用于ABAQUS软件中,以模拟复杂的焊接过程。这种模型考虑了焊枪形状、熔池形成以及热量分布的非线性特性,为理解和预测焊接过程中材料的温度场、...

    激光熔覆数值模拟 COMSOL仿真 双椭球热源 采用双椭球热源模型,考虑材料热物性参数、相变、马兰戈尼效应、布辛涅斯克近似等,动

    激光熔覆数值模拟 COMSOL仿真 双椭球热源 采用双椭球热源模型,考虑材料热物性参数、相变、马兰戈尼效应、布辛涅斯克近似等,动网格模拟熔覆层,计算瞬态温度场和流场。

    椭球拟合代码和数据,使用matlab写的

    在IT领域,尤其是在数据分析、信号处理以及几何建模中,椭球拟合是一项重要的技术。本文将详细讨论基于MATLAB实现的椭球拟合过程,同时结合“Magnetometer-calibration”这一主题,探讨如何利用椭球拟合来校准磁力计...

    HFSS中画椭球

    在实际应用中,比如天线设计,椭球可能被用来模拟地球或者其他不规则形状的物体。在光学领域,椭球面镜片也是一种常见的元素。因此,熟练掌握在HFSS中绘制椭球的方法对于解决这些问题至关重要。 在提供的文件...

    椭球表面积计算器

    椭球是一种几何形状,类似于地球或其他行星的形状,它在三个轴上都有不同的半径。在数学和物理学中,椭球表面积的计算是相当重要的,特别是在地理学、工程学和计算机图形学等领域。本计算器专门针对椭球的表面积进行...

    DoubleEllipsoid_ansysworkbench_workbench椭球_

    在ANSYS Workbench环境中,双椭球热源模型是一种用于模拟复杂热问题的高级方法,尤其在处理非均匀热量分布时非常有用。本模型适用于分析设备内部或周围存在多个独立热源的情况,如电子元器件、机械设备或者热流体...

    15_TIG_20180507_焊接_fat5rq_双椭球_双椭球热源_abaqusDFLUX_

    在本主题中,我们主要关注的是使用ABAQUS软件进行焊接过程模拟,特别是涉及双椭球热源模型的实现。ABAQUS是一款强大的有限元分析软件,广泛应用于各种工程问题,包括热传导、结构力学以及复杂的耦合问题,如焊接过程...

    运用Matlab讨论椭球面性质.pdf

    此外,在Matlab中还可以实现椭球面的切平面和法线的设计。通过在特定点上的导数计算和梯度向量,可以得到椭球面在任意一点处的切平面方程和法线向量。这在进行几何分析或优化问题求解时非常有用。 利用Matlab的数据...

    ArcEngine椭球面积计算代码

    本篇将详细探讨"ArcEngine椭球面积计算代码"这一主题,以及如何基于ArcGIS进行多边形椭球面积的计算。 首先,理解“椭球面积”是关键。地球不是一个完美的球体,而是一个椭球体,因此在处理地理数据时,采用椭球...

    基于椭球平移法的高斯平面坐标转换(new)

    这些资源可以帮助我们理解并实现椭球平移法在高斯平面坐标转换中的应用。对于学习GIS技术或者进行地理信息处理的人来说,掌握这种转换方法是必不可少的技能,因为它能够确保在不同坐标系统之间准确无误地进行数据...

Global site tag (gtag.js) - Google Analytics