title: "DOLFIN：C++/Python有限元库》变分形式/组装/边界条件【翻译】"
date: 2021-01-11T09:27:35+08:00
draft: true
categories: ["科学计算"]
tags: ["FEniCS","有限元法","偏微分方程"]

【[[docs/fem/_index#第十章 DOLFIN：C＋＋／Python有限元库|章节目录]]】

10·3·7 变分形式

DOLFIN依靠FEniCS工具链FIAT–UFL–FFC / SFC–UFC来计算有限元变分形式。 使用形式编译器FFC或SFC（第11章和第15章）之一来编译以UFL形式语言表示的变分形式（第17章），DOLFIN使用生成的UFC代码（第16章）来计算（组装）变分形式 。

UFL形式语言允许广泛的变分形式以接近数学符号的语言来表达，作为例证，下面的表达式定义（或部分定义）了线性弹性问题离散化的双线性和线性形式：

# UFL code
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx

这应该与相应的数学符号进行比较：

$$ a(u, v) = \int_\Omega \sigma(u) : \epsilon(v)dx \tag{10.3} $$

$$ L(v) = \int_\Omega f \cdot v dx \tag{10.4} $$

在此，$\epsilon(v) = (\mathrm{grad} \\ v + (\mathrm{grad} \\ v)^T)/2$表示对称梯度，$\sigma(v) = 2\mu \epsilon(v) + \lambda \mathrm{tr} \\ \epsilon(v) I$是应力张量。 有关UFL形式语言的详细介绍，请参阅第17章。

C++接口的用户必须通过在命令行上调用形式编译器来显式处理代码生成过程。 为了求解选择了特定参数值（Lamé常数$\mu$和$\lambda$）的线性弹性问题，用户可以在名为Elasticity.ufl的文件中键入以下代码：

# UFL code

V = VectorElement("Lagrange", tetrahedron, 1)

u = TrialFunction(V)
v = TestFunction(V)
f = Coefficient(V)

E = 10.0
nu = 0.3

mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)

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

可以使用符合UFL/UFC的形式编译器来编译此代码，以生成UFC C++代码。 例如，使用FFC：

# Bash code
ffc -l dolfin Elasticity.ufl

这将生成一个名为Elasticity.h的C++头文件（包括实现），该头文件可以包含在C++程序中，并用于实例化$a$和$L$这两种形式：

// C++ code
#include <dolfin.h>
#include "Elasticity.h"

using namespace dolfin;

int main()
{
    UnitSquare mesh(8, 8);
    Elasticity::FunctionSpace V(mesh);
    Elasticity::BilinearForm a(V, V);
    Elasticity::LinearForm L(V);
    MyExpression f; // code for the definition of MyExpression omitted
    L.f = f;

    return 0;

}

形式的实例化涉及其定义所在的FunctionSpace实例化。 在创建形式之后，任何出现在形式定义中的系数（此处为右侧f）必须被加上。 Python用户可依赖于自动代码生成，并直接将变分形式定义为Python脚本的一部分：

# Python code

from dolfin import *

mesh = UnitSquare(8, 8)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

u = TrialFunction(V)
v = TestFunction(V)
f = MyExpression() # code emitted for the definition of f

E = 10.0
nu = 0.3

mu = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)

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

这个脚本将触发自动代码生成，以定义FunctionSpace V。 两种形式$a$和$L$的代码生成被推迟到相应的离散算符（矩阵和向量）被组装的那一刻。

10·3·8 有限元组装

DOLFIN的核心功能是有限元变分形式的组装。 给定变分形式（$a$），DOLFIN组装出相应的离散算符（$A$）。 离散算符的组装服从第6章中描述的一般算法。 以下代码说明了如何分别从泛函（$M$），线性形式（$L$）和双线性形式（$a$）组装出标量（$m$），向量（$b$）和矩阵（$A$）：

// C++ code
Vector b;
Matrix A;

double m = assemble(M);
assemble(b, L);
assemble(A, a);

# Python code
m = assemble(M)
b = assemble(L)
A = assemble(a)

从Python接口组装变分形式会在运行时自动触发代码生成，编译和链接。 生成的代码将自动实例化并发送到DOLFIN C++编译器。 结果，从Python接口进行有限元组装与从C++接口进行组装同样高效，而处理自动代码生成的开销却很小。 生成的代码将被缓存以供以后重用，因此，相同形式的重复组装或两次运行相同程序不会重新触发代码生成。 相反，先前生成的代码会从缓存中自动加载。

DOLFIN提供了一种通用的组装算法，用于以任何形式组装任何等级的张量（标量，向量，矩阵等）。 这是可能的，因为组装算法依赖于GenericTensor接口，组装算法依赖于变体形式及其特定离散化的部分是在组装之前生成的，并且网格接口是维数无关的。 组装算法接受许多可选参数，这些参数控制在组装之前是否应重置已组装张量的稀疏性，以及在组装之前是否应将张量清零。 如果在特定子域（使用dx(0)， dx(1)等）上定义了形式，参数也可以支持网格的特定子域。

除了assemble函数外，DOLFIN还提供assemble_system函数，该函数可组装由双线性和线性形式组成的一对形式，并在组装过程中应用基本边界条件。 边界条件的应用作为assemble_system的调用的一部分，可保持所组装矩阵的对称性（请参见第6章）。

组装算法已针对使用MPI的分布式内存架构（集群）和使用OpenMP的共享内存架构（多核）进行了并行化。 第10.4节对此进行了更详细的讨论。

10·3·9 边界条件

DOLFIN可以处理Neumann（自然）边界和Dirichlet（基本）边界条件的应用。 自然边界条件通常通过问题的变分陈述来应用，而基本边界条件通常应用于离散的方程组。

自然边界条件。 自然边界条件通常是作为边界项出现的，这是将微分方程的部分乘以测试函数再积分得到的结果。 作为一个简单的例子，我们考虑线性弹性变分问题。 支配弹性体位移的偏微分方程可以表示为

$$ \begin{aligned} - \mathrm{div} \\ \sigma(u) &= f \qquad \mathrm{in} \\ \Omega \\ \sigma \cdot n &= g \qquad \mathrm{on} \\ \Gamma_N \subset \partial \Omega \\ u &= u_0 \qquad \mathrm{on} \\ \Gamma_D \subset \partial \Omega \end{aligned} \tag{10.5} $$

其中，$u$是要计算的未知位移场，$\sigma(u)$是应力张量，$f$是给定的力密度，$g$是在边界的一部分$\Gamma_N$上的给定牵引，$u_0$是在边界的一部分$\Gamma_D$上的给定位移。 乘以测试函数v并分部积分，我们得到

$$ \int_\Omega \sigma(u) : \epsilon(v)dx − \int_{\partial \Omega}(\sigma \cdot n) \cdot v ds = \int_\Omega f \cdot v dx \tag{10.6} $$

这里，我们利用$\sigma(u)$的对称性，可将$\mathrm{grad} \\ v$替换为对称梯度$\epsilon(v)$。 由于位移$u$在Dirichlet边界$\Gamma_D$上是已知的，因此让$\Gamma_D$上的$v=0$。 此外，我们用边界$\Gamma_N$其余部分（Neumann）上的给定牵引$g$替换$\sigma \cdot n$以获得

$$ \int_\Omega \sigma(u) : \epsilon(v)dx = \int_\Omega f \cdot v dx + \int_{\partial \Omega}(\sigma \cdot n) \cdot v ds\int_{\partial \Omega}g \cdot v ds \tag{10.7} $$

以下代码演示了，作为.ufl文件的一部分或Python脚本的一部分，如何以UFL形式语言实现此变分问题：

# UFL code
a = inner(sigma(u), sym(grad(v)))*dx
L = dot(f, v)*dx + dot(g, v)*ds

要指定边界积分dot(g, v)*ds只沿Neumann边界$\Gamma_N$求值，必须指定边界的哪一部分包括在$ds$积分中。 如果只有一个Neumann边界，则可以简单地将$ds$积分写为，包括Dirichlet边界在内的，整个边界的积分，因为将沿着Dirichlet边界将测试函数v设置为零。

如果存在多个Neumann边界条件，则必须使用FacetFunction来指定Neumann边界。 此FacetFunction必须为Mesh的每个维面指定其所属的边界部分。 对于当前示例，一种适当的策略是用0标记Neumann边界上的每个维面，并用1标记所有其他维面（包括域内部维面）。 这可以通过许多不同的方式来完成。 一种简单的方法是使用程序MeshBuilder，并以图形方式标记Mesh的各个维面。 另一个选择是通过DOLFIN类SubDomain。 以下代码说明了如何将$x = 0.5$左侧的所有边界维面标记为第一Neumann边界，而将所有其他边界维面标记为第二Neumann边界。 注意，用on_boundary作为参数，提供给DOLFIN的inside函数。 此参数表明一个点是否位于$\Omega$的边界$\partial \Omega$上，这使我们可以仅标记$x = 0.5$左侧边界维面。 还要注意，使用DOLFIN_EPS可以确保由于有限精度算术而包含的点可能恰好位于$x = 0.5$的右侧。

// C++ code

class NeumannBoundary : public SubDomain
{
    bool inside(const Array<double>& x, bool on_boundary) const
    {
        return x[0] < 0.5 + DOLFIN_EPS && on_boundary;
    }
};

NeumannBoundary neumann_boundary;
FacetFunction<uint> exterior_facet_domains(mesh);
exterior_facet_domains.set_all(1);
neumann_boundary.mark(exterior_facet_domains, 0);

# Python code
class NeumannBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] < 0.5 + DOLFIN_EPS and on_boundary

neumann_boundary = NeumannBoundary()
exterior_facet_domains = FacetFunction("uint", mesh)
exterior_facet_domains.set_all(1)
neumann_boundary.mark(exterior_facet_domains, 0)

当与使用ds(0)和ds(1)定义的积分组合时，这些积分将分别对应于x = 0.5左侧和x = 0.5右侧所有维面边界上的积分。

一旦将边界指定为FacetFunction，就可以使用该对象来定义相应的积分域。 这在C++和Python中是不同的。 对C++，必须将ds成员变量分配给相应形式的：

// C++ code
a.ds = exterior_facet_domains;
L.ds = exterior_facet_domains;

除了根据ds成员变量指定的external_facet_domains之外，可以类似地使用dx成员变量指定cell_domains和使用dS变量指定interior_facet_domains。 请注意，不同形式可能会使用不同的边界定义。 对Python，可以通过下标简单地将边界定义连接到相应的测度：

# Python code
dss = ds[neumann_boundary]
a = ... + g*v*dss(0) + h*v*dss(1) + ...

是否正确指定边界是常见的错误来源。 为了调试边界条件的详情，可通过将FacetFunction写入VTK文件（请参阅文件I/O部分）或使用plot命令来绘制指定边界标记的FacetFunction可能会有所帮助。 使用plot命令时，绘图会显示插值到Mesh顶点的维面值。 因此，在这种情况下，必须小心解释靠近区域边界（角）的图。 在VTK输出中不存在该问题。

基本边界条件。 基本边界条件的应用是由DirichletBC类处理的。 使用此类，可以根据FunctionSpace，Function或Expression以及子域来指定Dirichlet边界条件。 可以根据SubDomain对象或FacetFunction的各项来指定子域。 DirichletBC指定的解应等于给定子域上的给定值。

以下代码示例说明了，对弹性问题 (10.5)，在Dirichlet边界$\Gamma_D$（此处假定为$x = 0.5$右侧的边界的一部分）上，如何使用SubDomain类来定义Dirichlet条件$u(x) = u_0(x) = \sin x$ 。 或者，可以使用FacetFunction指定子域。

// C++ code
class DirichletValue : public Expression
{
    void eval(Array<double>& values, const Array<double>& x) const
    {
        values[0] = sin(x[0]);
    }
};

class DirichletBoundary : public SubDomain
{
    bool inside(const Array<double>& x, bool on_boundary) const
    {
        return x[0] > 0.5 - DOLFIN_EPS && on_boundary;
    }
};

DirichletValue u_0;
DirichletBoundary Gamma_D;

DirichletBC bc(V, u_0, Gamma_D);

# Python code
class DirichletValue(Expression):
    def eval(self, value, x):
        values[0] = sin(x[0])

class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] > 0.5 - DOLFIN_EPS and on_boundary

u_0 = DirichletValue()
Gamma_D = DirichletBoundary()

bc = DirichletBC(V, u_0, Gamma_D)

Python用户还可以使用以下紧凑语法：

# Python code
u_0 = Expression("sin(x[0])")
bc = DirichletBC(V, u_0, "x[0] > 0.5 && on_boundary")

为了加快Dirichlet边界条件的应用，Python接口的用户还可以使用函数compile_subdomains。 有关详细信息，请参阅DOLFIN程序员参考。

可以将Dirichlet边界条件应用于线性系统或与Function相关的自由度矢量，如以下代码示例所示：

// C++ code
bc.apply(A, b);
bc.apply(u.vector());

# Python code
bc.apply(A, b)
bc.apply(u.vector())

将Dirichlet边界条件应用于线性系统将识别，应设为给定值的，所有自由度，并对线性系统进行修改，以使其解满足边界条件。 这是通过将与Dirichlet值相对应的矩阵行对角线归零并插入1，并将Dirichlet值插入右侧矢量的相应项中来实现的。 边界条件的这种应用不能保持对称性。 如果需要对称性，则可以选择使用assemble_system函数，它在组装过程中对称地应用Dirichlet边界条件。

多个边界条件可以应用于单个系统或向量。 如果将两个不同的边界条件应用于相同的自由度，则最后应用的值将覆盖任何先前设置的值。

【[[docs/fem/_index#第十章 DOLFIN：C＋＋／Python有限元库|章节目录]]】