FEniCS

7 minute read

既然之前保证过的,那么就现在有点时间介绍下我现在用的有限元计算软件 FEniCS。由于设计数值计算,所以把这篇博文归类到数学里面去了。

首先要注意,我说的是广义上的有限元计算,并不仅限于我的本行力学。这个软件所做的,就是用有限元方法数值求解一个偏微分方程(组)。也就是说,你需要做的,就是得告诉程序

  1. 你要求解的方程是什么。要注意,必须采用 Formulation faible/variationnelle 而不能采用强形式。这个大家做有限元的都知道。
  2. 定义几何,更具体的说就是定义一个 Maillage。
  3. 你想用什么有限元求解方程,也就是说定义 Espace de discrétisation。

至于其他琐碎的事情比如 Assemblage 之类的 FEniCS 会自动帮你做,这就是它强大的地方。此外你可操控的自由度也很大,你可以自己选择 Schéma d’intégration,选择线性方程组的求解方法之类的。

下面简单举个例子,就拿最最常见的 Poisson 下手吧。

其中我们假设外力 等于

此外边界条件为最简单的 上。

第一步就是把强形式变成弱形式了,得到寻找 ,使得对于所有 ,有

然后就可以打开你们喜欢的编辑器开始 Python 啦!第一步从加载 FEniCS 包开始

from dolfin import *

然后我们定义几何空间并网格化,FEniCS 里面有一个函数可以定义二维单位平方,所以我们直接使用就好了。在一般情况下你可能需要自己在 FEniCS 里面定义几何然后用 CGAL 去 Mailler,或者一开始从几何开始就用其他的 Mailleur,比如 GMSH。

n = 100
mesh = UnitSquareMesh(n, n)

这里我们把 用三角形单元进行划分,每边有100个三角形的边。如果需要 Raffiner 那么只需要变大 就好了。

然后我们定义有限元函数空间。还记得当初在巴黎 CDD 的时候来过一个小教授讲了有限差分法和有限元法的区别,他讲得很好,就是说 DF 去近似的是导数,把各种导数用有限差分去代替;而有限元则是直接去近似函数,即把无线维函数空间用一个有限维函数空间去代替。

V = FunctionSpace(mesh, 'CG', 1)
u, v = TrialFunction(V), TestFunction(V)
U = Function(V)

这里 V 就是我们的有限元空间,CG 代表我们采用最普遍的 Lagrange 有限元,阶数是1,即每个三角形单元上我们用线性函数去近似解。其次的 uv 定义了弱形式中的解和测试函数,而 U 定义了一个 上的函数,之后用来放置解。

现在开始定义外力 ,在 FEniCS 中很简单,只需要做

f = Expression('sin(pi*x[0])*sin(pi*x[1])')

即可,在 Expression 中使用 C++ 语法定义一个函数,x[0] 代表第1个坐标 x[1] 代表第2个坐标

下面就要定义要求解的问题了,具体来说就要定义弱形式中的双线性形式 和线性形式

a = inner(grad(u), grad(v))*dx
L = f*v*dx

是不是很简单,连解释都免了。之后就差定义边界条件了,我们有

def bord(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, Constant(0.0), bord)

其中 bord 定义了边界,on_boundary 是 FEniCS 提供的,意思就是在网格的边界上,由于我们这里的边界条件就是完全定义在边界上的,所以直接用就好了(否则需要自己定义一个在边界上对 x 的条件)。然后的 DirichletBC 函数定义一个 Dirichelet 的边界条件,也很直观吧。

最后就是求解了,我们有

solve(a == L, U, bc)

并把解存到我们之前定义的 U 中。之后如果想简单地看求解结果的话,只需要

plot(U, interactive=True)

即可,是不是很简单!

Tags:

Updated:

Comments