科学计算可视化报告.docx
- 文档编号:8908631
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:20
- 大小:168.91KB
科学计算可视化报告.docx
《科学计算可视化报告.docx》由会员分享,可在线阅读,更多相关《科学计算可视化报告.docx(20页珍藏版)》请在冰豆网上搜索。
科学计算可视化报告
《科学计算可视化》实验报告
等值线生成
tjuwar
天津大学计算机科学与技术学院
2011年4月28日
一、实验目的
编程实现等值线生成的网格序列法和单元剖分法
二、实验内容
数据可采用随机方法生成,代表网格中每个点的值。
根据网格序列法和单元剖分法对网格的点进行扫描,并由给出的值生成等值线。
实现方法:
1.网格序列法:
假设网格单元都是矩形,其等值线生成算法的主要步骤如下:
1)逐个计算每一个网格单元与等值线的交点;
2)连接该单元内等值线的交点,生成该单元内的等值线线段;
3)由一系列单元内的等值线线段构成该网格中的等值线;
网格单元与等值线的交点计算主要计算各单元边与等值线的交点,可采用顶点判定,边上插值的方法计算。
设等值线的值为Ft,若Fij(Ft,则记顶点为‘-’;若Fij>Ft,则记顶点为‘+’。
若单元的四个顶点全为‘+’或‘-’,则网格单元内无等值线;否则对两个顶点分别为’+’‘-‘的单元边插值计算等值线的交点,并在单元内连线,连线情况见图2-2。
图1单元内等值线连接情况
在图1(d)的情况下,实际上存在着两种连接方式的二义情况,不可能判断哪种连接情况是正确的。
这里的二义性问题通用中点的值加以判断来解决。
2.单元剖分法:
为了避免网格序列法的二义性情况,可采用单元剖分法,算法的基本思想是利用对角线将矩形单元分成四个三角形单元,求出中心点的函数值,等值线的抽取直接在三角单元中进行。
三角单元中至多只包含一条等值线,从而避免了二义性问题,但处理单元数目增加了四倍。
图2剖分法连接情况
三、实验结果
1.网格序列法
图3网格宽30,网格数10
图4网格宽20,网格数15
2.单元剖分法
图5网格宽30,网格数10
图6网格宽20,网格数15
四、实验分析和总结
网格序列法计算量小,有二义性,这里用中间点的值判断来解决,处理效果较粗糙。
而单元剖分法运算点多,不存在二义性问题,但计算量大。
在实现算法的过程中熟练了插值的方法。
五、源代码
注:
运行说明:
要切换序列法和单元剖分法只需在
voiddisplay(void)
{
//display_1();//网格序列法
display_2();//单元剖分法
}中切换。
/*
王安然3008216120
等值线生成
网格序列法和单元剖分法
2011年4月28日
*/
#include
#include
#include"math.h"
#include
#include
#include
#include
#include
#include
#include"math.h"
usingnamespacestd;
#definemax100//最大网格数
floatA[max][max];
intFlag[max][max];
constfloatB=5;
intN=10;//实际网格数
intgrid=30;//网格宽度
//初始化OpenGL
voidinit(void)
{
glClearColor(0.0f,0.0f,0.0f,0.0f);//设置背景颜色
glShadeModel(GL_FLAT);//设置明暗处理
}
//主要的绘制过程
//网格法
voiddisplay_1(void)
{
glClear(GL_COLOR_BUFFER_BIT);//清除颜色缓存
glBegin(GL_LINES);//开始画直线
glColor3f(1.0f,1.0f,1.0f);//设置颜色为白色
inti,j;
floatx1,y1,x2,y2,x3,y3,x4,y4;
floatX0,X1,Y0,Y1;
intstartx=100;
intstarty=100;
for(i=0;i for(j=0;j { X0=startx+i*grid; Y0=starty+j*grid; X1=startx+(i+1)*grid; Y1=starty+(j+1)*grid; if(j! =N-1) //画点上方的线 { glVertex2f(X0,Y0); glVertex2f(X0,Y1); } if(i! =N-1) //画点右方的线 { glVertex2f(X0,Y0); glVertex2f(X1,Y0);} //左下角点与其他三点符号不同 glColor3f(1.0f,0.0f,0.0f);//设置颜色为红色 if( ( (Flag[i][j]&&! Flag[i][j+1]&&! Flag[i+1][j]&&! Flag[i+1][j+1]) ||(! Flag[i][j]&&Flag[i][j+1]&&Flag[i+1][j]&&Flag[i+1][j+1]) ) &&(j+1) ) { x1=X0; y1=(Y0*(A[i][j+1]-B)+Y1*(B-A[i][j]))/(A[i][j+1]-A[i][j]); y2=Y0; x2=(X0*(A[i+1][j]-B)+X1*(B-A[i][j]))/(A[i+1][j]-A[i][j]); glVertex2f(x1,y1); glVertex2f(x2,y2); } //左上角点与其他三点符号不同 elseif( ( (! Flag[i][j]&&Flag[i][j+1]&&! Flag[i+1][j]&&! Flag[i+1][j+1]) ||(Flag[i][j]&&! Flag[i][j+1]&&Flag[i+1][j]&&Flag[i+1][j+1]) ) &&(j+1) ) { x1=X0; y1=(Y0*(A[i][j+1]-B)+Y1*(B-A[i][j]))/(A[i][j+1]-A[i][j]); y2=Y1; x2=(X0*(A[i+1][j+1]-B)+X1*(B-A[i][j+1]))/(A[i+1][j+1]-A[i][j+1]); glVertex2f(x1,y1); glVertex2f(x2,y2); } //右下角点与其他三点符号不同 elseif( ( (! Flag[i][j]&&! Flag[i][j+1]&&Flag[i+1][j]&&! Flag[i+1][j+1]) ||(Flag[i][j]&&Flag[i][j+1]&&! Flag[i+1][j]&&Flag[i+1][j+1]) ) &&(j+1) ) { x1=X1; y1=(Y0*(A[i+1][j+1]-B)+Y1*(B-A[i+1][j]))/(A[i+1][j+1]-A[i+1][j]); y2=Y0; x2=(X0*(A[i+1][j]-B)+X1*(B-A[i][j]))/(A[i+1][j]-A[i][j]); glVertex2f(x1,y1); glVertex2f(x2,y2); } //右上角点与其他三点符号不同 elseif( ( (! Flag[i][j]&&! Flag[i][j+1]&&! Flag[i+1][j]&&Flag[i+1][j+1]) ||(Flag[i][j]&&Flag[i][j+1]&&Flag[i+1][j]&&! Flag[i+1][j+1]) ) &&(j+1) ) { x1=X1; y1=(Y0*(A[i+1][j+1]-B)+Y1*(B-A[i+1][j]))/(A[i+1][j+1]-A[i+1][j]); y2=Y1; x2=(X0*(A[i+1][j+1]-B)+X1*(B-A[i][j+1]))/(A[i+1][j+1]-A[i][j+1]); glVertex2f(x1,y1); glVertex2f(x2,y2); } //左右两边顶点符号不同 elseif( ( (! Flag[i][j]&&! Flag[i][j+1]&&Flag[i+1][j]&&Flag[i+1][j+1]) ||(Flag[i][j]&&Flag[i][j+1]&&! Flag[i+1][j]&&! Flag[i+1][j+1]) ) &&(j+1) ) { x1=(X0*(A[i+1][j+1]-B)+X1*(B-A[i][j+1]))/(A[i+1][j+1]-A[i][j+1]); y1=Y1; y2=Y0; x2=(X0*(A[i+1][j]-B)+X1*(B-A[i][j]))/(A[i+1][j]-A[i][j]); glVertex2f(x1,y1); glVertex2f(x2,y2); } //上下两边顶点符号不同 elseif( ( (! Flag[i][j]&&Flag[i][j+1]&&! Flag[i+1][j]&&Flag[i+1][j+1]) ||(Flag[i][j]&&! Flag[i][j+1]&&Flag[i+1][j]&&! Flag[i+1][j+1]) ) &&(j+1) ) { y1=(Y0*(A[i][j+1]-B)+Y1*(B-A[i][j]))/(A[i][j+1]-A[i][j]); x1=X0; x2=X1; y2=(Y0*(A[i+1][j+1]-B)+Y1*(B-A[i+1][j]))/(A[i+1][j+1]-A[i+1][j]); glVertex2f(x1,y1); glVertex2f(x2,y2); } //四个顶点符号交叉,即左上与右下相同,右上与左下相同 elseif( ( (! Flag[i][j]&&Flag[i][j+1]&&Flag[i+1][j]&&! Flag[i+1][j+1]) ||(Flag[i][j]&&! Flag[i][j+1]&&! Flag[i+1][j]&&Flag[i+1][j+1]) ) &&(j+1) ) { y1=(Y0*(A[i][j+1]-B)+Y1*(B-A[i][j]))/(A[i][j+1]-A[i][j]); x1=X0; x2=X1; y2=(Y0*(A[i+1][j+1]-B)+Y1*(B-A[i+1][j]))/(A[i+1][j+1]-A[i+1][j]); x3=(X0*(A[i+1][j+1]-B)+X1*(B-A[i][j+1]))/(A[i+1][j+1]-A[i][j+1]); y3=Y1; y4=Y0; x4=(X0*(A[i+1][j]-B)+X1*(B-A[i][j]))/(A[i+1][j]-A[i][j]); //判断中间的点符号,并做相应处理 intflag=0; if((A[i][j]+A[i+1][j]+A[i][j+1]+A[i+1][j+1])/4>=B) flag=1; elseflag=0; if(flag==Flag[i][j]) { glVertex2f(x1,y1); glVertex2f(x3,y3); glVertex2f(x2,y2); glVertex2f(x4,y4); } else { glVertex2f(x3,y3); glVertex2f(x2,y2); glVertex2f(x4,y4); glVertex2f(x1,y1); } } glColor3f(1.0f,1.0f,1.0f);//设置颜色为白色 } glEnd(); glFlush();//开始绘制 } //单元剖分法 //点的结构体 structPoint{ floatx,y; Point(float_x,float_y): x(_x),y(_y){} Point(){} voidprint(constcharname[]){ printf("%s%f%f\n",name,x,y); } }; //结点的结构体,包含点坐标,符号等信息 structNode{ PointP; floatF; intFlag; Node(Point_P,float_F,int_Flag): P(_P.x,_P.y),F(_F),Flag(_Flag){} Node(){} }; //插值函数 Pointinterperlate(Pointp0,floatf0,Pointp1,floatf1){ Pointp; p.x=(p0.x*(f1-B)+p1.x*(B-f0))/(f1-f0); p.y=(p0.y*(f1-B)+p1.y*(B-f0))/(f1-f0); returnp; } //在三角形内判断是否插值和插值位置 voidcalcTriangle(NodeN0,NodeN1,NodeN2) { Pointp1,p2; //如果三角形三个点符号相同那么不画线 if(N0.Flag==0&&N1.Flag==0&&N2.Flag==0)return; if(N0.Flag==1&&N1.Flag==1&&N2.Flag==1)return; if(N0.Flag==0&&N1.Flag==1&&N2.Flag==1) { p1=interperlate(N0.P,N0.F,N1.P,N1.F); p2=interperlate(N0.P,N0.F,N2.P,N2.F); } //如果其中两个点符号相同,那么标记两点,并画线 elseif(N0.Flag==0&&N1.Flag==0&&N2.Flag==1) { N0.P.print("N0.P");printf("N0.F%f",N0.F); N1.P.print("N1.P");printf("N1.F%f",N1.F); N2.P.print("N2.P");printf("N2.F%f",N2.F); p1=interperlate(N0.P,N0.F,N2.P,N2.F); p2=interperlate(N1.P,N1.F,N2.P,N2.F); } elseif(N0.Flag==0&&N1.Flag==1&&N2.Flag==0) { p1=interperlate(N0.P,N0.F,N1.P,N1.F); p2=interperlate(N2.P,N2.F,N1.P,N1.F); } elseif(N0.Flag==1&&N1.Flag==0&&N2.Flag==0) { p1=interperlate(N1.P,N1.F,N0.P,N0.F); p2=interperlate(N2.P,N2.F,N0.P,N0.F); } elseif(N0.Flag==1&&N1.Flag==0&&N2.Flag==1) { p1=interperlate(N1.P,N1.F,N0.P,N0.F); p2=interperlate(N1.P,N1.F,N2.P,N2.F); } elseif(N0.Flag==1&&N1.Flag==1&&N2.Flag==0) { p1=interperlate(N2.P,N2.F,N0.P,N0.F); p2=interperlate(N2.P,N2.F,N1.P,N1.F); } else { printf("error"); } glVertex2f(p1.x,p1.y);glVertex2f(p2.x,p2.y); printf("p1%f%fp2%f%f\n",p1.x,p1.y,p2.x,p2.y); } voiddisplay_2(void) { glClear(GL_COLOR_BUFFER_BIT);//清除颜色缓存 glBegin(GL_LINES);//开始画直线 glColor3f(1.0f,1.0f,1.0f);//设置颜色为白色 inti,j; floatx1,y1,x2,y2,x3,y3,x4,y4; floatX0,X1,Y0,Y1; PointP00,P01,P10,P11,PMid; floatMidX,MidY; floatF00,F01,F10,F11,FMid; NodeN00,N01,N10,N11,NMid; intFlagMid; intstartx=100; intstarty=100; for(i=0;i for(j=0;j { X0=startx+i*grid; Y0=starty+j*grid; X1=startx+(i+1)*grid; Y1=starty+(j+1)*grid; printf("x0y0%f%f%f%f\n",X0,Y0,X1,Y1); MidX=(X0+X1)/2; MidY=(Y0+Y1)/2; F00=A[i][j];F01=A[i][j+1];F10=A[i+1][j];F11=A[i+1][j+1]; FMid=(F00+F01+F10+F11)/4; if(FMid>=B) { FlagMid=1; }else { FlagMid=0; } glVertex2f(X0,Y0);glVertex2f(X0,Y1); glVertex2f(X0,Y1);glVertex2f(X1,Y1); glVertex2f(X1,Y1);glVertex2f(X1,Y0); glVertex2f(X1,Y0);glVertex2f(X0,Y0); //网格对角线 glVertex2f(X0,Y0);glVertex2f(X1,Y1); glVertex2f(X0,Y1);glVertex2f(X1,Y0); //对每个点和结点赋值 glColor3f(1.0f,0.0f,0.0f);//设置颜色为红色 P00=Point(X0,Y0);P01=Point(X0,Y1);P10=Point(X1,Y0);P11=Point(X1,Y1);PMid=Point(MidX,MidY); N00=Node(P00,F00,Flag[i][j]); N01=Node(P01,F01,Flag[i][j+1]); N10=Node(P10,F10,Flag[i+1][j]); N11=Node(P11,F11,Flag[i+1][j+1]); NMid=Node(PMid,FMid,FlagMid); //对格内4个三角形分别进行计算并画线 calcTriangle(N00,N01,NMid); calcTriangle(N01,N11,NMid); calcTriangle(N11,N10,NMid); calcTriangle(N10,N00,NMid); glColor3f(1.0f,1.0f,1.0f);//设置颜色为白色 } glEnd(); glFlush();//开始绘制 } voiddisplay(void) { //display_1();//网格法 display_2();//单元剖分法 } //在窗口改变大小时调用 voidreshape(intwidth,intheight) { glViewport(0,0,width,height);//设置视口 glMatrixMode(GL_PROJECTION);//设置当前为投影变换模式 glLoadIdentity();//用单位矩阵替换当前变换矩阵 gluOrtho2D(0.0,width,0.0,height);//设置正交投影视图体 } intmain(intargc,char
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 科学 计算 可视化 报告