Gridap 简介

Gridap (Github, 官方文档) 是一个基于 Julia 编程语言 的求解偏微分方程 (PDE) 的网格剖分逼近方法的函数库, 涵盖了线性/非线性、标量/向量型 PDE 的协调/非协调元等众多类型的有限元逼近方法. 可使用 Julia 的包管理工具 Pkg 直接安装.

  • 安装方法一 (REPL):

    julia> using Pkg
julia> Pkg.add("Gridap")

  • 安装方法二 (在 Julia REPL 中输入 ] 进入包管理模式)

    pkg> add Gridap

Poisson 方程的弱形式与有限元格式

考虑如下齐次 Dirichlet 边界条件的 Poisson 方程:

$$ \begin{equation}\label{eq} -\Delta u = f, \quad\text{ in }\; \Omega,\qquad u = 0,\quad \text{ on }\;\partial\Omega, \end{equation} $$

其中 $\Omega=(0,1)\times (0,1)$ 是一个求解区域, $u\in H^2(\Omega)\cap H_0^1(\Omega)$ 是待求函数, $f\in L^2(\Omega)$ 是一个给定的函数.
方程两边作用于实验函数 $v\in H_0^1(\Omega)$, 就得到了该方程的弱形式: 求解 $u\in H_0^1(\Omega)$, 满足

$$ \begin{equation} a(u,v) = b(v),\quad \forall v\in H_0^1(\Omega), \end{equation} $$

这里, 若使用 $(\cdot,\cdot)$ 表示 $L^2(\Omega)$ 内积, 那么双线性型 $a: H_0^1(\Omega)\times H_0^1(\Omega)\longrightarrow \mathbb{R}$ 为

$$ \begin{equation} a(u,v)=(\nabla u,\nabla v) = \int_{\Omega}\nabla u\cdot\nabla v, \end{equation} $$

线性泛函 $b:H_0^1(\Omega)\longrightarrow\mathbb{R}$ 即为内积 $(f, v)$.
设 $V_h\subset H_0^1(\Omega)$ 为有限元空间, 则 Poisson 方程 $\eqref{eq}$ 的 (协调) 有限元离散格式为: 求解 $u_h\in V_h$, 满足

$$ \begin{equation} a(u_h,v_h) = b(v_h), \quad \forall v_h\in V_h. \end{equation} $$

网格剖分

这一部分介绍网格剖分的两种方法.

Gmsh API

Gmsh 是一款著名的免费开源的有限元网格生成器, 不仅提供了适用于 Windows/macOS/Linux 等系统的图形化操作界面, 还提供了这些操作系统的 C, C++, Fortran, Python 以及 Julia 等编程语言的 SDK 以进行 API 接入. 下面就通过 Gmsh 的 Julia API 来生成网格.
首先, 下载 SDK 使用 API. 进入下载页面选择相对应操作系统的 SDK 压缩包进行下载, 然后解压到你认为合适的位置, 这里假设安装路径为 ~/gmsh/. 然后设置环境变量(Julia 的运行加载路径) JULIA_LOAD_PATH 的位置为 SDK 下 lib 文件夹, 在这里即为 ~/gmsh/lib. 此时, 即可在任意位置的 Julia 程序中使用 import 来引入 gmsh 模块

import gmsh

在开始之前, 初始化以打开 API 接口

gmsh.initialize()

紧接着, 使用 gmsh.model.add 新建空白模型, 参数为模型名称

gmsh.model.add("square")

得到了这个空白模型之后, 就要对这个模型添加一系列的基本几何特征来一步一步地构建出所研究的具体几何图形. 第一种要添加的几何特征就是"点"了, 所使用的函数为 gmsh.model.geo 下的 addPoint. 该函数有 3 个必选参数和 2 个可选参数:

  • 必选参数: 即点的 $x, y, z$ 坐标, 若考虑的问题是 2 维的, 把 $z$ 轴坐标设为 $0$ 即可.
  • 第一个可选参数为该点附近网格的尺寸, 若无, 将基于模型整体大小给出一个适合的值.

  • 第二个可选参数为该点的标号/签(正整数, 不能重复), 如果没有显式给出的话, 函数会自动生成一个标签并且作为函数的返回值输出. 我们可以利用这个返回值自定义一个引用, 在程序中个性化点的名称 (下面的 $p_i$ 即为引用).

    lc = 0.05
p₁ = gmsh.model.geo.addPoint(0,0,0,lc)
p₂ = gmsh.model.geo.addPoint(1,0,0,lc)
p₃ = gmsh.model.geo.addPoint(1,1,0,lc)
p₄ = gmsh.model.geo.addPoint(0,1,0,lc)

    第二种要给定的几何特征就是点所构成的"线"了, 对于我们这个模型来说, 是最简单的两个点所确定的直线, 对应的函数为 addLine, 使用方法与上面的 addPoint 类似, 前两个必选参数为起点和终点的标号, 可选参数为这条线的标签(当然, 也不能重复). 需要注意的, 这些线是有方向的, 这个方向用来确定形成面的闭合曲线.

    l₁ = gmsh.model.geo.addLine(p₁, p₂)
l₂ = gmsh.model.geo.addLine(p₂, p₃)
l₃ = gmsh.model.geo.addLine(p₃, p₄)
l₄ = gmsh.model.geo.addLine(p₁, p₄)

    形成闭合曲线的函数是 addCurveLoop, 需传入两个参数, 第一个参数为按照一定顺序排列的线的标号所构成的向量, 第二个(可选)参数是标签, 在确定顺序的时候, 可以在线的标号前面加上负号来表示本来方向的反方向. 最后就是用 addPlaneSurface 添加"面", 传入的参数与闭合曲线类似, 即形成这个面的闭合曲线的标号所组成的向量和这歌面的自定义标号.

    curve = gmsh.model.geo.addCurveLoop([l₁, l₂, l₃, -l₄])
surface = gmsh.model.geo.addPlaneSurface([curve])

    至此, 这个正方形的几何特征就全部确定了, 接下来就该生成网格了. 不过, 在剖分之前, 还需要进行两个操作, 分别叫做同步和分组, 其中同步是为了生成模型所对应的数据结构, 对应的函数为 synchronize

    gmsh.model.geo.synchronize()

    分组是为了在生成有限元空间时指定哪部分边界对应的是什么,比如说指定强制性的 Dirichlet 边界, 这个功能用函数 setPhysicalGroup 来实现. 该函数有三个参数, 第一个参数用来限制参加分组的几何维度, 比如说值为 0 时表示对点分组, 值为 1 是对线分组, 值为 2 是为面分组, 第二个参数是参与分组的点/线/面/...所对应的标号, 第三个参数(可选)为该分组的编号(正整数). 以 Dirichlet 边值问题为例, 四条边均为同类边界条件, 因此将这四条边统归为一组, 然后还可用 setPhysicalName 为这个分组设置一个别名, 用于后续使用, setPhysicalName 同样具有三个参数, 第一个参数为要设置别名的分组所属的维度, 即 0(点)、1(线)、2(面)等, 第二个参数为该分组的原始编号, 十个正整数, 第三个参数是要设置的别名, 为字符串.

    boundary = gmsh.model.addPhysicalGroup(1, [l₁, l₂, l₃, l₄])   # bc 是个引用, 用来替代因没有显式给出二自动生成的标号
gmsh.model.addPhysicalGroup(2, [surface], 1)
gmsh.model.setPhysicalName(1, boundary, "bc")

    最后, 使用 gmsh.model.mesh.generate(2) 来生成二维图形的三角剖分(如图所示), 使用 gmsh.write("文件名") 保存该剖分的单元信息, 用 gmsh.finalize() 来关闭 API 接口.

    gmsh.model.mesh.generate(2)
gmsh.write("square.msh")
gmsh.finalize()


    请注意, 一旦建立对应维度的分组, Gmsh 只会默认到处该维度中属于某个分组的元素, 而没有分组的将不会导出, 为了让 Gmsh 强制导出所有元素, 可以使用 gmsh.option.setNumber("Mesh.SaveAll", 1).

Cartesian 网格

针对于一维区间/二维矩形/高维笛卡尔积形式的区域, Gridap 还内置了用来生成一致剖分(正规)笛卡尔网格的函数 CartesianDiscreteModel. 这个函数使用两个元组.$(a_1,b_2,a_2,b_2,\cdots,a_n,b_n)$ 和 $(J_1,J_2,\cdots,J_n)$ 作为传入参数, 分别是区域对应笛卡尔积 ($\prod_{k=1}^n (a_k,b_k)$) 中各方向的左右端点和各个方向的一致剖分单元数. 对于正方形 $\Omega=(0,1)\times (0,1)$, 代码及其生成的剖分如下.

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain, partition)

有限元空间

这一部分介绍在得到网格剖分之后怎样生成相对应的有限元空间.

导入 Gmsh 网格剖分

如果上一步网格剖分是由 Gmsh 生成的, 那么我们还需要安装库/插件 GridapGmsh 中的 GmshDiscreteModel 来导入 square.msh 文件使得 Gridap 能够处理. 直接使用 using Pkg; Pkg.add("GridapGmsh") 安装即可.

using GridapGmsh
model = GmshDiscreteModel("square.msh")

协调 Lagrange 有限元空间

有了离散模型之后, 就可以生成其所对应的有限元空间了. 对于通常的齐次 Dirichlet 边界条件的有限元空间, 首先用 ReferenceFE(lagrangian, Float64, order) 指定参考单元上的基函数类型, 具体含义类型为 order 次的 Lagrange 基函数, Float64 表示存储类型为 64 位浮点数的标量函数/空间, 与之相对应的为向量函数/空间, 比如说 Navier-Stokes 方程所对应的函数空间. 然后通过构造器 TestFESpace(model, reffe; conformity=:H1, dirichlet_tags="") 得到(齐次边界条件)测试函数空间, 这里, 前两个参数的含义是生成以参考单元上的基函数类型 reffe 为模板的离散模型 model 上的测试函数空间; 关键字参数 conformity 用来指定空间的正则性, 即生成的有限元空间为 $H^1(\Omega)$ 的子空间; 因为 Dirichlet 边界条件是强制性边界条件, 需在相应的有限元空间中显式给出, 而关键字参数 dirichlet_tags 就是用来指定 Dirichlet 边界是哪一部分的, 具体传入的参数为 Gmsh 中利用 setPhysicalName 为相应边界设定的别名或者 CartesianDiscreteModel 默认的全部边界 boundary. 这一部分的完整代码为

order = 1
reffe = ReferenceFE(lagrangian, Float64, order)
V = TestFESpace(model, reffe; conformity=:H1, dirichlet_tags="bc")

有了 V 之后, 就可以利用该测试函数空间得到试用函数空间, 相应的代码为 TrialFESpace(V, x->0), 第一个参数为测试函数空间, 第二个参数为一个函数的引用/句柄, 用来指定边值, 这里使用的是匿名函数的语法, 当然也可以写作 TrialFESpace(V, 0) 或者 g(x)=0; TrialFESpace(V,g).

U = TrialFESpace(V, x->0)

三角化和数值积分

在 Gridap 中, 是使用类型 Triangulation 将离散模型 model (三角化)得到积分区域(记为 Ω)的, 进而, 使用 Measure(Ω, degree) 得到数值积分中积分单元的 Lebesgue 测度, degree 用来指定所用数值积分的精度.

degree = 2  # 对于分片线性多项式组成的有限元空间, 2阶精度的数值积分已经够用了
Ω = Triangulation(model)
dΩ = Measure(Ω, degree)

弱形式与有限元问题

经过上述一系列准备工作后, 我们可以正式定义有限元格式的弱形式及其所对应的线性系统了. 在 Gridap 中, 弱形式的写法是很数学化的:

f(x) = 2 * π^2 * sin(π*x[1]) * sin(π*x[2])  # 源函数
a(u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ  # 双线性型
b(v) = ∫(v*f) * dΩ  #  右端项

这里的右端项

$$ \begin{equation} f(x) = 2\pi^2\sin{(\pi x_1)}\sin{(\pi x_2)}. \end{equation} $$

然后通过 AffineFEOperator(a, b, U, V) 组合成线性系统, 其中四个参数分别为双线性型、右端项、试用函数空间和测试函数空间. 下面就通过求解算法为 LU 分解 (LUSolver()) 线性有限元求解器 LinearFESolver(LUSolver()) 来得到有限元格式的离散解 uₕ. 当然, 也可以用 get_matrix(op) 和 get_vector(op) 来得到对应的刚度矩阵和右端向量, 使用自己的求解器求解.

op = AffineFEOperator(a,b,U,V)
solver = LinearFESolver(LUSolver())
uₕ = solve(solver, op)

得到的数值解为

最终代码

import gmsh
using Gridap, GridapGmsh

"""
    generate_mesh(height::Float64=1.0, width::Float64=1.0, lc::Float64=0.05, 
                dirichlet_tags::String="boundary", path::String="./", name::String="model") :: String 
调用 Gmsh API 生成第一象限中原点为左下顶点、高度为 `height`、宽度为 `width` 的矩形区域的三角网格, 单元尺寸为 `lc`, 
Dirichlet 边界的标签为 `dirichlet_tags`, 并且将网格数据保存在路径 `path` 下的名为 `name` 且后缀名为 `.msh` 的文件中.
"""
function generate_mesh(height::Float64=1.0, width::Float64=1.0, lc::Float64=0.05, dirichlet_tags::String="boundary", path::String="./", name::String="model") :: String 
    # 初始化
    gmsh.initialize()
    # 新建名为 name 的模型
    gmsh.model.add(name)
    # 添加节点
    p₁ = gmsh.model.geo.addPoint(0.0, 0.0, 0.0, lc)
    p₂ = gmsh.model.geo.addPoint(width, 0.0, 0.0, lc)
    p₃ = gmsh.model.geo.addPoint(width , height, 0.0, lc)
    p₄ = gmsh.model.geo.addPoint(0.0, height, 0.0, lc)
    # 按逆时针将节点连接成线
    l₁ = gmsh.model.geo.addLine(p₁, p₂)
    l₂ = gmsh.model.geo.addLine(p₂, p₃)
    l₃ = gmsh.model.geo.addLine(p₃, p₄)
    l₄ = gmsh.model.geo.addLine(p₄, p₁)
    # 将线按逆时针组合成闭合曲线形成一个面
    curve = gmsh.model.geo.addCurveLoop([l₁, l₂, l₃, l₄])
    surface = gmsh.model.geo.addPlaneSurface([curve])
    # 同步化, 生成数据结构来存储模型中点线面的信息
    gmsh.model.geo.synchronize()
    # 指定物理边界
    dirichlet = gmsh.model.addPhysicalGroup(1, [l₁, l₂, l₃, l₄])
    gmsh.model.addPhysicalGroup(2, [surface], 1)
    gmsh.model.setPhysicalName(1, dirichlet, dirichlet_tags)
    # 生成(2d)网格
    gmsh.model.mesh.generate(2)
    # 存储网格
    file_path = string(path, name, ".msh")
    gmsh.write(file_path)
    gmsh.finalize()
    return file_path
end

"""
    PoissonSolver(mesh::String, f::Function, dirichlet_tags::String="boundary")
网格剖分为 `mesh`, 右端项为 `f` 且 Dirichlet 边界为 `dirichlet_tags` 的齐次 Dirichlet 边界的 Poisson 方程的求解器. 数值解存储在文件 `result.vtu` 中.
"""
function PoissonSolver(mesh::String, f::Function, dirichlet_tags::String="boundary")
    # 导入网格
    model = GmshDiscreteModel(mesh)
    # 基函数的次数
    order = 1
    # 参考单元上的基函数类型
    reffe = ReferenceFE(lagrangian, Float64, order)
    # 测试函数空间
    V = TestFESpace(model, reffe; conformity=:H1, dirichlet_tags=dirichlet_tags)
    # 试用函数空间
    U = TrialFESpace(V, x->0)
    # 数值积分的精度
    degree = 2
    # 将离散模型三角化、可测化
    Ω = Triangulation(model)
    # 积分单元的测度
    dΩ = Measure(Ω, degree)
    # 双线性型
    a(u, v) = ∫(∇(u) ⊙ ∇(v)) * dΩ 
    # 右端项
    b(v) = ∫(v*f) * dΩ  
    # 组装得到线性系统
    op = AffineFEOperator(a,b,U,V)
    # 求解器
    solver = LinearFESolver(LUSolver())
    # 数值解
    uₕ = solve(solver, op)
    writevtk(Ω, "result", cellfields=["uh"=>uₕ])
    return uₕ
end

f(x) = 2 * π^2 * sin(π*x[1]) * sin(π*x[2])
mesh = generate_mesh()
uₕ= PoissonSolver(mesh, f)
最后修改：2023 年 09 月 14 日
© 允许规范转载