文章目录
前言
只是为方便学习,不做其他用途,课程理论学习来自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==101
N_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==101
number_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==101
h_basis=h_partition;
end
maxerror=get_maximum_error_1D(result,N_basis,left,h_basis);
maximum_error_at_all_nodes_of_FE = maxerror
end
二、一维有限元求解程序-框架
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==101
N_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==101
number_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==101
h_basis=h_partition;
end
maxerror=get_maximum_error_1D(result,N_basis,left,h_basis);
maximum_error_at_all_nodes_of_FE = maxerror
end
三、组装刚度矩阵
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;
end
end
end
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==2
standard_GaussWeight=[1,1];
standard_GaussPoint=[-1/sqrt(3),1/sqrt(3)];
elseif Gauss_point_number==4
standard_GaussWeight=[0.3478548451,0.3478548451,0.6521451549,0.6521451549];
standard_GaussPoint=[0.8611363116,-0.8611363116,0.3399810436,-0.3399810436];
elseif Gauss_point_number==8
standard_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];
end
end
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)/2
local_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下载
最后非常感谢何晓明老师无私的讲解