CN112883609B - 实时模拟物体运动或形变的方法 - Google Patents

实时模拟物体运动或形变的方法 Download PDF

Info

Publication number
CN112883609B
CN112883609B CN202110153901.4A CN202110153901A CN112883609B CN 112883609 B CN112883609 B CN 112883609B CN 202110153901 A CN202110153901 A CN 202110153901A CN 112883609 B CN112883609 B CN 112883609B
Authority
CN
China
Prior art keywords
matrix
constraint
rigid body
internal force
energy
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110153901.4A
Other languages
English (en)
Other versions
CN112883609A (zh
Inventor
施若愚
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Shanghai Suoyi Intelligent Technology Co ltd
Original Assignee
Shanghai Suoyi Intelligent Technology Co ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Shanghai Suoyi Intelligent Technology Co ltd filed Critical Shanghai Suoyi Intelligent Technology Co ltd
Priority to CN202110153901.4A priority Critical patent/CN112883609B/zh
Publication of CN112883609A publication Critical patent/CN112883609A/zh
Application granted granted Critical
Publication of CN112883609B publication Critical patent/CN112883609B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Complex Calculations (AREA)

Abstract

一种实时模拟物体运动或形变的方法，包括：根据约束能计算约束能的梯度为该能量产生的内力为Fc＝‑GradCi×kiCi，‑Ci×ki为内力大小，GradCi为内力方向；将所有约束能的梯度组合成雅可比矩阵；令
Figure DDA0002933755830000011
得到Fc＝GradCi×λi，Ci＝αi×λi，物体处一点的应变能εtx E＝εt×D×ε，离散化到体积四面体后，得到一个体积四面体产生的约束能
Figure DDA0002933755830000012
将D分解为L‑1NL，得出
Figure DDA0002933755830000013
确定系统所有外力的和Fout；离散化后的物体运动微分方程求解公式(JM‑1Jt+Δt‑2*α)λ＝‑Δt‑1Jv‑JM‑1Fout‑C(Xn)Δt‑2，求解内力λ；更新系统的质点或刚体的位置、速度；上述方法实现几何约束能和有限元应变能的统一，同时引入的正则化方法使得物体的刚性系数到无穷大，兼容刚体模拟和软体模拟，使得模拟在统一框架下，达到很棒的视觉效果和稳定性。

Description

实时模拟物体运动或形变的方法
技术领域
本发明涉及模拟方法，特别涉及一种实时模拟物体运动或形变的方法。
背景技术
在计算机图形学领域，实时物理模拟，包括软体，刚体和流体是一个非常重要的方向。由于对运算速度和稳定性的极高要求，大量实时模型被提出：
主流的方法有：
各种非线性应变能的有限元分析方法O.C.Zienkiewicz and R.L.Taylor.<<TheFinite Element Method>>
基于极分解的应力变换的线性应变能有限元方法:M.Muller,<<InteractiveVirtual Materials>>
基于约束能的方法:M.Teschner<<A Versatile and Robust Model forGeometrically Complex Deformable Solids>>
基于约束的方法：M.Müller,B.Heidelberger,<<Position Based Dynamics>>
基于弹簧-质点的方法：<<Deformation constraints in a mass-spring modelto describe rigid cloth behavior>>
以上的模型又包括使用显式积分，隐式积分方法，半隐式积分方法。
使用显式积分的优点是速度快，计算复杂度低，缺点是系统容易进入不稳定状态。
隐式，半隐式积分方法则可以无条件稳定系统，但是却会引入而外的阻尼到系统中，缺点是公式推导复杂运算量也大。
由于稳定性的要求，现在都采用隐式积分或者半隐式积分方法。但是在这两种框架下以上的方法(有限元，约束能，质点弹簧)并不能互相兼容，例如基于有限元的方法虽然模拟效果更加真实，却很难引入使用几何约束能的力(抓取力，碰撞力，和其他刚体物件的交互等等)。
这是由于离散化后上述方法无法得到统一的线性方程组来进行整个解算过程，因此在引入碰撞约束和几何约束后并不能很好的融合如现有的有限元解算系统中。
此外，由于刚体的刚性系数是无穷大，而所有的上述方法都无法兼容无穷大的刚性系数，因此，刚体模拟算法很难和上述方法进行融合。
通常的解决方法是解算完软体后在结算刚体和其他约束。这样解算器分别在孤立的进行解算。使得模拟刚体-软体。或者软体-几何约束，碰撞等的时候效果非常不好。
发明内容
基于此，有必要提供一种可提高视觉效果的实时模拟物体运动或形变的方法。
同时，提供另一种可提高视觉效果的实时模拟物体运动或形变的方法。
一种实时模拟物体运动或形变的方法，包括：
确定约束内力方向：每个约束对系统产生的约束能为
Figure GDA0003467817310000021
根据约束能计算约束能的梯度为该能量产生的内力为Fc＝-GradCi×kiCi,其中，ki为约束Ci的刚性系数，k为约束C的刚性系数，GradCi为约束Ci关于整个系统坐标的梯度，所有约束势能之和产生的内力合力所有产生内力的约束能内力之和
Figure GDA0003467817310000022
变换成矩阵形式为
Figure GDA0003467817310000023
中J是为所有GradCi作为行向量所组成的雅可比阵，K是由每个约束的刚性系数组成的对角阵，
Figure GDA0003467817310000024
是系统所有约束组成的列向量，-Ci×ki为约束Ci的内力大小，GradCi为约束Ci的内力方向，J的每一列为一个内力的方向，
Figure GDA0003467817310000025
为内力大小；
确定所有约束能的梯度：将所有约束能的梯度组合成雅可比矩阵J(C1,C2....Ck)，其中C1、C2……Ck代表系统中的k个约束；
引入新变量：令
Figure GDA0003467817310000026
得到单个约束i作用在系统所有质点的内力Fc＝GradCi×λi，Ci＝αi×λi，系统所有质点或刚体受到的内力合力Fin＝∑GradCi×(-(Ci×ki))＝∑GradCi×λi＝Jtλ，
Figure GDA0003467817310000027
其中，GradCi代表约束Ci的梯度，ki为约束Ci的刚性系数，J为所有GradCi作为行向量所组成的雅可比阵，Jt为J的转置矩阵，K为由每个约束的刚性系数组成的对角矩阵，
Figure GDA0003467817310000031
为系统所有约束组成的列向量，λ为所有约束产生的内力大小、为所有λi组成的列向量，Matα由所有αi的对角矩阵，是个对角阵或者分块对角阵，它的对角元素为每个约束刚性系数ki的倒数αi＝1/ki，其他元素都为0的矩阵
Figure GDA0003467817310000032
离散化后的约束能：格林应变张量作为应变系数，应力张量E＝D×ε，其中D为由杨氏模量和泊松系数组成的弹性模量阵、为正定对阵矩阵，物体处一点的应变能εtxE＝εt×D×ε，离散化到体积四面体后，得到一个体积四面体产生的约束能是
Figure GDA0003467817310000033
其中V为四面体的体积，ε为格林应变张量，将D分解为L-1NL，得出
Figure GDA0003467817310000034
其中，N为对角矩阵，N是D的相似矩阵，N为D的特征值组成的对角矩阵，L为D的特征向量组成的矩阵；
确定外力：确定系统所有外力的和Fout；
求解内力：离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解内力大小
Figure GDA0003467817310000035
其中，v为系统当前所有质点速度组成的列向量，M-1代表所有质点的质量倒数组成的对角矩阵，Δt表示本模拟时间步长，C(Xn)为tn时刻系统的质点位置的3n维向量，或刚体速度和角速度组成的6n维向量，α为刚性系数矩阵；
更新位置、速度：根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt＝Vn+
M-1*(Jtλ+Fout)*Δt求出本时间步长模拟后的物体速度Vn+1，根据Xn+1＝Xn+Vn+1Δt求出本时间步长模拟完成后系统新的速度，根据求解出的速度、位置更新系统的质点或刚体的位置、速度，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置，Vn为本时间步长模拟前系统的速度，Xn本时间步长模拟前系统新的速度和位置，Fin代表系统由约束产生的内力。
在优选的实施例中，所述离散化后的约束能步骤中U等价于包含了6个约束表达式的约束能，N为6维对角矩阵、它是矩阵D对角化后得到的一个对角矩阵，所述矩阵
Figure GDA0003467817310000041
Figure GDA0003467817310000042
其中E为杨氏模量，μ为泊松系数。
在优选的实施例中，离散化后的约束能步骤中U等价于包含了6个约束表达式的约束能，N为6维对角矩阵、它是矩阵D对角化后得到的一个对角矩阵。
在优选的实施例中，若为软体则C为
Figure GDA0003467817310000043
所有约束能的梯度组合成雅可比矩阵J：
Figure GDA0003467817310000044
其中C1、C2……Ck代表系统中的所有的k个约束，xn为系统第n个质点的坐标、为一个三维向量。
在优选实施例中，刚体的约束C为C(P,Q)其中P为刚体的位置、为3维向量，Q为刚体的旋转四元数、为4维向量，刚体位置和速度的关系、以及刚体旋转四元数和角速度的关系如下：
Figure GDA0003467817310000051
对C(P,Q)求关于速度v、角速度ω的梯度值计算出作用在刚体上的力。
在优选实施例中，刚体上的点和系统中软体上的一个质点的固联约束为
C(P,Q,X)＝P+Q*R–X，其中P为刚体的位置，Q为刚体的旋转四元数，R为固联点在刚体本地坐标系下的向量，X为软体中某个质点得到位置，位置和速度的关系、以及旋转四元数和角速度的关系
Figure GDA0003467817310000052
Figure GDA0003467817310000053
求C(P,Q,X)＝P+Q*R–X关于速度v、角速度ω、以及质点位置的梯度值计算作用在刚体和质点上的力。
在优选实施例中，约束C(P,Q,X)＝||P+Q*R–X||在当前位置P(t)、Q(t)处做线性展开，由位置和速度的关系、以及旋转四元数和角速度的关系得到
Figure GDA0003467817310000054
Figure GDA0003467817310000055
Figure GDA0003467817310000056
Figure GDA0003467817310000057
Figure GDA0003467817310000058
写成向量形式：
Figure GDA0003467817310000059
同理
Figure GDA0003467817310000061
Figure GDA0003467817310000062
Figure GDA0003467817310000063
Figure GDA0003467817310000064
写成向量形式：
Figure GDA0003467817310000065
其中Pn是当前刚体位移，Qn是当前刚体旋转四元数，dt是模拟步长，ν是刚体线性速度，ω是角速度，对与上面的约束表达式，求它的梯度值计算作用在刚体和质点上的力，即求它关于速度ν、角速度ω、以及质点位置的梯度值计算式子作用在刚体和质点上的力。
在优选实施例中，对于离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，若系统有n个不同的约束C则
Figure GDA0003467817310000066
若系统有n个不同的约束C和m个应变能约束，则
Figure GDA0003467817310000067
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure GDA0003467817310000071
Figure GDA0003467817310000072
ki为约束Ci的刚性系数，系统有m个应变能约束对应有m个D对角块，
Figure GDA0003467817310000073
为第一个应变能约束对应的D的逆，
Figure GDA0003467817310000074
为第m个应变能约束对应的D的逆。
一种实时模拟物体运动或形变的方法，包括：
确定外力：确定系统所有外力的和Fout；
离散化的应变张量：对于任意一点的应变张量
Figure GDA0003467817310000075
G为形变梯度
Figure GDA0003467817310000076
其中x,y,z为物质空间某点的位置即未形变坐标，X、Y、Z为该点当前的位置坐标，将G离散化表示到四面体的四个顶点，得到G的离散表达式：G＝AB-1，其中矩阵A中的
Figure GDA0003467817310000077
Figure GDA0003467817310000078
为四面体第k个顶点的当前坐标，k＝0,1,2,3，
Figure GDA0003467817310000079
对于矩阵B，其中
Figure GDA00034678173100000710
Figure GDA00034678173100000711
Figure GDA00034678173100000712
为四面体第k个顶点的未形变坐标，k＝0,1,2,3，
Figure GDA00034678173100000713
确定应变能：有限元模拟中的应变能
Figure GDA00034678173100000714
令：
Figure GDA00034678173100000715
对Eu求梯度，所形成的内力为
Figure GDA00034678173100000716
VD*Cu为内力大小，V为四面体体积，ε＝(ε11 ε22 ε33 ε12 ε13 ε23)为6维的格林应变张量，D为杨氏模量和泊松系数组成的弹性模量阵，VD刚性系统矩阵为
Figure GDA00034678173100000717
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure GDA0003467817310000081
ki为约束Ci的刚性系数，系统有m个应变能约束对应有m个D对角块，
Figure GDA0003467817310000082
为第一个应变能约束对应的D的逆，
Figure GDA0003467817310000083
为第m个应变能约束对应的D的逆；
确定应变能约束的梯度：应变能约束的梯度向量为雅克比矩阵
Figure GDA0003467817310000084
求解内力：根据离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解内力大小λ＝VD*Cu，其中，v为系统当前所有质点速度组成的列向量，M-1代表所有质点的质量倒数组成的对角矩阵，Δt表示本模拟时间步长，α＝VD为刚性系统矩阵,C(Xn)为tn时刻系统的质点位置的3n维向量，或刚体速度和角速度组成的6n维向量；
更新位置、速度：根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt求出本时间步模拟后的物体速度Vn+1，根据Xn+1＝Xn+Vn+1Δt，其中，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置，Vn为本时间步长模拟前系统的速度，Xn本时间步长模拟前系统新的速度和位置，Fin代表系统由约束产生的内力，根据求解出的速度、位置更新系统的质点或刚体的位置、速度。
在优选的实施例中，所述应变张量
Figure GDA0003467817310000085
G为四面体的形变梯度的离散化表示G＝AB-1＝(A1A2A3)(c1c2c3)，B由四面体的顶点在未形变的时候的位置组成的矩阵，其逆矩阵B-1三个列向量为c1c2c3，由于四面体的体积总是>0，故B的行列式值等于B的混合积等于四面体体积，始终不等于0，故B-1始终存在，G＝(g1，g2，g3)，所以
Figure GDA0003467817310000086
其中δij＝1(i＝j)否则＝0；
根据G的表达式可知gi＝P*ci，gj＝P*cj
所以
Figure GDA0003467817310000091
Figure GDA0003467817310000092
同理
Figure GDA0003467817310000093
....
Figure GDA0003467817310000094
Figure GDA0003467817310000095
其中，ci为B-1的第i个列向量，cik为B-1的第i行第k列的元素，gi为G的第i个列向量。
在优选实施例中，
Figure GDA0003467817310000096
为六维应变张量的六个内力方向组成的矩阵，所述矩阵
Figure GDA0003467817310000097
Figure GDA0003467817310000098
其中E为杨氏模量，μ为泊松系数。
上述实时模拟物体运动或形变的方法，采用半隐式积分模拟方法，实现几何约束能和有限元应变能的统一，同时引入的正则化方法可以使得物体的刚性系数到无穷大，从而兼容刚体模拟和软体模拟。使得模拟在一个统一的框架下，达到了很棒的视觉效果和稳定性，同时并没有增加任何的运算开销。
具体实施方式
本发明一实施例的实时模拟物体运动或形变的方法，包括：
确定约束内力方向：每个约束对系统产生的约束能为
Figure GDA0003467817310000101
根据约束能计算约束能的梯度为该能量产生的内力为Fc＝-GradCi×kiCi,其中，ki为约束Ci的刚性系数，k为约束C的刚性系数，GradCi为约束Ci关于整个系统坐标的梯度，所有约束势能之和产生的内力合力是所有产生内力的约束能内力之和
Figure GDA0003467817310000102
变换成矩阵形式为
Figure GDA0003467817310000103
中J是为所有GradCi作为行向量所组成的雅可比阵，K是由每个约束的刚性系数组成的对角阵，
Figure GDA0003467817310000104
是系统所有约束组成的列向量，-Ci×ki为约束Ci的内力大小，GradCi为约束Ci的内力方向，J的每一列为一个内力的方向，
Figure GDA0003467817310000105
为内力大小；
确定所有约束能的梯度：将所有约束能的梯度组合成雅可比矩阵J(C1,C2....Ck)，其中C1、C2……Ck代表系统中的k个约束；
引入新变量：令
Figure GDA0003467817310000106
得到单个约束i作用在系统所有质点的内力Fc＝GradCi×λi，Ci＝αi×λi，系统所有质点或刚体受到的内力合力Fin＝∑GradCi×(-(Ci×ki))＝∑GradCi×λi＝Jtλ，
Figure GDA0003467817310000107
其中，GradCi代表约束Ci的梯度，ki为约束Ci的刚性系数，J为所有GradCi作为行向量所组成的雅可比阵，Jt为J的转置矩阵，K为由每个约束的刚性系数组成的对角矩阵，
Figure GDA0003467817310000108
为系统所有约束组成的列向量，λ为所有约束产生的内力大小、为所有λi组成的列向量，每个约束施加在系统所有质点上的力为
Figure GDA0003467817310000109
Matα由所有αi的对角矩阵，是个对角阵或者分块对角阵，它的对角元素为每个约束刚性系数ki的倒数αi＝1/ki，其他元素都为0的矩阵
Figure GDA00034678173100001010
离散化后的约束能：格林应变张量作为应变系数，应力张量E＝D×ε，其中D为由杨氏模量和泊松系数组成的弹性模量阵、为正定对阵矩阵，物体处一点的应变能εtxE＝εt×D×ε，离散化到体积四面体后，得到一个体积四面体产生的约束能是
Figure GDA0003467817310000111
其中V为四面体的体积，ε为格林应变张量，将D分解为L-1NL，其中，N为对角矩阵，得出
Figure GDA0003467817310000112
U等价于包含了6个约束表达式的约束能，N为6维对角矩阵它是矩阵D对角化后得到的一个对角矩阵；
确定外力：确定系统所有外力的和Fout；外力Fout可以根据力检测元件进行检测计算，如力传感器等；
求解内力：离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解内力大小
Figure GDA0003467817310000113
其中，v为系统当前所有质点速度组成的列向量，M-1代表所有质点的质量倒数组成的对角矩阵，Δt表示本模拟时间步长，C(Xn)为tn时刻系统的质点位置的3n维向量，或刚体速度和角速度组成的6n维向量；
更新位置、速度：根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt＝Vn+M-1*((Jtλ+Fout)*Δt求出本时间步长模拟后的物体速度Vn+1，根据Xn+1＝Xn+Vn+1Δt求出本时间步长模拟完成后系统新的速度，根据求解出的速度、位置更新系统的质点或刚体的位置、速度，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置，Vn为本时间步长模拟前系统的速度，Xn本时间步长模拟前系统新的速度和位置，Fin代表系统由约束产生的内力。
进一步，求解内力步骤前还包括确定刚性系数矩阵α，对于离散化后的物体运动微分方程求解公式：(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，若系统有n个不同的约束C则
Figure GDA0003467817310000121
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure GDA0003467817310000122
Figure GDA0003467817310000123
进一步，本实施例的矩阵D为一个6×6的对称矩阵，其中元素由杨氏模量和泊松系数组成：
矩阵
Figure GDA0003467817310000124
Figure GDA0003467817310000125
其中E为杨氏模量，μ为泊松系数。
若为软体则C为
Figure GDA0003467817310000126
所有约束能的梯度组合成雅可比矩阵J：
Figure GDA0003467817310000127
其中C1、C2……Ck代表系统中的所有的k个约束，xn为系统第n个质点的坐标、为一个三维向量。
每个约束对系统产生的约束能为
Figure GDA0003467817310000128
其中
Figure GDA0003467817310000129
是约束的表达式，k是该约束的刚性系数。约束能的梯度乘以-1就是该约束对系统中所有的质点产生的内力。所以对约束能
Figure GDA00034678173100001210
梯度(偏导数)
Figure GDA0003467817310000131
其中
Figure GDA0003467817310000132
称为内力的方向,
Figure GDA0003467817310000133
为内力的大小，
Figure GDA0003467817310000134
内力的方向就是约束的梯度向量，而大小就是刚性系数k乘以约束当前的值
Figure GDA0003467817310000135
若为刚体，可将刚体的约束C为写为C(P,Q)其中P为刚体的位置、为3维向量，Q为刚体的旋转四元数、为4维向量，刚体位置和速度的关系、以及刚体旋转四元数和角速度的关系如下：
Figure GDA0003467817310000136
对C(P,Q)求关于速度v、角速度ω的梯度值计算作用在刚体上的力。
具体的约束C是根据自己的需求定义，没有固定的形式。
如可将刚体上的点和系统中软体上的一个质点的固联约束为C(P,Q,X)＝P+Q*R–X，其中P为刚体的位置，Q为刚体的旋转四元数，R为固联点在刚体本地坐标系下的向量，X为软体中某个质点得到位置，位置和速度的关系、以及旋转四元数和角速度的关系
Figure GDA0003467817310000137
求C(P,Q,X)＝P+Q*R–X关于速度v、角速度ω以及质点位置的梯度值计算作用在刚体和质点上的力。
约束C(P,Q,X)＝||P+Q*R–X||在当前位置P(t)、Q(t)处做线性展开，由位置和速度的关系、以及旋转四元数和角速度的关系可以得到
Figure GDA0003467817310000138
Figure GDA0003467817310000139
Figure GDA0003467817310000141
Figure GDA0003467817310000142
写成向量形式：
Figure GDA0003467817310000143
同理
Figure GDA0003467817310000144
Figure GDA0003467817310000145
Figure GDA0003467817310000146
Figure GDA0003467817310000147
写成向量形式：
Figure GDA0003467817310000148
其中，Pn是当前刚体位移，Qn是当前刚体旋转的四元数，dt是模拟步长，v是刚体线性速度，ω是角速度，对与上面的约束表达式，求它的梯度值计算式子作用在刚体和质点上的力；P为刚体的位置，Q为刚体的旋转四元数。
针(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，推导如下：
GradCi代表约束Ci的梯度ki为该约束的刚性系数；
J为所有GradCi作为行向量所组成的雅可比阵；
Δt代表本模拟时间步长；
M-1代表所有质点的质量倒数组成的对角矩阵；
对于刚体部分来说就是刚体质量的倒数和转动惯量阵的逆矩阵形成的分块对角部分，Fin代表系统由约束产生的内力，Fout代表系统外部施加的力，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置；
Vn为本时间步长模拟前系统新的速度，Xn为本时间步长模拟前系统新的位置；
根据内力的计算公式，系统所有质点或者刚体受到的内力合力为Fin＝∑GradCi*-(Ci*ki)＝∑GradCi*λi＝Jtλ，其中J为所有GradCi组成的雅可比阵，其中
Figure GDA0003467817310000151
λi为每个约束施加在系统所有质点上的力，λ是所有λi组成的列向量、是所有约束产生的内力；
根据牛顿运动公式整个系统的运动速度公式为V＝V0+M-1*(Fin+Fout)*Δt------(1)
位置公式为X＝X0+V*Δt------(2)
将位置的表达式代入Jt中Jt(X0)同时代入C中并且将C在X0处一阶泰勒展开得到C＝C0+J*(X-X0)，Fin＝Jtλ Jtα-1*C＝-Jtα-1(C0+J*(X-X0))
Fin＝Jtλ代入速度表达式式1得到Vn+1＝Vn+M-1*(Fin+Fout)*Δt＝Vn+M-1*(Jtλ+Fout)*Δt-----(3)
Jtα-1*C＝-Jtα-1(C0+J*(X-X0))
由λ＝-α-1C(Xn+1)得到-αλ＝C(Xn+1)
上式在Xn处泰勒展开-αλ＝C(Xn)+J(X-Xn)＝C(Xn)+J*V*Δt-----(4)
将式3代入式4可得到(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2
其中λ为系统内力的大小也就是前面说的-k*C(X)＝-α-1*C(X)
求出λ后可以代入公式3求出本步模拟后的物体速度Vn+1从而求出位置的更新量Xn+1＝Xn+Vn+1Δt。
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，可以用各种已有的线性方程组求解方法。也可以使用传统的Guess迭代方法来行求解。
离散化后四面体可能一部分四面体表示肌肉，一部分表示脂肪，分别有各自的杨氏模量。
本实施例中，对于拥有无限大刚性系数的约束来说Ci得值应该始终为0，而刚性系数ki为无穷大，由于
Figure GDA0003467817310000161
在ki＝∞,Ci＝0时无法确定最终结果，引入新的变量
Figure GDA0003467817310000162
代入Fc的表达式子可以得到：
Fc＝GradCi×λi (5)
Ci＝αi×λi (6)
用新的变量α代替原来的刚性系数，使得可以表示：ki→∞无穷大时候，αi→0，从而根据式6：Ci＝0满足刚性系数无穷大的时候约束始终＝0的要求，同时由于将→∞转换成了→0。所以式5和式6在刚性系数为无穷大的时候都有了意义以及计算机可以表示，同样写成所有合力的矩阵形式：
Figure GDA0003467817310000163
Figure GDA0003467817310000164
其中Matα是由所有αi的对角矩阵。
本实施例中，N是D的相似矩阵D＝L-1NL,先求出D的特征值，再求出特征向量，特征值组成的对角阵就是N,特征向量组成的矩阵就是L。可以采用Eigen软件库都能计算特征值和向量比如Eigen。
Figure GDA0003467817310000171
U等价于包含了6个约束表达式的约束能，N为6维对角矩阵它是矩阵D对角化后得到的一个对角矩阵。
Figure GDA0003467817310000172
为λ的向量表示。D为正定对阵矩阵，实际使用中，并不对D进行分解，采用直接将其作为对角块纳入系统，更加方便推导应变张量的梯度。
最后由于无限大刚性系数的存在，方便的引入刚体的约束，只要定义刚体的约束为C刚体的α＝0即可在现有模型下同一计算。其中刚体的关于约束的梯度和质点的略有不同，它是关于质心位置，和角度旋转的偏导数，将约束在当前时间不做一阶泰勒展开，使得约束变为，线性速度和角速度的表达式，能够更加方便的推导出梯度。
克服了当前半隐式积分，无法实现无穷大刚性系数的问题，以及无法在一个统一的框架下模拟有限元和几何约束能，以及刚体的问题，使得软体和刚体在同一个解算器内完成，大大提高了模拟的视觉效果，同时在保留了有限元的物理真实性的同时，可以很方便的扩展各种几何约束能所产生的力。
若本发明一实施例的系统中有n个不同的约束C，还有m个应变能约束。n个不同的约束C对应n个独立的刚性系数k。
本发明克服了当前半隐式积分，无法实现无穷大刚性系数的问题，以及无法在一个统一的框架下模拟有限元和几何约束能，以及刚体的问题，使得软体和刚体在同一个解算器内完成，大大提高了模拟的视觉效果，同时在保留了有限元的物理真实性的同时，可以很方便的扩展各种几何约束能所产生的力。
本实施例的实时模拟物体运动或形变的方法包括：
计算所有约束的梯度，从而确定系统约束内力的方向；
计算所有外力的和；
将所有约束的梯度组合成雅克比矩阵J；
对于离散化的有限元应变能
Figure GDA0003467817310000181
求出矩阵D的逆矩阵，和刚性系数的逆组成α；
将J、α、质量的倒数带入离散化的物体运动微分方程公式(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解出内力的值；
根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt＝Vn+M-1*(Jtλ+Fout)*Δt，Xn+1＝Xn+Vn+1Δt，更新系数或刚体的位置和速度。
针对离散化的的物体运动微分方程公式(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，系统有n个不同的约束C和m个应变能约束，则
Figure GDA0003467817310000182
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure GDA0003467817310000183
Figure GDA0003467817310000184
系统有m个应变能约束对应有m个D对角块，
Figure GDA0003467817310000185
为第一个应变能约束对应的D的逆，
Figure GDA0003467817310000186
为第m个应变能约束对应的D的逆。
本发明另一实施例的实时模拟物体运动或形变的方法，包括：
确定外力：确定系统所有外力的和Fout；
离散化的应变张量：对于任意一点的应变张量
Figure GDA0003467817310000187
G为形变梯度
Figure GDA0003467817310000188
其中x,y,z为物质空间某点的位置即未形变坐标，X、Y、Z为该点当前的位置坐标，将G离散化表示到四面体的四个顶点，得到G的离散表达式：G＝AB-1，其中矩阵A中的
Figure GDA0003467817310000191
Figure GDA0003467817310000192
为四面体第k个顶点的当前坐标，k＝0,1,2,3，
Figure GDA0003467817310000193
对于矩阵B，其中
Figure GDA0003467817310000194
Figure GDA0003467817310000195
Figure GDA0003467817310000196
为四面体第k个顶点的未形变坐标，k＝0,1,2,3，
Figure GDA0003467817310000197
确定应变能：有限元模拟中的应变能
Figure GDA0003467817310000198
令：
Figure GDA0003467817310000199
对Eu求梯度，所形成的内力为
Figure GDA00034678173100001910
V为四面体体积，VD*Cu为内力的大小，
Figure GDA00034678173100001911
为内力方向，ε＝(ε11 ε22 ε33 ε12 ε13 ε23)为6维的格林应变张量，D为杨氏模量和泊松系数组成的弹性模量阵，VD刚性系统矩阵为
Figure GDA00034678173100001912
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure GDA00034678173100001913
ki为约束Ci的刚性系数，系统有m个应变能约束对应有m个D对角块，
Figure GDA00034678173100001914
为第一个应变能约束对应的D的逆，
Figure GDA00034678173100001915
为第m个应变能约束对应的D的逆；D为正定对阵矩阵，实际使用中，并不对D进行分解，采用直接将其作为对角块纳入系统，更加方便推导应变张量的梯度；
确定应变能约束的梯度：应变能约束的梯度向量为雅克比矩阵
Figure GDA00034678173100001916
求解内力：根据离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解内力大小λ＝VD*Cu，其中，v为系统当前所有质点速度组成的列向量，M-1代表所有质点的质量倒数组成的对角矩阵，Δt表示本模拟时间步长，α＝VD为刚性系统矩阵,C(Xn)为tn时刻系统的质点位置的3n维向量，或刚体速度和角速度组成的6n维向量；
更新位置、速度：根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt求出本时间步模拟后的物体速度Vn+1，根据Xn+1＝Xn+Vn+1Δt，其中，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置，Vn为本时间步长模拟前系统的速度，Xn本时间步长模拟前系统新的速度和位置，Fin代表系统由约束产生的内力，根据求解出的速度、位置更新系统的质点或刚体的位置、速度。
进一步，
Figure GDA0003467817310000201
为六维应变张量的六个内力方向组成的矩阵。
进一步，本实施例的应变张量
Figure GDA0003467817310000202
G为四面体的形变梯度的离散化表示G＝AB-1＝(A1A2A3)(c1c2c3)，B由四面体的顶点在未形变的时候的位置组成的矩阵，其逆矩阵B-1三个列向量为c1 c2 c3，由于四面体的体积总是>0，故B的行列式值等于B的混合积等于四面体体积，始终不等于0，故B-1始终存在，G＝(g1，g2，g3)，所以
Figure GDA0003467817310000203
其中δij＝1(i＝j)否则＝0；
根据G的表达式可知gi＝P*ci，gj＝P*cj
所以
Figure GDA0003467817310000204
Figure GDA0003467817310000205
同理
Figure GDA0003467817310000211
Figure GDA0003467817310000212
Figure GDA0003467817310000213
Figure GDA0003467817310000214
其中，ci为B-1的第i个列向量，cik为B-1的第i行第k列的元素。
D本身是一个6×6的对称矩阵，其中元素由杨氏模量和泊松系数组成。将它的逆放入刚性矩阵α对应的对角块的位置。对于非有限元的约束，刚性系数简单的是个标量放入对角线元素对应的位置。对具有有限元应变能约束的刚性矩阵α，将D的逆矩阵放在对角块和刚性系数的逆组成。
本实施例的矩阵
Figure GDA0003467817310000215
Figure GDA0003467817310000216
对于离散化的有限元应变能，求出矩阵D的逆矩阵D-1，根据矩阵D的逆矩阵D-1和刚性系数的倒数组成α。
本发明从约束能的力学模型开始逐步展开，依次实现，无穷大刚性系数的实现；传统有限元方法统一到本模型；刚体模拟统一到本模型中。
本发明提出的半隐式积分模拟方法，可以实现几何约束能和有限元应变能的统一，同时引入的正则化方法可以使得物体的刚性系数到无穷大。从而兼容刚体模拟和软体模拟。实得模拟在一个统一的框架下，达到了很棒的视觉效果和稳定性，同时并没有增加任何的运算开销。
本发明保留了软体刚体模拟中半隐式积分的稳定性，同时在引入新的方法后实现了无穷大刚性系数的算法同时实现了有限元模拟和几何能力学实时模拟算法的统一性。在保留有限元方法的高真实性的同时极大的增加了自由度和可扩展性，同时可以在一个解算器内模拟刚体和软体，使得视觉效果大大增强。
以上述依据本申请的理想实施例为启示，通过上述的说明内容，相关工作人员完全可以在不偏离本项申请技术思想的范围内，进行多样的变更以及修改。本项申请的技术性范围并不局限于说明书上的内容，必须要根据权利要求范围来确定其技术性范围。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和/或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和/或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能。

Claims (10)

1.一种实时模拟物体运动或形变的方法，其特征在于，包括：
确定约束内力方向：每个约束对系统产生的约束能为
Figure FDA0003467817300000011
根据约束能计算约束能的梯度为该能量产生的内力Fc＝-GradCi×kiCi,其中，ki为约束Ci的刚性系数，k为约束C的刚性系数，GradCi为约束Ci关于整个系统坐标的梯度，所有约束势能之和产生的内力合力是所有产生内力的约束能内力之和
Figure FDA0003467817300000012
变换成矩阵形式为
Figure FDA0003467817300000013
中J是为所有GradCi作为行向量所组成的雅可比阵，K是由每个约束的刚性系数组成的对角矩阵，
Figure FDA0003467817300000014
是系统所有约束组成的列向量，-Ci×ki为约束Ci的内力大小，GradCi为约束Ci的内力方向，J的每一列为一个内力的方向，
Figure FDA0003467817300000015
为内力大小；
确定所有约束能的梯度：将所有约束能的梯度组合成雅可比矩阵J(C1,C2....Ck)，其中C1、C2……Ck代表系统中的k个约束；
引入新变量：令
Figure FDA0003467817300000016
得到单个约束i作用在系统所有质点的内力Fc＝GradCi×λi，Ci＝αi×λi，系统所有质点或刚体受到的内力合力Fin＝∑GradCi×(-(Ci×ki))＝∑GradCi×λi＝Jtλ，
Figure FDA0003467817300000017
其中，GradCi代表约束Ci的梯度，ki为约束Ci的刚性系数，J为所有GradCi作为行向量所组成的雅可比阵，Jt为J的转置矩阵，K为由每个约束的刚性系数组成的对角矩阵，
Figure FDA0003467817300000019
为系统所有约束组成的列向量，λ为所有约束产生的内力大小、为所有λi组成的列向量，Matα由所有αi的对角矩阵，是个对角阵或者分块对角阵，它的对角元素为每个约束刚性系数ki的倒数αi＝1/ki，其他元素都为0的矩阵
Figure FDA0003467817300000018
离散化后的约束能：格林应变张量作为应变系数，应力张量E＝D×ε，其中D为由杨氏模量和泊松系数组成的弹性模量阵、为正定对阵矩阵，物体处一点的应变能εt x E＝εt×D×ε，离散化到体积四面体后，得到一个体积四面体产生的约束能是
Figure FDA0003467817300000021
其中V为四面体的体积，ε为格林应变张量，将D分解为L-1NL，得出
Figure FDA0003467817300000022
其中，N为对角矩阵，N是D的相似矩阵，N为D的特征值组成的对角矩阵，L为D的特征向量组成的矩阵；
确定外力：确定系统所有外力的和Fout；
求解内力：离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解内力大小
Figure FDA0003467817300000025
其中，v为系统当前所有质点速度组成的列向量，M-1代表所有质点的质量倒数组成的对角矩阵，Δt表示本模拟时间步长，C(Xn)为tn时刻系统的质点位置的3n维向量，或刚体速度和角速度组成的6n维向量，α为刚性系数矩阵；
更新位置、速度：根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt＝Vn+M-1*(Jtλ+Fout)*Δt求出本时间步长模拟后的物体速度Vn+1，根据Xn+1＝Xn+Vn+1Δt求出本时间步长模拟完成后系统新的速度，根据求解出的速度、位置更新系统的质点或刚体的位置、速度，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置，Vn为本时间步长模拟前系统的速度，Xn本时间步长模拟前系统新的速度和位置，Fin代表系统由约束产生的内力。
2.根据权利要求1所述的实时模拟物体运动或形变的方法，其特征在于，所述离散化后的约束能步骤中U等价于包含了6个约束表达式的约束能，N为6维对角矩阵、它是矩阵D对角化后得到的一个对角矩阵；
所述矩阵
Figure FDA0003467817300000023
Figure FDA0003467817300000024
其中E为杨氏模量，μ为泊松系数。
3.根据权利要求1所述的实时模拟物体运动或形变的方法，其特征在于，若为软体则C为
Figure FDA0003467817300000036
所有约束能的梯度组合成雅可比矩阵J：
Figure FDA0003467817300000031
其中C1、C2……Ck代表系统中的所有的k个约束，xn为系统第n个质点的坐标、为一个三维向量。
4.根据权利要求1所述的实时模拟物体运动或形变的方法，其特征在于，刚体的约束C为C(P，Q)其中P为刚体的位置、为3维向量，Q为刚体的旋转四元数、为4维向量，刚体位置和速度的关系、以及刚体旋转四元数和角速度的关系如下：
Figure FDA0003467817300000032
对C(P，Q)求关于速度v、角速度ω的梯度值计算出作用在刚体上的力。
5.根据权利要求1至4任意一项所述的实时模拟物体运动或形变的方法，其特征在于，刚体上的点和系统中软体上的一个质点的固联约束为C(P，Q，X)＝P+Q*R-X，其中P为刚体的位置，Q为刚体的旋转四元数，R为固联点在刚体本地坐标系下的向量，X为软体中某个质点得到位置，位置和速度的关系、以及旋转四元数和角速度的关系
Figure FDA0003467817300000033
Figure FDA0003467817300000034
求C(P，Q，X)＝P+Q*R-X关于速度v、角速度ω、以及质点位置的梯度值计算作用在刚体和质点上的力。
6.根据权利要求1至4任意一项所述的实时模拟物体运动或形变的方法，其特征在于，刚体上的点和系统中软体上的一个质点的固联约束C(P，Q，X)＝||P+Q*R-X||在当前位置做线性展开，由位置和速度的关系、以及旋转四元数和角速度的关系得到
Figure FDA0003467817300000035
Figure FDA0003467817300000041
Figure FDA0003467817300000042
Figure FDA0003467817300000043
写成向量形式：
Figure FDA0003467817300000044
同理
Figure FDA0003467817300000045
Figure FDA0003467817300000046
Figure FDA0003467817300000047
Figure FDA0003467817300000048
写成向量形式：
Figure FDA0003467817300000049
其中，Pn是当前刚体位移，Qn是当前刚体旋转四元数，dt是模拟步长，v是刚体线性速度，ω是角速度，求它关于速度v、角速度ω以及质点位置的梯度值计算作用在刚体和质点上的力，P为刚体的位置，Q为刚体的旋转四元数，R为固联点在刚体本地坐标系下的向量，X为软体中某个质点得到位置。
7.根据权利要求1至4任意一项所述的实时模拟物体运动或形变的方法，其特征在于，所述求解内力步骤前，还包括：确定刚性系数矩阵α，若系统有n个不同的约束C，则
Figure FDA0003467817300000051
若系统有n个不同的约束C和m个应变能约束，则
Figure FDA0003467817300000052
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure FDA0003467817300000053
Figure FDA0003467817300000054
ki为约束Ci的刚性系数，系统有m个应变能约束对应有m个D对角块，
Figure FDA0003467817300000055
为第一个应变能约束对应的D的逆，
Figure FDA0003467817300000056
为第m个应变能约束对应的D的逆。
8.一种实时模拟物体运动或形变的方法，其特征在于，包括：
确定外力：确定系统所有外力的和Fout；
离散化的应变张量：对于任意一点的应变张量
Figure FDA0003467817300000057
G为形变梯度
Figure FDA0003467817300000058
其中x，y，z为物质空间某点的位置即未形变坐标，X、Y、Z为该点当前的位置坐标，将G离散化表示到四面体的四个顶点，得到G的离散表达式：G＝AB-1，其中矩阵A中的
Figure FDA0003467817300000059
Figure FDA00034678173000000510
为四面体第k个顶点的当前坐标，k＝0，1，2，3，
Figure FDA00034678173000000511
对于矩阵B，其中
Figure FDA00034678173000000512
Figure FDA00034678173000000513
Figure FDA00034678173000000514
为四面体第k个顶点的未形变坐标，k＝0，1，2，3，
Figure FDA00034678173000000515
Figure FDA00034678173000000516
确定应变能：有限元模拟中的应变能
Figure FDA0003467817300000061
令
Figure FDA0003467817300000062
对Eu求梯度，所形成的内力为
Figure FDA0003467817300000063
VD*Cu为内力大小，V为四面体体积，ε＝(ε11 ε22 ε33 ε12ε13 ε23)为6维的格林应变张量，D为杨氏模量和泊松系数组成的弹性模量阵，VD刚性系统矩阵
Figure FDA0003467817300000064
系统n个不同的约束C就对应有n个独立的刚性系数k，
Figure FDA0003467817300000065
ki为约束Ci的刚性系数，系统有m个应变能约束对应有m个D对角块，
Figure FDA0003467817300000066
为第一个应变能约束对应的D的逆，
Figure FDA0003467817300000067
为第m个应变能约束对应的D的逆；
确定应变能约束的梯度：应变能约束的梯度向量为雅克比矩阵
Figure FDA0003467817300000068
求解内力：根据离散化后的物体运动微分方程求解公式：
(JM-1Jt+Δt-2*α)λ＝-Δt-1Jv-JM-1Fout-C(Xn)Δt-2，求解内力大小λ＝VD*Cu，其中，v为系统当前所有质点速度组成的列向量，M-1代表所有质点的质量倒数组成的对角矩阵，Δt表示本模拟时间步长，α＝VD为刚性系统矩阵，C(Xn)为tn时刻系统的质点位置的3n维向量、或刚体速度和角速度组成的6n维向量；
更新位置、速度：根据Vn+1＝Vn+M-1*(Fin+Fout)*Δt求出本时间步模拟后的物体速度Vn+1，根据Xn+1＝Xn+Vn+1Δt，其中，Vn+1为本时间步长模拟完成后系统新的速度，Xn+1为本时间步长模拟完成后系统新的位置，Vn为本时间步长模拟前系统的速度，Xn本时间步长模拟前系统新的速度和位置，Fin代表系统由约束产生的内力，根据求解出的速度、位置更新系统的质点或刚体的位置、速度。
9.根据权利要求8所述的实时模拟物体运动或形变的方法，其特征在于，所述应变张量
Figure FDA0003467817300000071
G为四面体的形变梯度的离散化表示G＝AB-1＝(A1A2A3)(c1c2c3)，B由四面体的顶点在未形变的时候的位置组成的矩阵，其逆矩阵B-1三个列向量为c1 c2 c3，由于四面体的体积总是＞0，故B的行列式值等于B的混合积等于四面体体积，始终不等于0，故B-1始终存在，G＝(g1，g2，g3)所以
Figure FDA0003467817300000072
其中δij＝1(i＝j)否则＝0；
根据G的表达式可知gi＝P*ci，gj＝P*cj
所以
Figure FDA0003467817300000073
Figure FDA0003467817300000074
同理
Figure FDA0003467817300000075
Figure FDA0003467817300000076
Figure FDA0003467817300000077
Figure FDA0003467817300000078
Figure FDA0003467817300000079
其中，ci为B-1的第i个列向量，cik为B-1的第i行第k列的元素，gi为G的第i个列向量。
10.根据权利要求8或9所述的实时模拟物体运动或形变的方法，其特征在于，所述
Figure FDA00034678173000000710
为六维应变张量的六个内力方向组成的矩阵，
所述矩阵
Figure FDA00034678173000000711
Figure FDA0003467817300000081
其中E为杨氏模量，μ为泊松系数。
CN202110153901.4A 2021-02-04 2021-02-04 实时模拟物体运动或形变的方法 Active CN112883609B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110153901.4A CN112883609B (zh) 2021-02-04 2021-02-04 实时模拟物体运动或形变的方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110153901.4A CN112883609B (zh) 2021-02-04 2021-02-04 实时模拟物体运动或形变的方法

Publications (2)

Publication Number Publication Date
CN112883609A CN112883609A (zh) 2021-06-01
CN112883609B true CN112883609B (zh) 2022-02-18

Family

ID=76057213

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110153901.4A Active CN112883609B (zh) 2021-02-04 2021-02-04 实时模拟物体运动或形变的方法

Country Status (1)

Country Link
CN (1) CN112883609B (zh)

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108986220A (zh) * 2018-07-16 2018-12-11 浙江大学 一种加速有限元求解物体网格模型弹性变形的方法
KR101948482B1 (ko) * 2018-08-31 2019-02-14 숭실대학교산학협력단 가변형 연체 조직의 물리 시뮬레이션 방법, 이를 수행하기 위한 기록매체 및 장치
CN110992456A (zh) * 2019-11-19 2020-04-10 浙江大学 一种基于位置动力学的雪崩模拟方法
CN113409443A (zh) * 2021-05-19 2021-09-17 南昌大学 一种基于位置约束和非线性弹簧的软组织建模方法

Family Cites Families (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US11200356B2 (en) * 2019-04-09 2021-12-14 Nvidia Corporation Using a computer to model the reactions of objects to simulated physical interactions

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN108986220A (zh) * 2018-07-16 2018-12-11 浙江大学 一种加速有限元求解物体网格模型弹性变形的方法
KR101948482B1 (ko) * 2018-08-31 2019-02-14 숭실대학교산학협력단 가변형 연체 조직의 물리 시뮬레이션 방법, 이를 수행하기 위한 기록매체 및 장치
CN110992456A (zh) * 2019-11-19 2020-04-10 浙江大学 一种基于位置动力学的雪崩模拟方法
CN113409443A (zh) * 2021-05-19 2021-09-17 南昌大学 一种基于位置约束和非线性弹簧的软组织建模方法

Also Published As

Publication number Publication date
CN112883609A (zh) 2021-06-01

Similar Documents

Publication Publication Date Title
Sifakis et al. FEM simulation of 3D deformable solids: a practitioner's guide to theory, discretization and model reduction
CN107330972B (zh) 模拟生物力学特性的实时软组织形变方法和系统
Roberson et al. Dynamics of multibody systems
Terzopoulos et al. Physically based models with rigid and deformable components
Metaxas et al. Dynamic deformation of solid primitives with constraints
US8190412B2 (en) Method of simulating deformable object using geometrically motivated model
Essahbi et al. Soft material modeling for robotic manipulation
Basdogan Real-time simulation of dynamically deformable finite element models using modal analysis and spectral lanczos decomposition methods
EP3296068A1 (en) Robot behavior generation method
Hess-Coelho et al. Modular modelling methodology applied to the dynamic analysis of parallel mechanisms
CN103426196B (zh) 一种流体环境下的关节动画建模方法
Aguirre et al. An implicit 3D corotational formulation for frictional contact dynamics of beams against rigid surfaces using discrete signed distance fields
Lee et al. Volumetric Object Modeling Using Internal Shape Preserving Constraint in Unity 3D.
Nesme et al. Physically realistic interactive simulation for biological soft tissues
Scaglioni et al. Closed-form control oriented model of highly flexible manipulators
CN112883609B (zh) 实时模拟物体运动或形变的方法
Lennartsson Efficient multibody dynamics
CN117131742A (zh) 基于可逆转有限元能量法的人体组织模拟在gpu下的高并发实现方法
Li et al. Soft articulated characters in projective dynamics
Huang et al. A survey on fast simulation of elastic objects
CN115455753A (zh) 一种软组织与刚性地面的碰撞仿真方法及装置
Han et al. A Convex Formulation of Frictional Contact between Rigid and Deformable Bodies
Lugrís et al. Efficient calculation of the inertia terms in floating frame of reference formulations for flexible multibody dynamics
Sase et al. Haptic rendering of contact between rigid and deformable objects based on penalty method with implicit time integration
Taves et al. Dwelling on the Connection Between SO (3) and Rotation Matrices in Rigid Multibody Dynamics–Part 2: Comparison Against Formulations Using Euler Parameters or Euler Angles

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant