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

news2025/1/25 1:31:43

文章目录

  • 前言
  • 一、主程序
  • 二、一维有限元求解程序-框架
  • 三、组装刚度矩阵
    • 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==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下载

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

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/413614.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

RT-Thread线程管理以及内核裁剪

RT-Thread线程管理以及内核裁剪 1. RTOS概述 1.1 RTOS的定义 实时操作系统(Real-time operating system, RTOS),又称即时操作系统,它会按照排序运行、管理系统资源,并为开发应用程序提供一致的基础。 实时操作系统与…

核心业务2:借款人申请借款额度

核心业务2:借款人申请借款额度 1.业务流程图 ------------截止提交个人信息部分-------- 2.借款人申请借款额度数据库设计 3.借款人申请额度流程 4.前端代码逻辑 5.后端代码逻辑 ------------截止提交个人信息部分-------- 核心业务2:借款人申请借…

Linux从命令行管理文件

目录 一、创建链接文件 二、目录操作命令 1. 创建目录(make directory) 2. 统计目录及文件的空间占用情况 3. 删除目录文件 三、创建、删除普通文件 文件命名规则: (1)不能使用/来当文件名,/是用来做…

【WCH】CH32F203软件I2C驱动SSD1306 OLED

【WCH】CH32F203软件I2C驱动SSD1306 OLED📌相关篇《【WCH】CH32F203硬件I2C驱动SSD1306 OLED》📺驱动显示效果: 🌿OLED屏幕:i2c ssd1306 oled🔖驱动单片机型号:CH32F203 ✨由于CH32F203主频为96…

wordpress下载插件,安装失败,无法创建目录问题

刚开始安装这个wordpress,在发表文章时候想要在其中加上图片,不想一个个手动上传媒体库,耽误时间,然后就去下了个imagepaste这个复制粘贴的插件,当我打开安装插件搜索到的时候准备安装,尼玛出现“安装失败&…

若一个单词被拆分成多少token, word_ids得到的序号是相同的?还是序号累加的?

目录 问题描述: 问题实现: 方法一: 方法二: 问题描述: 在使用tokenizer进行编码的时候,经常会存在word被拆分成多个token的情况,不同的参数设置,会得到不同的结果。总的来说&…

redis——使用

session缓存缓存更新方式删除缓存vs更新缓存缓存和数据库操作原子性缓存和数据库操作顺序结论缓存问题缓存穿透缓存雪崩缓存击穿全局唯一ID数据并发线程安全单体分布式redis分布式锁的问题redis消息队列listpubsubstream消息推送session 问题:session存在tomcat服务…

【Linux驱动开发】023 platform设备驱动

一、前言 驱动分离目的:提高Linux代码重用性和可移植性。 二、驱动的分隔与分离 百度看了很多,大多都没讲清楚为什么使用platform驱动,为什么驱动分隔与分离可以提高代码重用性,只是在讲实现的结构体、函数接口等等&#xff0c…

npm、pnpm、yarn的常用命令

npm、pnpm、yarn的常用命令 文章目录npm、pnpm、yarn的常用命令一、常用命令1、npm命令2、pnpm命令:3、yarn命令二、对比一、常用命令 1、npm命令 npm init: 初始化一个新的npm包。 npm install: 安装项目依赖项。 npm install : 安装指定的包。 npm install --sa…

【Java数据结构】链表(Linked List)-双向链表

双向链表(Linked List)是一种常用的数据结构,它允许在所有节点中快速添加或删除元素,并且可以有效地实现反向遍历。本篇文章将介绍双向链表的基础知识,并提供使用Java语言实现该数据结构的示例代码。 一、双向链表的基…

mysql数据库事务脏读、不可重复度、幻读详解

文章目录1 事务隔离级别2 脏读3 不可重复度3.1 解决了脏读的问题。3.2 有不可重复度的问题4 幻读4.1 没有脏读和不可重复读的问题4.2 有幻读的问题5 serializable1 事务隔离级别 read-uncommitted:脏读、不可重复度、幻读,均可出现。安全性低&#xff0…

HBase架构篇 - Hadoop家族的天之骄子HBase

HBase的基本组成结构 表(table) HBase 的数据存储在表中。表名是一个字符串。表由行和列组成。 行(row) HBase 的行由行键(rowkey)和 n 个列(column)组成。行键没有数据类型&…

《花雕学AI》06:抢先体验ChatGPT的九个国内镜像站之试用与综合评测

最近ChatGPT持续大火,大家们是不是在网上看到各种和ChatGPT有趣聊天的截图,奈何自己实力不够,被网络拒之门外,只能眼馋别人的东西。看别人在体验,看别人玩,肯定不如自己玩一把舒服的啊。 上一期&#xff0…

2.5d风格的游戏模式如何制作

文章目录一、 介绍二、 绘制瓦片地图三、 添加场景物体,添加碰撞器四、 创建玩家五、 创建玩家动画六、 玩家脚本七、 2d转换成2.5d八、 “Q”键向左转动视角、“E”键向右转动视角九、 下载工程文件一、 介绍 制作一个类似饥荒风格的2.5d游戏模板。 2.5D游戏是指以…

Spring之循环依赖

什么事循环依赖 很简单的定义就是就如有两个对象A类,B类,其中两个类中的属性都有对方。 A类 public class A{private B b;}B类 public class B{ private A a; }在Spring中,什么情况下会出现循环依赖 如果要了解循环依赖,首先…

基于matlab进行雷达信号模拟

一、前言此示例说明如何将基本工具箱工作流应用于以下方案:假设有一个工作频率为 4 GHz 的各向同性天线。假设天线位于全局坐标系的原点。有一个目标,其非波动雷达横截面为0.5平方米,最初位于(7000,5000,0&…

Linux下使用ClamAV病毒查杀

一、介绍Clam AntiVirus 是一款 UNIX 下开源的 (GPL) 反病毒工具包,专为邮件网关上的电子邮件扫描而设计。该工具包提供了包含灵活且可伸缩的监控程序、命令行扫描程序以及用于自动更新数据库的高级工具在内的大量实用程序。该工具包的核心在于可用于各类场合的反病…

CompletableFuture使用详解(IT枫斗者)

CompletableFuture使用详解 简介 概述 CompletableFuture是对Future的扩展和增强。CompletableFuture实现了Future接口,并在此基础上进行了丰富的扩展,完美弥补了Future的局限性,同时CompletableFuture实现了对任务编排的能力。借助这项能力…

2023最新快速单机创建三主三从Redis集群

单机搭建Redis集群 本次采用Redis的5.0.14版本在单机centos8上搭建Redis三主三从集群. 1.创建6个文件夹 一个文件夹代表一个节点,同时也代表每个节点的端口号. 2.下载Redis文件并解压 使用命令: #下载Redis 可以将5.0.14替换成自己想要的版本 wget http://download.redis…

JavaScript面向对象编程再讲

JavaScript面向对象编程再讲 JavaScript支持的面向对象比较复杂,和其他编程语言又有其独特之处。本文是对以前博文 JavaScript的面向对象编程 https://blog.csdn.net/cnds123/article/details/109763357 补充。 概述 这部分是JavaScript面向对象的概括&#xff0c…