UFL：有限元形式语言》表达式【翻译】

2021-02-01

| 科学计算 | | FEniCS , 有限元法 , 偏微分方程 |

17·4 定义表达式

UFL的大多数内容涉及如何声明表达式，例如公式17.13中的积分表达式。最基本的表达式是终值，它们不依赖于其他表达式。其他表达式称为算符，将在17.4.2-17.4.5节中讨论。

UFL中的终值类型包括形式参数（这是第17.4.1节的主题），几何量和字面常量。字面常量中有标量整数和浮点值，以及dxd单位矩阵I = Identity(d)。要获得单位向量，只需使用单位矩阵的行或列，例如e0 = I[0,:]。同样，I[i,j]表示Kronecker $\delta$ 函数 $\delta_{ij}$ （有关索引记号的详细信息，请参见第17.4.2节）。可用的几何值是空间坐标x ↔ cell.x和维面法线n ↔ cell.n。几何维数可由cell.d获得。

17·4·1 形式的参数

基函数和系数函数分别由参数和系数表示。形式的参数顺序由UFL代码中声明的形式参数的顺序决定。每个基函数参数代表其有限元空间的基底中的任何函数

\phi^j \in \{\phi_k^j\}, \qquad V_h^j = \mathrm{span} \left\{φ_k^j\right\} \tag{17.14}

目的是稍后针对所有 $\phi_k$ 计算如公式（17.12）所示的形式。每个系数函数 $$w$$ 代表有限元空间 $$V_h$$ 中的离散函数；通常是带系数 $$w_k$$ 的基函数 $\phi_k\in V_h$ 加权和

w =\sum^{|V_h|}_{k=1}w_k\phi_k \tag{17.15}

例外情况是系数函数只能逐点求值，这些函数使用带有“正交”族的有限元声明。基函数被声明为任意单元，如下所示：

# UFL code
# phi = Argument(element) 
phi = Argument(element, 0) # 较新的版本中至少要有两个输入； 0表示关联首个子单元
v = TestFunction(element)
u = TrialFunction(element)

通过使用TestFunction和TrialFunction声明代替Argument，您可以忽略它们的相对顺序。对于元数为r > 2的形式，只能使用Argument。

对任意单元，可类似地声明系数函数，并且存在用于声明常系数的简写记号：

# UFL code
w = Coefficient(element)
c = Constant(cell)
v = VectorConstant(cell)
M = TensorConstant(cell)

如果需要混合有限元空间 $V_h = V_h^0 \times V_h^1$ 中的形式参数 $$u$$ ，但是使用子函数 $u_0\in V_h^0$ 和 $u_1\in V_h^1$ 更容易表示此形式，则可以使用split的通用方式，将混合函数或基函数分解到其子函数中：

# UFL code
V = V0*V1
# V = element * element  # 由于V0和V1前面没给赋值，不妨用这个V代替测试
u = Coefficient(V)
u0, u1 = split(u)

split函数可以处理任意混合单元。另外，一个方便的速记符号表示法是在参数声明后加上分隔符

# UFL code
v0, v1 = TestFunctions(V)
u0, u1 = TrialFunctions(V)
f0, f1 = Coefficients(V)

17·4·2 索引记号

UFL允许同时使用张量代数和索引记号来处理任意秩的张量表达式。假设大家基本熟悉张量代数和索引记号。这里的重点是如何在UFL中表示索引记号。

假设标准正交欧几里德基底 $\{e_k \in \mathbb{R}^d\}^d_{k=1}$ ，那么向量可表示成此基底下的标量分量。二阶张量可以用标量分量来表示，对应的基底是 $\left\{e_i \otimes e_j\right\}^d_{i, j=1}$ 。任意阶张量可以用相同的方式表示，如此所示。

v =\sum^d_{k=1}v_k e_k \tag{17.16}

A =\sum^d_{i=1}\sum^d_{j=1}A_{ij}e_i \otimes e_j \tag{17.17}

C =\sum^d_{i=1}\sum^d_{j=1}\sum^d_{k=1} C_{ijk}e_i \otimes e_j \otimes e_k \tag{17.18}

在这里， $$v$$ ， $$A$$ 和 $$C$$ 分别是1、2和3阶张量。如果索引没有赋值（例如 $$v_i$$ 中的 $$i$$ ），则称其为自由的；如果索引具有固定值（例如 $$v_1$$ 中的 $$1$$ ），则称其为固定的。带有自由索引的表达式表示您可以通过将固定值分配给索引来获得任意表达式。表达式 $A_{ij}$ 是标量值，表示欧几里德基底下的张量 $$A$$ 的任意分量 $$ (i, j)$$ 。在纸上工作时，很容易在张量记号（ $$A$$ ）和索引记号（ $A_{ij}$ ）之间切换，因为知道张量及其组成部分是相同物理量的不同表示。在编程语言中，我们必须表示成一个操作，从张量到标量分量的映射，并显式返回。从张量到其分量的映射，对于2阶张量可定义为

A_{ij} = A : (e_i \otimes e_j) \tag{17.19}

这使用索引记号A[i,j]来做到。从分量值 $A_{ij}$ 定义张量 $$A$$ 为

A = A_{ij}e_i \otimes e_j \tag{17.20}

并使用函数as_tensor(Aij, (i,j))来做到。为了说明这点，考虑两个向量的外积 $A = u \otimes v = u_i v_j e_i \otimes e_j$ ，以及相应的标量分量 $A_{ij}$ 。一种实现方法是

# UFL code
A = outer(u, v)
Aij = A[i, j]

或者，可以使用索引记号直接表示A的分量，例如 $A_{ij} = u_i v_j$ 。然后可以通过以下方式将 $A_{ij}$ 映射到 $$A$$ ：

# UFL code
Aij = v[j]*u[i]
A = as_tensor(Aij, (i, j))

这两对代码在数学上是等效的，每对代码的结果都是变量A代表张量 $$A$$ ，变量Aij代表张量 $A_{ij}$ 。请注意，自由索引没有顺序，因此它们在表达式v[j]*u[i]中的出现顺序微不足道。可以使用专用函数as_vector和as_matrix代替as_tensor。尽管上面的示例使用了2阶张量，但是映射可推广到任意阶张量。

在为表达式建立索引时，还可以使用固定索引，例如在表示单个标量分量的A[0,1]中。固定索引也可以与自由索引混合，例如A[0,i]。另外，可以使用切片代替索引。使用切片的示例是A[0,:]，它是表示A的第0行的向量表达式。要创建新索引，您可以创建一个索引或一次创建多个索引：

# UFL code
i = Index()
j, k, l = indices(3)

有一组预定义的索引i, j, k, l和p, q, r, s，这些索引对于大多数应用来说应该是足够的。

如果您的分量不是表示为具有自由索引的表达式，而是表示为独立的不相关的标量表达式，那么可以使用as_tensor及其对等的函数来构建张量。例如，让我们定义一个2D旋转矩阵，可将向量表达式旋转 $\frac{\pi}{2}$ ：

# UFL code
th = pi/2
A = as_matrix([[ cos(th), -sin(th)],
               [ sin(th), cos(th)]])
u = A*v

当在一项中出现重复索引时，根据爱因斯坦约定，隐含对这些索引的求和。特别是，索引化二阶或更高阶张量（A[i,i]），微分一个带自由索引的表达式 (v[i].dx(i))，或者将两个具有共享自由索引的表达式相乘(u[i]*v[i])，这三种情况都可出现重复指标。

A_{ii} \equiv \sum_i A_{ii}, \qquad v_i u_i \equiv \sum_i v_i u_i, \qquad v_{i, i} \equiv \sum_i v_{i, i} \tag{17.21}

表达式Aij = A[i,j]在内部使用Indexed类表示。 Aij是A的引用，并保持原始张量表达式A的表示不变。隐式求和被显式表示成使用IndexSum类的表达式。由于这种显式表示，许多算法变得更容易实现，因为Product实例永远不能隐式表示成求和。有关表示类的更多详细信息，请参见第17.6节。

17·4·3 代数算符和函数

UFL定义了可用于组成表达式的一组全面的算符。基本的代数算符+, -, *, /可以在大多数UFL表达式之间使用，但有一些限制。除法要求分母是没有自由索引的标量表达式。求和操作各项必须具有相同的形状和一组自由索引。

乘法算符*只能在两个标量，一个标量和任何张量，一个矩阵和一个向量以及两个矩阵之间才有效。对于一些极少数情况，我们可以明确地使用张量代数算符和索引记号来定义其他乘积。两个具有共享自由索引表达式的乘积意味着对这些索引求和，有关索引记号的更多信息，请参见第17.4.2节。

三个常用算符是dot(a, b)，inner(a, b)和outer(a, b)。两个任意阶张量的点积是第1个张量的最后一个索引与第2个张量的第一个索引的求和。有一些例子：

v \cdot u = v_i u_i \tag{17.22}

A \cdot u = A_{ij} u_j e_i \tag{17.23}

A \cdot B = A_{ik} B_{kj} e_i e_j \tag{17.24}

C \cdot A = C_{ijk} A_{kl}e_i e_j e_l \tag{17.25}

内积是对所有索引求和，例如 $v : u = v_i u_i \tag{17.26}$

A : B = A_{ij} B_{ij}, \tag{17.27}

C : D = C_{ijkl}D_{ijkl} \tag{17.28}

外积的一些示例 $v \otimes u = v_i u_j e_i e_j \tag{17.29}$

A \otimes u = A_{ij} u_k e_i e_j e_k \tag{17.30}

A \otimes B = A_{ij} B_{kl} e_i e_j e_k e_l \tag{17.31}

其他常见的张量代数算符有cross(u,v)，transpose(A)（或A.T），tr(A)，det(A)，inv(A)，ofac(A)，dev（A），skew(A)和sym(A)。这些张量代数算符中的大多数，都要求作用到没有自由索引的张量。这些算符的详细定义可以在手册中找到。

作用于无自由索引标量表达式的一组通用基本函数算符，比如：abs(f)，pow(f, g)，sqrt(f)，exp(f)，ln(f)，cos(f)， sin(f)，tan(f)，acos(f)，asin(f)，atan(f)和sign(f)。任何携带标量参数的算符，都可逐元地应用到张量上，比如 elem_op(sin, A)。

17·4·4 微分算符

UFL实现了关于三种不同的变量的导数。最常用的一种是空间导数。表达式也可以求关于任意用户定义变量的微分。并且最后一种导数是形式或函数关于离散函数系数的导数；即关于系数或常数。形式导数在第17.5.1节中说明。

请注意，导数在声明时是不会立即计算的。有关如何计算导数的讨论，请参见第17.7节。

空间导数 基本空间导数 $\frac{\partial f}{\partial x_i}$ 可以用两种等效的方式表示：

# UFL code
df = Dx(f, i)
df = f.dx(i)

在此，df表示 $$f$$ 在空间方向 $$x_i$$ 上的导数。索引 $$i$$ 可以是代表一个固定空间方向 $$x_i$$ 上以示区别的整数，也可以是代表以示方向差异的自由索引。符号f.dx(i)旨在反映索引记号 $f_{,i}$ ，它是 $\frac{\partial f}{\partial x_i}$ 的简写。重复索引意味着求和，这样向量值表达式v的散度可以写成 $v_{i, i}$ 或v[i].dx(i)。

定义了几种常见的复合空间导数算符，即梯度(gradient)，散度(divergence)和旋度算符(curl)。这些算符分别命名为grad，div，nabla_grad，nabla_div，curl和rot（rot是curl的同义词）。请注意，定义梯度和散度有两种常用方法，并且UFL同时支持之。

令 $$s$$ 为标量表达式， $$v$$ 为向量表达式， $$M$$ 为r阶张量表达式。在UFL中，算符的grad被明确定义为

(\mathrm{grad}(s))_i = s_{,i} \tag{17.32}

(\mathrm{grad}(v))_{ij} = v_{i,j} \tag{17.33}

(\mathrm{grad}(M))_{i_1\dots i_r k} = M_{i_1\dots i_r, k} \tag{17.34}

并且算符div相应地定义为

\mathrm{div}(v) = v_{i, i} \tag{17.35}

(\mathrm{div}(M))_{i_1 \dots i_{r−1}} = M_{i_1\dots i_r, i_r} \tag{17.36}

相反，nabla_*算符的定义依赖 $\nabla$ 算符

∇ \equiv e_k \frac{\partial}{\partial x_k} \tag{17.37}

算符nabla_grad是 $\nabla$ 与其操作对象的外积：

(\nabla s)_i = s_{,i} \tag{17.38}

(\nabla v)_{ij} = v_{j,i} \tag{17.39}

(\nabla M)_{k,i_1 \dots i_r} = M_{i_1 \dots i_r, k} \tag{17.40}

同样，算符nabla_div是 $\nabla$ 与其操作对象的点积：

\nabla \cdot v = v_{i, i} \tag{17.41}

(\nabla \cdot M)_{i_2 \dots i_r} = M_{i_1 \dots i_r, i_1} \tag{17.42}

从值形状的角度考虑，grad算符将一个轴附加到其操作对象的形状的末尾，而nabla_grad算符则是预先就选定了轴。对于标量梯度，结果是相同的。相应地，div算符是与其操作对象关于最后一个索引求和，而nabla_div算符是与其操作对象关于第一个索引求和。对于向量的散度，也是相同的结果。

算符curl和rot按惯例是没有区别的。对3维向量表达式，可以由 $\nabla$ 算符和叉积来定义旋度：

\begin{aligned}\mathrm{curl}(v) & \equiv \nabla \times v \\ &=e_0(v_{2,1} − v_{1,2}) − e_1(v_{2,0} − v_{0,2}) + e_2(v_{1,0} − v_{0,1})\end{aligned} \tag{17.43}

对于2维向量和标量表达式，定义为：

\mathrm{curl}(v) \equiv v_{1,0} − v_{0,1} \tag{17.44}

\mathrm{curl}(f) \equiv f_{,1}e_0 − f_{,0}e_1 \tag{17.45}

用户定义变量 第二类微分变量是用户定义的变量，可以表示任意表达式。关于任意量的自动求导，对很多任务都是有用的，比如：从物性定律微分到计算灵敏度。可以将任意表达式 $$g$$ 指定为一个变量 $$v$$ 。被定义为 $$v$$ 的函数的表达式 $$f$$ ，可对 $$f$$ 作关于 $$v$$ 的微分:

v = g \tag{17.46}

f = f (v) \tag{17.47}

h(v) = \frac{\partial f (v)}{\partial v} \tag{17.48}

设 $g = \sin(x_0)$ 和 $f = e^{v^2}$ ，得到 $h = 2v e^{v^2} = 2 \sin(x_0)e^{\sin^2(x_0)}$ ，可以实现如下：

# UFL code
# 在较新的版本中，cell没有x这个成员，改用随后的两行代码
# g = sin(cell.x[0])  
x = SpatialCoordinate(cell)
g = sin(x[0])
v = variable(g)
f = exp(v**2)
h = diff(f, v)

尝试在Python会话中运行此代码并打印表达式。结果是

# Python code
»>print(v)
var0(sin((x)[0]))
»>print(h)
d/d[var0(sin((x)[0]))] (exp((var0(sin((x)[0]))) ** 2))

注意，该变量有标签“var0”，并且h仍表示抽象导数。第17.7节解释了如何计算导数。

17·4·5 其他算符

还提供了一些算符来实现不连续Galerkin方法。基本概念是将表达式限制在内部维面的正侧或负侧，分别简单地表示为v("+")或v("-")。最重要的是，实现了算符avg和jump，被定义为

\mathrm{avg}(v) =\frac{1}{2}(v^+ + v^−) \tag{17.49}

\mathrm{jump}(v) = v^+ − v^− \tag{17.50}

这些算符只能被用于内部维面上的积分（*dS）。

UFL中包含的唯一控制流结构是条件表达式。条件表达式取两个值之一，具体取决于布尔逻辑表达式的结果。语法是

# UFL code
f = conditional(condition, true_value, false_value)

解释为

f = \left\{\begin{aligned}\mathrm{true\_value}, &\qquad 如果 \mathrm{condition}为真 \\ \mathrm{false\_value} , &\qquad 否则 \end{aligned} \right. \tag{17.51}

条件可以是以下之一


`lt(a, b)` $\leftrightarrow (a < b)$	`gt(a, b` $\leftrightarrow (a > b)$
`le(a, b)` $\leftrightarrow (a \le b)$	`ge(a, b)` $\leftrightarrow (a \ge b)$
`eq(a, b)` $\leftrightarrow(a = b)$	`ne(a, b)` $\leftrightarrow (a \ne b)$
`And(P, Q)` $\leftrightarrow (P \wedge Q)$	`Or(P, Q)` $\leftrightarrow (P \vee Q)$
`Not(P)` $\leftrightarrow (\neg P)$

【章节目录】