MATLAB有限元二维编程(三角单元)
- 首先,原理我是参考了这位博主,博主讲的很细致,但是有一些细节没有详细的提及,在此把自己在参照该博主的matlab代码后自己编程的过程纪录一下,如果有人看博主的文章就已经懂了,那么我这篇文章可能就对这些游客多余了。https://blog.csdn.net/lusongno1/article/details/81125167
- 博主是自己画的网格,COMSOL有导出网格的功能,所以我是借助了COMSOL里画了一个简单的模型然后导出网格再用MATLAB处理的,自我感觉网格的大小可以自定义,所以选择这个方法。
- 处理网格,因为COMSOL导出的网格有自己的格式,且其边界单元有一个缺点,不会区分给出的节点是那条边上的,本问题中由于边界都取0,所以没什么影响,但是如果狄式边界取值为非零值,则COMSOL不会告诉你哪条边为非零值。当然,把边界都设为0,作为有限元二维编程入门足够了。网格的数据格式如下
- 进入正题,求解的题目为:
- 此时我有必要说明一下,希望大家不会犯我一样的错误,以前我以为借助编程求解函数,那么很多步骤其实自己只要了解就好,不用从头到尾的把函数求解一遍,感觉那样跟自己手动求解没什么两样。这种想法有缺陷,编程实际上就是把你整个推导过程写成程序,如果你连每一步怎么来的,接下来该怎么走都不清楚,又何谈把它编成程序,所以先把方程整个推导一遍是基础。在我理解看来,编程也就是在算积分,单元总装,求解部分有用,其他都是这几步的铺垫。
- 推导过程:
%读入网格 clear all; filename='mm1.txt'; %%节点坐标 [node_num]=textread(filename,'%n',1);%节点个数 n_coor=zeros(node_num,2);%节点坐标 m=1; [n1,n2]=textread(filename,'%n%n','headerlines',1); n_coor(:,1)=n1(1:node_num); n_coor_x=n1(1:node_num); n_coor_y=n2(1:node_num); n_coor(:,2)=n2(1:node_num); %%边界点 bou_num=n1(node_num+1); boundary=zeros(bou_num,2); a1=node_num+2; a2=a1+bou_num-1; boundary(:,1)=n1(a1:a2); boundary(:,2)=n2(a1:a2); %%单元索引 ele_num=n1(node_num+2+bou_num);%单元个数 ele_index=zeros(ele_num,3);%单元索引 a3=node_num+5+bou_num; [ele_index(:,1),ele_index(:,2),ele_index(:,3)]=textread(filename,'%n%n%n','headerlines',a3);%%计算单个矩阵 K=sparse(node_num,node_num); F=sparse(node_num,1); for i=1:ele_numke=caculate1(i,ele_index,n_coor);f=caculate2(i,ele_index,n_coor); q1=ele_index(i,1)+1;q2=ele_index(i,2)+1;q3=ele_index(i,3)+1;q=zeros(1,3);q(1)=q1;q(2)=q2;q(3)=q3;K(q,q)=K(q,q)+ke;F(q,1)=F(q,1)+f; end %施加边界条件 b=zeros(node_num,1); u_b=0; K(1,1)=1; for i=1:bou_numw1=boundary(i,1)+1;w2=boundary(i,2)+1;if(b(w1)==0)F(w1,1)=u_b;K(w1,:)=0;K(:,w1)=0;K(w1,w1)=1.0;b(w1)=1;endif(b(w2)==0)F(w2,1)=u_b; K(w2,:)=0;K(:,w2)=0;K(w2,w2)=1;b(w2)=1;end end u=K\F; subplot(1,2,1); plot3(n_coor_x,n_coor_y,u); %scatter3(n_coor(:,1),n_coor(:,2),u); L =1; nsamp = 1001; xsamp = linspace(0,L,nsamp);%100等分区间中间有100个数 [X,Y] = meshgrid(xsamp,xsamp); uexact = sin(pi*X).*sin(pi*Y); uexact_re = reshape(uexact,nsamp,nsamp); subplot(1,2,2); mmm=mesh(xsamp,xsamp,uexact_re);%%%%%function [k] = caculate1(i,ele_index,n_coor) %UNTITLED11 此处显示有关此函数的摘要 % 此处显示详细说明 a1=ele_index(i,1)+1; a2=ele_index(i,2)+1; a3=ele_index(i,3)+1; x1=n_coor(a1,1); y1=n_coor(a1,2); x2=n_coor(a2,1); y2=n_coor(a2,2); x3=n_coor(a3,1); y3=n_coor(a3,2); A=0.5*((x2*y3-x3*y2)-(y3-y2)*x1+y1*(x3-x2)); A=abs(A); A=1/(2*A); J=(x3-x1)*(y2-y3)-(y3-y1)*(x2-x3); J1=A*[(y2-y3) (x3-x2);(y3-y1) (x1-x3)];a11=J.*([1 0]*J1*(J1'*[1;0])).*0.5; a12=J.*([0 1]*J1*(J1'*[1;0])).*0.5; a13=J.*([-1 -1]*J1*(J1'*[1;0])).*(0.5); a22=J.*([0 1]*J1*(J1'*[0;1])).*0.5; a23=J.*([-1 -1]*J1*(J1'*[0;1])).*(0.5); a33=J.*(([-1 -1]*J1)*(J1'*[-1;-1]))*(0.5); k=[a11 a12 a13; a12 a22 a23;a13 a23 a33]; endfunction [f] = caculate2(i,ele_index,n_coor) %UNTITLED12 此处显示有关此函数的摘要 % 此处显示详细说明 a1=ele_index(i,1)+1; a2=ele_index(i,2)+1; a3=ele_index(i,3)+1; x1=n_coor(a1,1); y1=n_coor(a1,2); x2=n_coor(a2,1); y2=n_coor(a2,2); x3=n_coor(a3,1); y3=n_coor(a3,2); J=(x3-x1)*(y2-y3)-(y3-y1)*(x2-x3); f1=@(lam1,lam2) fun(x1*lam1+x2*lam2+x3*(1-lam2-lam1),y1*lam1+y2*lam2+y3*(1-lam2-lam1)).*lam1.*J; f2=@(lam1,lam2) fun(x1*lam1+x2*lam2+x3*(1-lam2-lam1),y1*lam1+y2*lam2+y3*(1-lam2-lam1)).*lam2.*J; f3=@(lam1,lam2) fun(x1*lam1+x2*lam2+x3*(1-lam2-lam1),y1*lam1+y2*lam2+y3*(1-lam2-lam1)).*(1-lam1-lam2).*J; lam=@(lam1) 1-lam1; g1=integral2(f1,0,1,0,lam); g2=integral2(f2,0,1,0,lam); g3=integral2(f3,0,1,0,lam); f=zeros(3,1); f(1)=g1; f(2)=g2; f(3)=g3; endfunction bx = fun(x,y) bx = (2*pi^2)*sin(pi*x).*sin(pi*y); end
- 真解和所求解对比图
MATLAB有限元二维编程(三角单元)相关推荐
- MATLAB编程(4)——MATLAB绘制二维高斯函数的三维图
本篇博文记录使用MATLAB绘制二维高斯函数的三维图. 用到的MATLAB函数--mesh()(绘制三维线框图)和surf()(绘制三维表面图). MATLAB命令窗口输入>> doc 函 ...
- matlab中的delaunay,基于MATLAB 实现二维delaunay 三角剖分
基于MATLAB 实现二维delaunay 三角剖分 刘锋涛凡友华 (哈尔滨工业大学深圳研究生院深圳518055) [摘要]在已知凸多边形的顶点坐标的前提情况下,利用MATLAB 中的meshgrid ...
- matlab建成二维数组,matlab绘制二维数组
hist 累计图 rose 极座标累计图 stairs 阶梯图 stem 针状图 fill 实心图 feather 羽毛图 compass 罗盘图 quiver 向量场图 Matlab 如何画出一个二 ...
- matlab生成二维服从高斯分布的数据
matlab生成二维服从高斯分布的数据 2015-12-30 21:31 1263人阅读 评论(0) 收藏 举报 分类: matlab(8) 由于实验需要,需要生成两类模式的数据,同时这两类数据要 ...
- matlab绘制二维曲线图
matlab绘制二维曲线图 今天,我们来讲一个用matlab绘制二维曲线图 下面直接上代码,会对代码一些部分进行一些讲解 %% 定义函数 x = 0:0.01:2*pi; y1 = sin(x); y ...
- MATLAB绘制二维曲线-fplot函数
MATLAB绘制二维曲线-fplot函数 fplot函数的基本用法 双输入函数参数的用法 fplot函数的基本用法 fplot(f,lims,选项) f代表一个函数,通常使用函数句柄的形式,lims为 ...
- Matlab:二维傅里叶变换
Matlab:二维傅里叶变换 二维傅里叶变换 二维衍射模式 fft2 函数将二维数据变换为频率空间.例如,您可以变换二维光学掩膜以揭示其衍射模式. 二维傅里叶变换 以下公式定义 m×n 矩阵 X 的离 ...
- 利用Matlab 解决二维矩阵问题
写在前面 Matlab是一款非常强大的数学计算工具,学习并使用它进行处理一些数据运算,将会非常之高效. 今天有同学问我了一道关于利用Matlab 解决二维矩阵问题,利用空闲时间给他解答,希望能帮助到他 ...
- matlab 极坐标 二维,matlab笔记二维绘图(极坐标隐函数等)008.docx
matlab笔记二维绘图(极坐标隐函数等)008.docx 008二维绘图(极坐标.隐函数等)一.极坐标图形调用格式为POLART,R,'选项'其中,T为极角,R为极径,选项的使用和PLOT类似.例1 ...
最新文章
- “干掉” Date,Java8 LocalDate 真香!
- c++ static用法,全局变量,与别的语言不一样
- 微信小程序image组件开发程序以及相关图片问题参考资料汇总
- SQLite学习手册(命令行工具)
- itstime后面跟什么_被父母当成摇钱树是种什么体验?
- 如何在ASP.NET Core程序启动时运行异步任务(1)
- 是可改写的随机存储器_PPT下载:磁阻随机存取存储器
- 无法定位软件包 docker-ce_自媒体!做自媒体账号需要注意什么?定位很重要
- 拍拍贷第三届“魔镜杯”启动 :10万美金邀你“秀出你的算法!”
- 开三个iframe不断刷访问量
- JMeter 修改字体大小
- 网件 无线打印机服务器,NETGEAR Genie让普通打印机实现Air print功能
- 生物信息学习——tophat使用手册
- VS code,Live Server更改默认浏览器
- 编程程序 runtime error
- 第三章:EB配置DIO输出(s32k144)
- C++打开网页,发起QQ对话,调用外部exe程序
- 大数乘法(快速傅立叶变换)下
- Windows畸形文件夹
- 源码学习之LAMMPS的一个时间步是如何工作的
热门文章
- 支付宝小程序-实名认证流程讲解
- Socks代理是什么?PC端怎么使用Socks5代理?
- JS基础----函数应用 案例
- html标签的message,Message 消息提示
- PEEKABOO——团队展示
- c#窗体程序生成错误_创建一个没有窗口的程序 (C#) | 学步园
- 解密-大象跳转如何实现微信中点击链接直接跳到默认浏览器(不是在微信内置浏览器打开)
- Synctoy2.1通过计划任务备份文件到网络驱动器注销不生效问题
- 超文本标记语言是指Java_超文本标记语言(HTML)
- w10取消自带杀毒服务器,win10家庭版关闭杀毒服务如何设置_win10家庭版怎么关闭自带杀毒系统-win7之家...