GSoC 2017:FEniCS 有限元工具箱的包装器

2017 年 9 月 1 日 | Yiannis Simillides、Bart Janssens、Chris Rackauckas

在这个 Google Summer of Code 项目中,我和我的导师一起努力创建了一个 FEniCS 有限元工具箱在 Julia 语言中的包装器。我们的工作成果可以在 FEniCS.jl 中找到。这将允许用户在 Julia 中直接执行 FEM 计算,利用我们的 PyCall.jl 包装功能。我们目前已经包装了主要功能,并提供了在必要时添加更多组件的必要说明。Julia 社区中与该项目没有直接关系的成员也在整个项目中贡献了一些小的修复和建议。大部分生成的代码已经合并到 GitHub 仓库(专门为该项目创建)。其中一项主要改进将极大地提高其使用率,那就是与 JuliaDiffEq 包的进一步集成。

  1. 什么是 FEniCS?
  2. 但是,它是如何工作的?
  3. 那么,内部呢?
  4. 演示
  5. 挑战
  6. 进一步改进
  7. 站在巨人的肩膀上
  8. 致谢

什么是 FEniCS?

FEniCS 自称是一个流行的开源 (LGPLv3) 计算平台,用于求解偏微分方程 (PDE)。FEniCS 使用户能够快速将科学模型转换为高效的有限元代码,即它允许我们描述(某些)复杂的数学方程并使用计算机模拟自动求解它们。

但是,它是如何工作的?

有限元方法是一种数值方法,它通过求解函数在某些稀疏基上的 Galerkin 近似的弱形式来求解偏微分方程。这是通过离散化域(将其分解成小的部分,通常是称为有限元的三角形)来完成的,通过三角形节点处的数值表示函数。这将导致一个方程组,它是在这些三角形的基底上对初始偏微分方程的解释。这几乎不可能手工完成,因此需要计算机软件来帮助求解它。我们的包装器提供了对网格划分功能、刚度矩阵的组装和变分问题的求解的调用。它还提供了一些辅助函数的访问权限,这些函数可以简化使用。下面可以演示一些示例网格

dolphin mesh circle mesh

那么,内部呢?

我们目前提供宏,用于从 FEniCS 类创建 Julia 类型,以帮助函数重载。除此之外,我们可以包装属性,并将这两者与 PyCall.jl 功能结合起来,我们可以直接从 FEniCS 包装必要的功能。我们定义了执行一些线性代数运算的方法,并定义了用于各种几何对象和 UFL 形式的运算符函数+,-,*。我们导出的 API 基本上与 Pythonic FEniCS 相同,只是略有变化。

演示

下面是一个小型演示,说明用户将如何使用我们的代码来求解具有狄利克雷条件的泊松方程。这直接反映了 教程 中 FEniCS 提供的一个教程

(u)=f-\bigtriangleup(u) = f 在单位正方形中

u=uDu = u_D 在边界上

uDu_D == 1+x2+2y21 + x^2 + 2y^2

f=6f = -6
using FEniCS
mesh = UnitSquareMesh(8,8)
V = FunctionSpace(mesh,"P",1)
u_D = Expression("1+x[0]*x[0]+2*x[1]*x[1]", degree=2)
u = TrialFunction(V)
bc1 = DirichletBC(V,u_D, "on_boundary")
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u),grad(v))*dx
L = f*v*dx
U = FEniCS.Function(V)
lvsolve(a,L,U,bc1) #linear variational solver
errornorm(u_D, U, norm="L2")
get_array(L) #this returns an array for the stiffness matrix
get_array(U) #this returns an array for the solution values
vtkfile = File("poisson/solution.pvd")
vtkfile << U.pyobject #exports the solution to a vtkfile

除了定义问题之外,我们还可以访问并保存与各种变分形式相对应的数组。这些将返回一个 数组类型 我们可以这样做:

a = dot(grad(u),grad(v))*dx #this sets up the variatonal form from the previous problem
variational_matrix = get_array(a)

我们还可以绘制解决方案(这依赖于 FEniCS 后端进行绘制)

FEniCS.Plot(mesh)
FEniCS.Plot(U)

Square Mesh

Solution

挑战

由于该项目依赖于 FEniCS,因此我们在整个夏季都面临着各种挑战。这些包括但不限于构建错误,其中该包的各个部分无法编译,以及代码可用性方面的意外错误。Chris 和 Bart 在指出这些问题以及帮助解决这些问题方面非常有帮助。在某些部分,文档略显零散,这也使得项目的某些部分变得复杂,因为一些函数对它们的预期用途模棱两可。

进一步改进

我希望能够维护和改进该包,并在我进一步的学习中尽可能地使用它。按难度顺序排列,一些可识别的改进包括

站在巨人的肩膀上

除了编码(这非常愉快并提供了独特的学习体验之外),参加这个暑期项目让我接触到一个很棒的社区。在与 Julia 短暂合作的时间里,我有机会访问了伦敦的 JuliaHub 办公室。之后,我获得了JuliaHubNumFOCUS 的资金,以参加JuliaCon2017 并展示海报。除了精彩的演讲外,我还借此机会与其他 GSoC 学生合租公寓,并与 Julia 社区的杰出成员一起享用午餐和饮品。我真正相信这是开源社区的奇妙之处之一。人们投入时间和精力来帮助其他人,并传播对所有人开放的科学发现。

JuliaCon

JuliaCon 2017 参会者

致谢

首先,我要感谢我的导师 Chris 和 Bart。Chris 尽管时差很大,但也一直都在帮助我回答(很多)问题并提供建议。Bart 找到了代码中的许多初始错误和不一致之处,提供了修复这些错误的必要信息。JuliaHub 以及 NumFOCUS 为我提供了参加 JuliaCon 2017 并展示海报的资金。最后感谢 Google 开源计划,他们为我提供了必要的资金,使我能够在整个夏季完成这个项目,并获得了一次美妙的体验。