有限元基础编程-何晓明老师课件-一维程序实现matlab

news/2024/11/20 4:45:23/

文章目录

  • 前言
  • 一、主程序
  • 二、一维有限元求解程序-框架
  • 三、组装刚度矩阵
    • assemble_matrix_from_1D_integral.m
    • 2.1 算法
    • 2.2 get_standard_gauss_1D.m
    • 2.3 get_Gauss_local_1D.m

前言

只是为方便学习，不做其他用途，课程理论学习来自b站视频有限元基础编程-何晓明

一、主程序

在这里插入图片描述

function result = main_poisson_solver_1D(left,right,h_partition,basis_type,Gauss_point_number)
% 算例为何老师课件第一章85页的例1
% 何晓明老师的程序改编而来
% 一维泊松方程的有限元求解
% 我们程序中用“FE”代替“finite element”
% 试探(trial)FE函数和测试(test)FE函数需要相同
% 到目前为止，该代码只能处理一维线性有限元和一维拉格朗日二次有限元。
% 对于其他类型的FE，我们需要在代码中添加相应的信息。
% 问题域是[left,right]
% h_partition是均匀网格划分的步长
% basis_type:FE基函数的类型
% basis_type=101:1D 线性有限元% N_basis: N代表FE基函数个数，而不是网格划分个数
% N_partition: N表示网格划分个数，而不是FE基函数个数
% N：子区间的个数
% M_partition,T_partition, M_basis,T_basis:参考程序“generate_M_T_1D.m”中的注释 
% function_a: 泊松方程左边的系数函数
% funciton_f: 泊松方程的右边函数 
% function_g: Dirichelet边界函数在何老师第一章的课件 章节1-1
% function_q_tilde: the name of the Neumann boundary function q(x,y) when p(x,y)=0 in my notes "Notes for tool box of standard triangular FE" section 1-1.
% function_q: the name of the Neumann boundary function q(x,y) when p(x,y) is nonzero in my notes "Notes for tool box of standard triangular FE" section 1-1.
% function_p: the name of the Robin coefficient function p(x,y) in my notes "Notes for tool box of standard triangular FE" section 1-1.% h_basis 有限元节点的步长.N_partition = (right-left)/h_partition;if basis_type==101N_basis=N_partition;
end% 划分网格信息和有限元基函数
[M_partition,T_partition]=generate_M_T_1D(left,right,h_partition,101);if basis_type==101  %1D 线性有限元M_basis=M_partition;T_basis=T_partition;
end % 高斯积分在参考区间[-1,1]上的高斯点和权重
[Gauss_coefficient_reference_1D,Gauss_point_reference_1D] = generate_Gauss_reference_1D(Gauss_point_number);% 组装刚度矩阵
number_of_elements=N_partition;
matrix_size=[N_basis+1 N_basis+1];
vector_size=N_basis+1;
if basis_type==101number_of_trial_local_basis=2;number_of_test_local_basis=2;
end%组装刚度矩阵
A=assemble_matrix_from_1D_integral('function_a',M_partition,T_partition,T_basis,T_basis,number_of_trial_local_basis,number_of_test_local_basis,number_of_elements,matrix_size,Gauss_coefficient_reference_1D,Gauss_point_reference_1D,basis_type,1,basis_type,1);%组装载荷矢量
b = assemble_vector_from_1D_integral('function_f',M_partition,T_partition,T_basis,number_of_test_local_basis,number_of_elements,vector_size,Gauss_coefficient_reference_1D,Gauss_point_reference_1D,basis_type,0);%得到边界节点的信息矩阵
boundary_nodes = generate_boundary_nodes_1D(N_basis);%处理狄利克雷边界条件
[A,b]=treat_Dirichlet_boundary_1D('function_g',A,b,boundary_nodes,M_basis);%计算数值解
result=A\b;%计算所有节点上的最大误差
if basis_type==101h_basis=h_partition;
end
maxerror=get_maximum_error_1D(result,N_basis,left,h_basis);
maximum_error_at_all_nodes_of_FE = maxerrorend

二、一维有限元求解程序-框架

function result = main_poisson_solver_1D(left,right,h_partition,basis_type,Gauss_point_number)
% 算例为何老师课件第一章85页的例1
% 何晓明老师的程序改编而来
% 一维泊松方程的有限元求解
% 我们程序中用“FE”代替“finite element”
% 试探(trial)FE函数和测试(test)FE函数需要相同
% 到目前为止，该代码只能处理一维线性有限元和一维拉格朗日二次有限元。
% 对于其他类型的FE，我们需要在代码中添加相应的信息。
% 问题域是[left,right]
% h_partition是均匀网格划分的步长
% basis_type:FE基函数的类型
% basis_type=101:1D 线性有限元% N_basis: N代表FE基函数个数，而不是网格划分个数
% N_partition: N表示网格划分个数，而不是FE基函数个数
% N：子区间的个数
% M_partition,T_partition, M_basis,T_basis:参考程序“generate_M_T_1D.m”中的注释 
% function_a: 泊松方程左边的系数函数
% funciton_f: 泊松方程的右边函数 
% function_g: Dirichelet边界函数在何老师第一章的课件 章节1-1
% function_q_tilde: the name of the Neumann boundary function q(x,y) when p(x,y)=0 in my notes "Notes for tool box of standard triangular FE" section 1-1.
% function_q: the name of the Neumann boundary function q(x,y) when p(x,y) is nonzero in my notes "Notes for tool box of standard triangular FE" section 1-1.
% function_p: the name of the Robin coefficient function p(x,y) in my notes "Notes for tool box of standard triangular FE" section 1-1.% h_basis 有限元节点的步长.N_partition = (right-left)/h_partition;if basis_type==101N_basis=N_partition;
end% 划分网格信息和有限元基函数
[M_partition,T_partition]=generate_M_T_1D(left,right,h_partition,101);if basis_type==101  %1D 线性有限元M_basis=M_partition;T_basis=T_partition;
end % 高斯积分在参考区间[-1,1]上的高斯点和权重
[Gauss_coefficient_reference_1D,Gauss_point_reference_1D] = generate_Gauss_reference_1D(Gauss_point_number);% 组装刚度矩阵
number_of_elements=N_partition;
matrix_size=[N_basis+1 N_basis+1];
vector_size=N_basis+1;
if basis_type==101number_of_trial_local_basis=2;number_of_test_local_basis=2;
end%组装刚度矩阵
A=assemble_matrix_from_1D_integral('function_a',M_partition,T_partition,T_basis,T_basis,number_of_trial_local_basis,number_of_test_local_basis,number_of_elements,matrix_size,Gauss_coefficient_reference_1D,Gauss_point_reference_1D,basis_type,1,basis_type,1);%组装载荷矢量
b = assemble_vector_from_1D_integral('function_f',M_partition,T_partition,T_basis,number_of_test_local_basis,number_of_elements,vector_size,Gauss_coefficient_reference_1D,Gauss_point_reference_1D,basis_type,0);%得到边界节点的信息矩阵
boundary_nodes = generate_boundary_nodes_1D(N_basis);%处理狄利克雷边界条件
[A,b]=treat_Dirichlet_boundary_1D('function_g',A,b,boundary_nodes,M_basis);%计算数值解
result=A\b;%计算所有节点上的最大误差
if basis_type==101h_basis=h_partition;
end
maxerror=get_maximum_error_1D(result,N_basis,left,h_basis);
maximum_error_at_all_nodes_of_FE = maxerrorend

三、组装刚度矩阵

assemble_matrix_from_1D_integral.m

function A = assemble_matrix_from_1D_integral(coefficient_function_name,P_mesh,T_mesh,T_basis_trial,T_basis_test,number_of_trial_local_basis,number_of_test_local_basis,number_of_elements,matrix_size,Gauss_point_number,trial_basis_type,trial_derivative_degree,test_basis_type,test_derivative_degree)
%% 注释
%{目的：计算刚度矩阵A从整个定义域上的一维积分组合出一个矩阵一维积分的被积函数必须为以下格式:a coefficient function * a trial FE basis function(or its derivatives) * a test FE basis function (or its derivatives).系数函数 *  试探函数（或其导数） * 测试函数（或其导数）
------------------------------------------------------------------------------------------------------------------------
Input:coefficient_function_name   ---------- 系数函数 c(x)P_mesh   ---------- 存储等分区间中所有节点的坐标T_mesh   ---------- 存储等分区间中每个单元节点的全局索引T_basis_trial   ---------- 试探(trial function)基函数的节点全局索引T_basis_test   ---------- 测试(test function)基函数的节点全局索引number_of_trial_local_basis   ---------- 局部试探函数的局部FE基函数的个数number_of_test_local_basis   ---------- 局部测试函数的局部FE基函数的个数number_of_elements   ---------- 网格划分的单元个数matrix_size   ---------- 总刚的大小  存储总刚的行数和列数Gauss_point_number   ---------- 高斯点数trial_basis_type   ---------- 试探FE基函数的类型trial_derivative_degree   ---------- 试探FE基函数的导数阶test_basis_type   ---------- 测试FE基函数的类型test_derivative_degree   ---------- 测试FE基函数的导数阶Output:A  ------- 总体刚度矩阵
------------------------------------------------------------------------------------------------------------------------------注：standard_GaussWeight,standard_GaussPoint:参考区间[-1,1]上的高斯系数和高斯点trial_basis_type:试用FE基函数的类型。trial_vertices: 试探函数单元的所有顶点的坐标trial_basis_type:试验FE基函数的类型trial_basis_index: 有限元试探基函数的指标，以指定我们要使用哪个有限元试探基函数trial_derivative_degree: the trial FE basis的导数阶test_vertices、test_basis_type、test_basis_index、test_derivative_degree:与trial_basis类似vertices: 一维单元上的两个顶点的坐标Gauss_coefficient_local_1D,Gauss_point_local_1D: 局部区间上的高斯系数和高斯点------------------------------------------------------------------------------------------------------------------------------ 孟伟, 大连理工大学- 1475207248@qq.com / mw21933005@mail.dlut.edu.cn
------------------------------------------------------------------------------------------------------------------------------
%}
%% 程序
%{
思路：遍历所有单元在每个元素上，计算试探和测验FE基函数的所有可能组合的非零积分将这些非零积分的值组合到刚度矩阵中
%}A=sparse(matrix_size(1),matrix_size(2)); %总刚 初始化
% 高斯积分在标准高斯区间[-1,1]上的高斯点和权重
[standard_GaussWeight,standard_GaussPoint] = get_standard_gauss_1D(Gauss_point_number);for n=1:number_of_elements    %遍历所有单元vertices = P_mesh(:,T_mesh(:,n));  %第n个单元的节点坐标   T_mesh(:,n)取出节点的全局编号矩阵T的第n列lower_bound=min(vertices(1),vertices(2));upper_bound=max(vertices(1),vertices(2));%  Gauss_coefficient_local_1D,Gauss_point_local_1D： 局部区间对应 区间[-1,1] 生成的新高斯系数和高斯点[local_GaussWeight_1D,local_GaussPoint_1D] = get_Gauss_local_1D(standard_GaussWeight,standard_GaussPoint,lower_bound,upper_bound);%  计算单元上非零积分 并将其组装到刚度矩阵A中for alpha=1:number_of_trial_local_basis    for beta=1:number_of_test_local_basis      % Gauss_quadrature_for_1D_integral_trial_test()函数   用高斯积分公式计算[xn,xn+1]上的一维积分temp = Gauss_quadrature_for_1D_integral_trial_test(coefficient_function_name,local_GaussWeight_1D,local_GaussPoint_1D,vertices,trial_basis_type,alpha,trial_derivative_degree,vertices,test_basis_type,beta,test_derivative_degree); %i = T_basis_test(beta,n);%j = T_basis_trial(alpha,n);%A(i,j) = A(i,j) + temp;A(T_basis_test(beta,n),T_basis_trial(alpha,n)) = A(T_basis_test(beta,n),T_basis_trial(alpha,n))+temp;endendend

2.1 算法

在这里插入图片描述

2.2 get_standard_gauss_1D.m

function [standard_GaussWeight,standard_GaussPoint] = get_standard_gauss_1D(Gauss_point_number)
%% 
%{
-------------------------------------------------------------------------------------------------得到标准高斯区间[-1,1]上的高斯点和高斯权重Input:Gauss_point_number ---------- 高斯点个数Output:standard_GaussWeight  ------- 标准高斯区间的高斯权重standard_GaussPoint   ------- 标准高斯区间的高斯点---------------------------------------- 孟伟, 大连理工大学- 1475207248@qq.com / mw21933005@mail.dlut.edu.cn
-------------------------------------------------------------------------------------------------
%}
%% 
if Gauss_point_number==2standard_GaussWeight=[1,1];standard_GaussPoint=[-1/sqrt(3),1/sqrt(3)];
elseif Gauss_point_number==4standard_GaussWeight=[0.3478548451,0.3478548451,0.6521451549,0.6521451549];standard_GaussPoint=[0.8611363116,-0.8611363116,0.3399810436,-0.3399810436];
elseif Gauss_point_number==8standard_GaussWeight=[0.1012285363,0.1012285363,0.2223810345,0.2223810345,0.3137066459,0.3137066459,0.3626837834,0.3626837834];standard_GaussPoint=[0.9602898565,-0.9602898565,0.7966664774,-0.7966664774,0.5255324099,-0.5255324099,0.1834346425,-0.1834346425];endend

2.3 get_Gauss_local_1D.m

function [local_GaussWeight_1D,local_GaussPoint_1D] = get_Gauss_local_1D(standard_GaussWeight_1D,standard_GaussPoint_1D,lower_bound,upper_bound)
%{目的： 使用仿射变换在任意区间[lower_bound,upper_bound]上生成高斯系数和高斯点由[lower_bound,upper_bound]区间 转换为 高斯积分标准区间[-1,1]上的高斯系数和高斯点
--------------------------------------------------------------------------------------Input:lower_bound  ----------------- 积分下限upper_bound  ----------------- 积分上限standard_GaussWeight  -------- 标准高斯区间[-1,1]的高斯权重standard_GaussPoint   -------- 标准高斯区间[-1,1]的高斯点Output:local_GaussWeight_1D   -------  任意区间上对应的高斯系数local_GaussPoint_1D    -------  任意区间上对应的高斯点
--------------------------------------------------------------------------------------注1：[a,b] 区间的gauss积分 ---> [-1,1] 区间的gauss积分   [-1,1] gauss积分系数 * (b-a)/2[-1,1] gauss积分点 * (b-a)/2 + (b+a)/2local_GaussWeight_1D：存储的是一个向量
--------------------------------------------------------------------------------------- 孟伟, 大连理工大学- 1475207248@qq.com / mw21933005@mail.dlut.edu.cn
--------------------------------------------------------------------------------------
%}local_GaussWeight_1D = (upper_bound-lower_bound) * standard_GaussWeight_1D / 2;
local_GaussPoint_1D =(upper_bound-lower_bound) * standard_GaussPoint_1D / 2 + (upper_bound + lower_bound)/ 2;
end

程序太多了，我就不一一放上去了，可以从我的博客有限元基础编程-何晓明老师课件-一维程序实现matlab下载

最后非常感谢何晓明老师无私的讲解


http://www.ppmy.cn/news/40919.html

相关文章

基于5G技术的智能导航机器人及AR巡逻应用开发项目实施方案（下）

基于5G技术的智能导航机器人及AR巡逻应用开发项目实施方案（下）

稳定性测试
阅读更多...
软件测试别再说简历项目不会写了，给你安排的明明白白

软件测试别再说简历项目不会写了，给你安排的明明白白

目录 个人信息 职业技能 工作经历 项目经历 工作经历 项目经历 教育经历 自我评价 个人信息 姓 名：xxx 性 别：女 手 机：xxxxxxxxxxxx 最高学历：统招硕士 工作年限：3 年 职…
阅读更多...
docker镜像深入学习，docker镜像发布公有云与私有云

docker镜像深入学习，docker镜像发布公有云与私有云

文章目录一、docker镜像概述1、什么是docker镜像2、镜像的分层3、什么是UnionFS（联合文件系统）4、 Docker镜像加载原理5、为什么 Docker 镜像要采用这种分层结构二、Docker镜像commit操作1、使用示例三、本地镜像发布到阿里云1、创建一个新的镜像2、在阿…
阅读更多...
fl studio 21打不开，FL工程文件也打不开怎么办？

fl studio 21打不开，FL工程文件也打不开怎么办？

FL Studio 21全称Fruity Loops Studio，就是大家熟悉的水果编曲软件，一个全能的音乐制作软件，包括编曲、录音、剪辑和混音等诸多功能，让你的电脑编程一个全能的录音室。FL Studio 21版本发布了，为我们带来了多种新功能&…
阅读更多...
正版软件 Microsoft 365 家庭版 1用户 58

正版软件 Microsoft 365 家庭版 1用户 58

Office 365 是一个基于云的订阅服务，汇集了当今人们工作中使用的最佳工具。通过将 Excel 和 Outlook 等一流应用与 OneDrive 和 Microsoft Teams 等强大的云服务相结合，Office 365 可让任何人使用任何设备随时随地创建和共享内容。 内附应用：…
阅读更多...
【ROS2指南-1】配置ROS2环境

【ROS2指南-1】配置ROS2环境

资料来源Configuring your ROS 2 environment — ROS 2 Documentation: Dashing documentationhttp://docs.ros.org/en/dashing/Tutorials/Configuring-ROS2-Environment.html 目标：本教程将向您展示如何准备 ROS 2 环境。 教程级别：初学者 时间&…
阅读更多...
【JAVA真的没出路了吗？】

【JAVA真的没出路了吗？】

2023年了，转行IT学习Java是不是已经听过看过很多次了。随之而来的类似学Java没出路、Java不行了、对Java感到绝望等等一系列的制造焦虑的话题也在网上层出不穷，席卷了一大片的对行业不了解的吃瓜群众或是正在学习中的人。如果是行外人真的会被这种言论轻…
阅读更多...
Flink CDC入门案例

Flink CDC入门案例

由于Flink CDC是基于日志的方式，因此需要开启MySQL的binlog日志。 开启binlog日志的配置如下 #1.编辑MySQL的配置文件 vim /etc/my.cnf #添加如下内容 [mysqld] log-binmysql-bin # 开启 binlog binlog-formatROW # 选择 ROW 模式 server_id1 # 配置 MySQL replact…
阅读更多...
最新文章

Copyright @ 2022~2023