Homework 1: 3D FEM (Explicit, Jacobi and conjugate gradient method)

实现了基于 Neohookean 弹性模型的 3D 线性有限元,并对比了三种数值方法:1.显式时间积分;隐式时间积分(2. Jacobi方法求解,3.共轭梯度法求解)

一周前就写完了,这几天忙着入职和租房,一直没有时间写分析……今天抽时间写了出来,顺便渲染了结果便于对比。

代码及分析-github

Introduction

An implementation of 3D linear FEM based on Neohookean elasticity model (No damping force).

Simulation Results

out-exp out-ja out-cg

Using Explicit, Jacobi and CG method respectively (Rendered by Blender).

Analysis

For Explicit method, there is no constraints between nodes when updating, so the system explodes if time step or Young’s modulus is large.

For Jacobi method, the convergence is determined by the spectral radius \rho of the matrix -D^{-1}(L+U) , Jacobi method converges if and only if \rho is less than 1.

For CG method, the system always converges since the coefficient matrix A is positive definite for (solid) FEM problems.

Same Scene

0-exp 0-ja 0-cg

Explicit and Jacobi method have to update many times using a small time step at each frame.

\rho = 0.0036

Large Time Step

1-exp 1-ja 1-cg

\rho = 1.09

Large Young’s Modulus

2-exp 2-ja 2-cg

\rho = 1.13

Problems of Conjugate Gradient Method

3-cg

  1. Object rotates slowly when time step is large.

This problem is caused by the lack of damping force, add damping force to object can avoid this problem.

4-cg

  1. CG doesn’t converge when both time step and Young’s modulus are large.

I still don’t know why CG doesn’t converge in this condition, I guess this is caused by numerical accuracy.

16 个赞

超级棒。
问一下,不同方法的性能有多大差距?

1 个赞

现在点数比较少,所以性能区别不太大,都是可以实时的。因为我的笔记本很差,每帧绘制得很慢,所以fps有点低……

2 个赞

这个工作对我很有启发。我对有限元的认识不是很深刻,想请教个问题:
对于隐式欧拉方法,需要计算dP/dF,如何计算 dP/dF,有相关书籍推荐嘛?

hi, @tianlajiangjun ,dP/dF的推导在我提到的参考文献 The classical FEM method and discretization methodology 这篇文章中有介绍。

1 个赞

好的,感谢你的帮助

我也遇到相同问题。当杨氏模量与时间步长都比较大时,CG方法是会模拟失败:frowning:,只有把时间步长调小才可以稳定运行,这点应该如何优化呢?

这个需要直接看一下求解 Ax=b 的时候,A 在调整参数的时候性质如何变化,主对角占优越明显,收敛越快,越不明显收敛越慢,最后甚至被数值误差淹没导致发散。想要改善这种情况,多可以使用预处理(pre-conditioner)来改善矩阵性质。关于性质可以参考:

1 个赞

好的,我好好看看,感谢你的帮助

这里只是介绍了条件数的基本原理和计算方法,没有提到条件数如何优化。我看有些说可以加阻尼来降低条件数,还有就是可以提高舍入误差精度也就是提高运算位数来延后发散时间,除此之外还有别的什么办法吗?