使用太极taichi写一个只有一个三角形的有限元
公式来源
https://blog.csdn.net/weixin_43940314/article/details/128935230
GAME103
https://games-cn.org/games103-slides/
初始化我们的三角形
全局的坐标范围为0-1
我们的三角形如图所示

@ti.kernel
def init():X[0] = [0.5, 0.5]X[1] = [0.5, 0.6]X[2] = [0.6, 0.5]x[0] = X[0] + [0, 0.01]x[1] = X[1]x[2] = X[2]
X是rest pos
x是current pos
这里给一个小的增量是为了看出来被拉了,否则产生不了弹性力
公式抄录
[f1f2]=−ArefFS[X10X20]−T\begin{bmatrix} \mathbf{f_1} & \mathbf{f_2} \end{bmatrix}= -A^{ref} \mathbf{F} \mathbf{S } \begin{bmatrix} \mathbf{X_{10}} & \mathbf{X_{20}} \end{bmatrix}^{-T} [f1f2]=−ArefFS[X10X20]−T
F=[x10x20][X10X20]−1F=\begin{bmatrix} x_{10} & x_{20} \end{bmatrix}\begin{bmatrix} X_{10} & X_{20} \end{bmatrix}^{-1} F=[x10x20][X10X20]−1
S=2μG+λtrace(C)IS = 2 \mu G + \lambda trace(C) I S=2μG+λtrace(C)I
G=12(FTF−I)G = \frac{1}{2} (F^T F -I) G=21(FTF−I)
0. 设定一下材料参数
dim=2
n_particles = 3
n_elements = 1
area = 0.1*0.1*0.5
dt = 1e-4
E, nu = 1e3, 0.33 # Young's modulus and Poisson's ratio
mu, lam = E / 2 / (1 + nu), E * nu / (1 + nu) / (1 - 2 * nu) # Lame parameters
1 计算F
根据上面的公式,我们要先算F
@ti.kernel
def substep():#compute deformation gradientfor i in range(n_elements):Dm =ti.Matrix([[x[1][0]-x[0][0], x[2][0]-x[0][0]], [x[1][1]-x[0][1], x[2][1]-x[0][1]]])Dm_inv[i] = Dm.inverse()Ds = ti.Matrix([[X[1][0]-X[0][0], X[2][0]-X[0][0]], [X[1][1]-X[0][1], X[2][1]-X[0][1]]])F[i] = Ds @ Dm_inv[i]
2 计算格林应变
#compute green strain
for i in range(n_elements):G[i] = 0.5 * (F[i].transpose() @ F[i] - ti.Matrix([[1, 0], [0, 1]]))
3 计算PK1
#compute second Piola Kirchhoff stress
for i in range(n_elements):S[i] = 2 * mu *G[i] + lam * (G[i][0,0]+G[i][1,1]) * ti.Matrix([[1, 0], [0, 1]])
4 计算粒子上的力
#compute force(先暂且就计算一个三角形的力,后面再考虑多个三角形的情况)
force_matrix = F[0] @ S[0] @ Dm_inv[0].transpose() * area
force[1] = ti.Vector([force_matrix[0, 0], force_matrix[1, 0]])
force[2] = ti.Vector([force_matrix[0, 1], force_matrix[1, 1]])
force[0] = -force[1] - force[2]
5 加个重力
#gravityfor i in range(n_particles):force[i][1] -= 0.1
6 时间积分 同时处理边界条件
#time integration(with boundary condition)eps = 0.01for i in range(n_particles):vel[i] += dt * force[i]#boundary conditioncond = (x[i] < eps) & (vel[i] < 0) | (x[i] > 1) & (vel[i] > eps)for j in ti.static(range(dim)):if cond[j]:vel[i][j] = 0 x[i] += dt * vel[i]
完整的程序
# ref: https://blog.csdn.net/weixin_43940314/article/details/128935230import taichi as ti
import numpy as npti.init(arch=ti.cpu, debug=True)dim=2
n_particles = 3
n_elements = 1
area = 0.1*0.1*0.5
# lam = 1
# mu = 1
dt = 1e-4
E, nu = 1e3, 0.33 # Young's modulus and Poisson's ratio
mu, lam = E / 2 / (1 + nu), E * nu / (1 + nu) / (1 - 2 * nu) # Lame parametersx = ti.Vector.field(dim, dtype=float, shape=n_particles) #deformed position
force = ti.Vector.field(dim, dtype=float, shape=n_particles)
vel = ti.Vector.field(dim, dtype=float, shape=n_particles)
X = ti.Vector.field(dim, dtype=float, shape=n_particles) #undeformed position
S = ti.Matrix.field(n=dim, m=dim, dtype=float, shape=n_elements) #Second Piola Kirchhoff stress
F = ti.Matrix.field(n=dim, m=dim, dtype=float, shape=n_elements) #deformation gradient
G = ti.Matrix.field(n=dim, m=dim, dtype=float, shape=n_elements) #green strain@ti.kernel
def init():X[0] = [0.5, 0.5]X[1] = [0.5, 0.6]X[2] = [0.6, 0.5]x[0] = X[0] + [0, 0.01]x[1] = X[1]x[2] = X[2]Dm_inv = ti.Matrix.field(n=dim, m=dim, dtype=float, shape=n_elements)
@ti.kernel
def substep():#compute deformation gradientfor i in range(n_elements):Dm =ti.Matrix([[x[1][0]-x[0][0], x[2][0]-x[0][0]], [x[1][1]-x[0][1], x[2][1]-x[0][1]]])Dm_inv[i] = Dm.inverse()Ds = ti.Matrix([[X[1][0]-X[0][0], X[2][0]-X[0][0]], [X[1][1]-X[0][1], X[2][1]-X[0][1]]])F[i] = Ds @ Dm_inv[i]# print(F[0])#compute green strainfor i in range(n_elements):G[i] = 0.5 * (F[i].transpose() @ F[i] - ti.Matrix([[1, 0], [0, 1]]))#compute second Piola Kirchhoff stressfor i in range(n_elements):S[i] = 2 * mu *G[i] + lam * (G[i][0,0]+G[i][1,1]) * ti.Matrix([[1, 0], [0, 1]])#compute force(先暂且就计算一个三角形的力,后面再考虑多个三角形的情况)force_matrix = F[0] @ S[0] @ Dm_inv[0].transpose() * areaforce[1] = ti.Vector([force_matrix[0, 0], force_matrix[1, 0]])force[2] = ti.Vector([force_matrix[0, 1], force_matrix[1, 1]])force[0] = -force[1] - force[2]# print(force[0])#gravityfor i in range(n_particles):force[i][1] -= 0.1#time integration(with boundary condition)eps = 0.01for i in range(n_particles):vel[i] += dt * force[i]#boundary conditioncond = (x[i] < eps) & (vel[i] < 0) | (x[i] > 1) & (vel[i] > eps)for j in ti.static(range(dim)):if cond[j]:vel[i][j] = 0 x[i] += dt * vel[i]def main():init()gui = ti.GUI('my', (1024, 1024))while gui.running:for e in gui.get_events():if e.key == gui.ESCAPE:gui.running = Falseelif e.key == 'r':init()for i in range(30):substep()vertices_ = np.array([[0, 1, 2]], dtype=np.int32)particle_pos = x.to_numpy()a = vertices_.reshape(n_elements * 3)b = np.roll(vertices_, shift=1, axis=1).reshape(n_elements * 3)gui.lines(particle_pos[a], particle_pos[b], radius=1, color=0x4FB99F)gui.circles(particle_pos, radius=5, color=0xF2B134)gui.show()if __name__ == '__main__':main()

相关文章:
使用太极taichi写一个只有一个三角形的有限元
公式来源 https://blog.csdn.net/weixin_43940314/article/details/128935230 GAME103 https://games-cn.org/games103-slides/ 初始化我们的三角形 全局的坐标范围为0-1 我们的三角形如图所示 ti.kernel def init():X[0] [0.5, 0.5]X[1] [0.5, 0.6]X[2] [0.6, 0.5]x[0…...
进程,线程
进程是操作系统分配资源的基本单位,线程是CPU调度的基本单位。 PCB:进程控制块,操作系统描述程序的运行状态,通过结构体task,struct{…},统称为PCB(process control block)。是进程管理和控制的…...
第03章_基本的SELECT语句
第03章_基本的SELECT语句 讲师:尚硅谷-宋红康(江湖人称:康师傅) 官网:http://www.atguigu.com 1. SQL概述 1.1 SQL背景知识 1946 年,世界上第一台电脑诞生,如今,借由这台电脑发展…...
干货 | 简单了解运算放大器...
运算放大器发明至今已有数十年的历史,从最早的真空管演变为如今的集成电路,它在不同的电子产品中一直发挥着举足轻重的作用。而现如今信息家电、手机、PDA、网络等新兴应用的兴起更是将运算放大器推向了一个新的高度。01 运算放大器简述运算放大器&#…...
C++定位new用法及注意事项
使用定位new创建对象,显式调用析构函数是必须的,这是析构函数必须被显式调用的少数情形之一!, 另有一点!!!析构函数的调用必须与对象的构造顺序相反!切记!!&a…...
【Android笔记75】Android之翻页标签栏PagerTabStrip组件介绍及其使用
这篇文章,主要介绍Android之翻页标签栏PagerTabStrip组件及其使用。 目录 一、PagerTabStrip翻页标签栏 1.1、PagerTabStrip介绍 1.2、PagerTabStrip的使用 (1)创建布局文件...
【Kafka】【二】消息队列的流派
消息队列的流派 ⽬前消息队列的中间件选型有很多种: rabbitMQ:内部的可玩性(功能性)是⾮常强的rocketMQ: 阿⾥内部⼀个⼤神,根据kafka的内部执⾏原理,⼿写的⼀个消息队列中间 件。性能是与Kaf…...
现代 cmake (cmake 3.x) 操作大全
cmake 是一个跨平台编译工具,它面向各种平台提供适配的编译系统配置文件,进而调用这些编译系统完成编译工作。cmake 进入3.x 版本,指令大量更新,一些老的指令开始被新的指令集替代,并加入了一些更加高效的指令/参数。本…...
how https works?https工作原理
简单一句话: https http TLShttps 工作原理:HTTPS (Hypertext Transfer Protocol Secure)是一种带有安全性的通信协议,用于在互联网上传输信息。它通过使用加密来保护数据的隐私和完整性。下面是 HTTPS 的工作原理:初始化安全会…...
Docker的资源控制管理
目录 一、CPU控制 1、设置CPU使用率上限 2、设置CPU资源占用比(设置多个容器时才有效) 3、设置容器绑定指定的CPU 二、对内存使用进行限制 1、创建指定物理内存的容器 2、创建指定物理内存和swap的容器 3、 对磁盘IO配额控制(blkio&a…...
MMSeg无法使用单类自定义数据集训练
文章首发及后续更新:https://mwhls.top/4423.html,无图/无目录/格式错误/更多相关请至首发页查看。 新的更新内容请到mwhls.top查看。 欢迎提出任何疑问及批评,非常感谢! 摘要:将三通道图像转为一通道图像,…...
Redis使用方式
一、Redis基础部分: 1、redis介绍与安装比mysql快10倍以上 *****************redis适用场合**************** 1.取最新N个数据的操作 2.排行榜应用,取TOP N 操作 3.需要精确设定过期时间的应用 4.计数器应用 5.Uniq操作,获取某段时间所有数据排重值 6.实时系统,反垃圾系统7.P…...
无主之地3重型武器节奏评分榜(9.25) 枪械名 红字效果 元素属性 清图评分 Boss战评分 泛用性评分 特殊性评分 最终评级 掉落点 掉率 图片 瘟疫传播
无主之地3重型武器节奏评分榜(9.25) 枪械名 红字效果 元素属性 清图评分 Boss战评分 泛用性评分 特殊性评分 最终评级 掉落点 掉率 图片 瘟疫传播者 发射巨大能量球,能量球会额外生成追踪附近敌人的伴生弹 全属性 SSS SSS SSS - T0 伊甸6号-…...
什么是编程什么是算法
1.绪论 编程应在一个开发环境中完成源程序的编译和运行。首先,发现高级语言开发环境,TC,Windows系统的C++,R语言更适合数学专业的学生。然后学习掌握编程的方法,在学校学习,有时间的人可以在网上学习,或者购买教材自学。最后,编写源程序,并且在开发环境中实践。 例如…...
【c++】函数
文章目录函数的定义函数的调用值传递常见样式函数的声明函数的分文件编写函数的作用: 将一段经常使用的代码封装起来,减少重复代码。 一个较大的程序,一般分为若干个程序块,每个模板实现特定的功能。 函数的定义 返回值类型 函数…...
[golang gin框架] 1.Gin环境搭建,程序的热加载,路由GET,POST,PUT,DELETE
一.Gin 介绍Gin 是一个 Go (Golang) 编写的轻量级 http web 框架,运行速度非常快,如果你是性能和高效的追求者,推荐你使用 Gin 框架.Gin 最擅长的就是 Api 接口的高并发,如果项目的规模不大,业务相对简单,这…...
【开源】祁启云网络验证系统V1.11
简介 祁启云免费验证系统 一个使用golang语言、Web框架beego、前端Naive-Ui-Admin开发的免费网络验证系统 版本 当前版本1.11 更新方法 请直接将本目录中的verification.exe/verification直接覆盖到你服务器部署的目录,更新前,请先关闭正在运行的验…...
震源机制(Focal Mechanisms)之沙滩球(Bench Ball)
沙滩球包含如下信息: a - 判断断层类型,可根据球的颜色快速判断 b - 判断断层的走向(strike),倾角(dip) c - 确定滑移角/滑动角(rake) 走向 ,倾角,滑移角 如不了解断层的定义,可以先阅读:震…...
C++入门:多态
多态按字面的意思就是多种形态。当类之间存在层次结构,并且类之间是通过继承关联时,就会用到多态。C 多态意味着调用成员函数时,会根据调用函数的对象的类型来执行不同的函数。1、纯虚函数声明如下: virtual void funtion1()0; 纯…...
华为OD真题_工位序列统计友好度最大值(100分)(C++实现)
题目描述 工位由序列F1,F2…Fn组成,Fi值为0、1或2。其中0代表空置,1代表有人,2代表障碍物。 1、某一空位的友好度为左右连续老员工数之和 2、为方便新员工学习求助,优先安排友好度高的空位 给出工位序列,求所有空位中友好度的最大值。 输入描述 第一行为工位序列:F1,F…...
泛微Ecology流程数据查询避坑指南:workflow_currentoperator表里isremark字段到底怎么用?
泛微Ecology流程数据查询实战:解密workflow_currentoperator表关键字段 在泛微Ecology系统的二次开发过程中,流程数据的精准查询往往是开发者面临的第一道门槛。特别是当需要对接第三方系统或构建定制化报表时,对workflow_currentoperator表中…...
深入TC397与TLF35584的SPI通信:从寄存器操作到汽车ECU低功耗状态管理实战
深入TC397与TLF35584的SPI通信:从寄存器操作到汽车ECU低功耗状态管理实战 在汽车电子领域,电源管理芯片的选择与配置直接关系到整车电子控制单元(ECU)的可靠性与能耗表现。英飞凌的TLF35584作为一款高集成度电源管理IC,…...
Ostrakon-VL-8B模型剪枝与量化入门:降低部署资源消耗
Ostrakon-VL-8B模型剪枝与量化入门:降低部署资源消耗 想让大模型在普通电脑上跑起来?这听起来像是个遥不可及的梦想,尤其是对于Ostrakon-VL-8B这种参数规模不小的视觉语言模型。它功能强大,但随之而来的就是对GPU显存和算力的高要…...
Lingbot-Depth-Pretrain-VitL-14处理复杂光照与反射场景效果展示
Lingbot-Depth-Pretrain-VitL-14处理复杂光照与反射场景效果展示 深度估计技术,简单来说就是让计算机像人眼一样,判断出画面中每个物体离我们有多远。这项技术在自动驾驶、机器人导航、增强现实等领域都扮演着关键角色。然而,当场景中出现一…...
OpenClaw调试技巧:nanobot镜像的日志分析与问题定位
OpenClaw调试技巧:nanobot镜像的日志分析与问题定位 1. 为什么需要关注OpenClaw日志 上周我在本地部署nanobot镜像时遇到一个诡异现象:OpenClaw能正常接收飞书消息,但执行自动化任务时总在"思考阶段"卡住。这个问题困扰了我两天&…...
Obsidian移动端深度评测:安卓/iOS同步技巧+5个必装生产力插件
Obsidian移动端深度评测:安卓/iOS同步技巧5个必装生产力插件 在移动办公场景下,Obsidian作为一款强大的知识管理工具,其跨平台能力与插件生态为商务人士和学生群体提供了独特的价值。本文将深入解析Obsidian在Android和iOS平台的核心差异&…...
跨平台软件兼容方案全解析:从痛点到完美体验的技术实践
跨平台软件兼容方案全解析:从痛点到完美体验的技术实践 【免费下载链接】deepin-wine 【deepin源移植】Debian/Ubuntu上最快的QQ/微信安装方式 项目地址: https://gitcode.com/gh_mirrors/de/deepin-wine 在数字化办公与娱乐日益融合的今天,跨平台…...
GIS小白必看!Global Mapper处理正射影像的5个高频问题解答(含奥维地图导入避坑指南)
GIS新手实战指南:Global Mapper正射影像处理全解析 第一次打开Global Mapper时,那些密密麻麻的工具栏和复杂的参数设置确实让人望而生畏。去年我刚接触GIS时,处理无人机航拍的正射影像就踩了不少坑——坐标系选错导致影像偏移几百米、导出分幅…...
程序员转行学习 AI 大模型: 提示词工程 | 附精选学习资料
本文是程序员转行学习AI大模型的第12个核心知识点笔记,笔记后附精选的提示词工程学习资料。 当前阶段:还在学习知识点,由点及面,从 0 到 1 搭建 AI 大模型知识体系中。 系列更新,关注我,后续会持续记录分享…...
实战指南:基于快马生成电商订单自动化n8n工作流,无缝衔接shopify与crm
实战指南:基于快马生成电商订单自动化n8n工作流,无缝衔接shopify与crm 最近在帮朋友优化他们电商业务的后台流程,发现手动处理订单实在太费时间了。特别是遇到大促期间,订单量暴增,人工操作不仅效率低还容易出错。于是…...
