有限元基础编程-何晓明老师课件-一维程序实现matlab
资源内容介绍
有限元基础编程-何晓明老师课件-一维程序实现matlab %function result=poisson_solver_1D(left,right,h_partition,basis_type,Gauss_point_number)%% 改写程序时的测试程序 %{ ------------------------------------------------------------------------------------------------- Input: left ----------------- 问题求解区间左边界 right ----------------- 问题求解区间右边界 h_partition ------------ 均匀网格划分的步长 basis_type ------------- FE(有限元)基函数的类型 Gauss_point_number ----- 高斯点数 Output: result ---------------- 求解出来 Ax = b 的 x------------------------------------------------------------------------------------------------- 改编何晓明老师的求解一维泊松方程的有限元程序 注释： basis_type=101:1D 线性有限元 程序中用“FE”代替“finite element” 试探(trial)FE函数和测试(test)FE函数需要相同 到目前为止，该代码只能处理一维线性有限元和一维拉格朗日二次有限元。对于其他类型的FE，我们需要在代码中添加相应的信息 问题求解域是[left,right] h_basis 有限元节点的步长 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. ------------------------------------------------------------------------------------------------- - 孟伟, 大连理工大学 - 1475207248@qq.com / mw21933005@mail.dlut.edu.cn-------------------------------------------------------------------------------------------------%}%% 初始输入物理量 % 初始输入物理量 % format：设置输出格式 format short e：5字长浮点数format short e % FE(有限元)基函数的类型：二维线性有限元basis_type=101;% 问题域为 [left,right]*[bottom,top].left=0;right=1;% 高斯点数Gauss_point_number=4;% 均匀网格划分的步长h_partition = 1/8;%% 划分网格和有限元基函数相关信息N_partition = (right-left)/h_partition; %网格划分个数if basis_type == 101 N_basis = N_partition; %FE基函数个数end% 网格节点的坐标矩阵P 和 每个单元的节点的全局索引矩阵 T[P_mesh,T_mesh] = get_P_T_1D(left,right,h_partition,basis_type);if basis_type==101 P_basis = P_mesh; T_basis = T_mesh;end % 相关矩阵的大小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; %局部trial基函数个数 number_of_test_local_basis = 2; %局部test基函数个数end% 高斯积分在标准高斯区间[-1,1]上的高斯点和权重[standard_GaussWeight,standard_GaussPoint] = get_standard_gauss_1D(Gauss_point_number);%% 组装刚度矩阵A = assemble_matrix_from_1D_integral('function_c',P_mesh,T_mesh,T_basis,T_basis,number_of_trial_local_basis,number_of_test_local_basis,number_of_elements,matrix_size,Gauss_point_number,basis_type,1,basis_type,1);%组装载荷矢量b = assemble_vector_from_1D_integral('function_f',P_mesh,T_mesh,T_basis,number_of_test_local_basis,number_of_elements,vector_size,Gauss_point_number,basis_type,0);%得到边界节点的信息矩阵boundary_nodes = generate_boundary_nodes_1D(N_basis);%处理狄利克雷边界条件[A,b] = treat_Dirichlet_boundary_1D('function_g',A,b,boundary_nodes,P_basis);%计算数值解result = A\b;%计算所有节点上的最大误差if basis_type == 101 h_basis = h_partition; % 基函数划分的步长endmaxerror=get_maximum_error_1D(result,N_basis,left,h_basis);maximum_error_at_all_nodes_of_FE = maxerror%%%{ 生成 P： 网格节点的坐标矩阵 T： 每个单元的节点的全局索引矩阵 ------------------------------------------------------------------------------ Input: left ----------------- 问题求解区间左边界 right ----------------- 问题求解区间右边界 h_partition ------------ 均匀网格划分的步长 basis_type ------------- FE(有限元)基函数的类型 Output: P ------- 网格节点的坐标矩阵 P：1*(N+1) T ------- 每个单元的节点的全局索引矩阵 T： 2*N ------------------------------------------------------------------------------ 注1： 在一维中，M只有一行 ------------------------------------------------------------------------------ - 孟伟, 大连理工大学 - 1475207248@qq.com / mw21933005@mail.dlut.edu.cn------------------------------------------------------------------------------%}