FE

# 使用 deal.II 解决有限元问题
<hr/>

# Remarks on FEM Programming with deal.II

</div>
</div>

# Part.1
# 从一个例子说起：
# 有限元方法求解 Poisson 方程

</div>
</div>

## Poisson 方程

我们先从一类特殊的 Poisson 方程入手：

> 在区域 $\Omega$ 上，求解：
> 
> $$
> -\Delta u = f \quad \text{in } \Omega
> $$
> 
> 边界条件：
> 
> $$
>  u = 0 \quad \text{on } \partial \Omega
> $$
>
> 这里，$f$ 是已知函数，$u$ 是我们要求解的函数。
>

我们的目的是找到 $u$ 在 $\Omega$ 上满足边界条件和上述方程的近似解。

## 弱解

**目标：**寻找函数 $u$（来自合适的函数空间，且边界值为零），使得上式对所有合适的测试函数 $\varphi$ 都成立。

这样就得到一个弱解了，因为我们不能直接比较函数，但可以比较对所有测试函数的投影（积分）是否相等。

例如：

$$
\int g = \int h , \quad \int x\cdot g = \int x\cdot h , \quad \int x^2 g = \int x^2 h , \quad \ldots \Rightarrow g = h
$$

如果对于所有的**测试函数 (Test Functions)** $\varphi(x)$，其积分都相等，则认为这是一个弱解。

## 弱形式

如果我们需要求解一个弱解，就将原方程 $-\Delta u = f$ 两边乘以一个测试函数 $\varphi$，并在求解域 $\Omega$ 上积分：

$$
-\int_{\Omega} \varphi \Delta u \, d\mathbf{x} = \int_{\Omega} \varphi f \, d\mathbf{x}
$$

**分部积分：**

对左侧项 $-\int_{\Omega} \varphi \Delta u \, d\mathbf{x}$ 使用分部积分法：

$$
\int_{\Omega} \nabla \varphi \cdot \nabla u \, d\mathbf{x} - \int_{\partial \Omega} \varphi \mathbf{n} \cdot \nabla u \, dS = \int_{\Omega} \varphi f \, d\mathbf{x}
$$

**边界项处理：**

由于原问题 $u$ 满足齐次 Dirichlet 边界条件 ($u=0$)，测试函数 $\varphi$ 也必须满足同样的条件 （$\varphi=0$ on $\partial \Omega$ ）。因此，边界积分项 $\int_{\partial \Omega} \varphi \mathbf{n} \cdot \nabla u \, dS$ 消失。

## 离散化

在计算机上无法找到无限维空间中的精确解 $u$。引入 [Galerkin 方法](#gotojlj)，要求测试函数空间与解函数空间完全相同，把测试函数缩小到有限维子空间中，用给定的基函数表示出解 $u$ 的近似。

**我们将其近似为有限维解 $u_h$：**

将解 $u$ 近似为一组 $N$ 个形函数（Shape Functions）$\{\varphi_{j=0\cdots N-1}\}$ 的线性组合：

$$
u_h(\mathbf{x}) = \sum_{j=0}^{N-1} U_j \varphi_j(\mathbf{x})
$$

其中 $U_j$ 是待定的展开系数，也就是我们要解的未知量。

（ 或在 Deal.II 中这个东西被称为自由度，Degrees of Freedom, DoFs ）

## Shape Functions

有限元形函数（Shape Functions）是定义在每个单元（Element）上的局部函数，通常具有以下特点：
- 每个形函数在其对应的单元内非零，在其他单元上为零。
- 在节点处，形函数的值为1，而在其他节点处为0。
- 形函数可以是线性的或高阶多项式等。

**一个直观理解：** $\varphi_j(x)$ 是像帐篷一样的形函数，在某个点为1，周围线性下降到0。整个近似解 $u_h(x)$ 就是对这些帐篷函数赋予不同高度的叠加。

</div>
<div class="col">

要确定 $u_h(x)$，我们只需要算出有限个系数（Coefficients）$U_j$ 。

</div>

## 建立线性系统

**记号：**$(f, g)= \int_{\Omega} f g \, d\mathbf{x}$

于是，离散解 $u_h$ 必须满足弱形式，但我们只需要它对有限个基函数 $\varphi_i$（作为测试函数）成立：
$$
(\nabla u_h, \nabla \varphi_i) = (\varphi_i, f), \quad i = 0, \dots, N-1
$$

将 $u_h$ 的表达式代入上式：
$$
\left(\nabla \left(\sum_{j=0}^{N-1} U_j \varphi_j\right), \nabla \varphi_i\right) = (\varphi_i, f)
$$
$$
\sum_{j=0}^{N-1} U_j (\nabla \varphi_j, \nabla \varphi_i) = (\varphi_i, f)
$$

## 建立线性系统

这可以写成一个线性代数方程组

$$
\mathbf{A} \mathbf{U} = \mathbf{F}
$$

刚度矩阵 (Stiffness Matrix) $A$:

$$
A_{ij} = (\nabla \varphi_j, \nabla \varphi_i) = \int_{\Omega} \nabla \varphi_j \cdot \nabla \varphi_i \, d\mathbf{x}
$$

右端向量 (Right Hand Side Vector) $F$:
$$
F_i = (\varphi_i, f) = \int_{\Omega} \varphi_i f \, d\mathbf{x}
$$

解向量 (Solution Vector) $U$:
$$
U = [U_0, U_1, \dots, U_{N-1}]^T
$$

## 求解线性系统

使用数值线性代数方法（如共轭梯度法、直接法等）求解线性系统 $\mathbf{A} \mathbf{U} = \mathbf{F}$，得到系数向量 $\mathbf{U}$。

一旦得到 $\mathbf{U}$，我们就可以构造出近似解 $u_h(\mathbf{x})$，并在整个域 $\Omega$ 上进行评估和可视化。

# Part.2
# 使用 deal.II 解决有限元问题
# 还是求解 Poisson 方程

</div>
</div>

我们试着解决这样的一个问题：

> 
> $$
> -\nabla\cdot(a(x,y) \nabla u(x,y)) = f(x,y), \quad (x,y) \in \Omega =[-2,2]^2
> $$
>
> $$
> u(x,y) = \frac{1}{4}\sin(x+y), \quad (x,y) \in \partial\Omega_D
> $$
>
> $$
> \frac{\partial u}{\partial n} = g(x,y), \quad (x,y) \in \partial\Omega_N
> $$
>
> $a(x,y)$ 是一个空间相关的系数，之后给定义。
> $f(x,y)$ 是一个函数，之后给定义。
>

## 给定问题使用 dealii 求解的大致步骤

> 对于给定的偏微分方程组，首先我们需要将其转换为弱形式（Weak Form）。这涉及到将原始方程乘以一个测试函数，并在求解域上进行积分。给出线性系统的矩阵形式。

对于给定的问题，我们首先将其转换为弱形式：

$$
\int_{\Omega} a(x,y) \nabla u \cdot \nabla \varphi \, d\mathbf{x} = \int_{\Omega} f \varphi \, d\mathbf{x}+ \int_{\partial\Omega_N} g \varphi \, dS
$$

其中 $\varphi$ 是测试函数.

离散化后，我们得到一个线性系统：

$$
\mathbf{A} \mathbf{U} = \mathbf{F}
$$

其中，矩阵 $\mathbf{A}$ 的元素为 $A_{ij} = \int_{\Omega} a(x,y) \nabla \varphi_j \cdot \nabla \varphi_i \, d\mathbf{x}$，右端向量 $\mathbf{F}$ 的元素为 $F_i = \int_{\Omega} f \varphi_i \, d\mathbf{x}+ \int_{\partial\Omega_N} g \varphi_i \, dS$。

## 现在开始可以写代码了

一般地，我们可以按照以下步骤来使用 deal.II 给出相应的定义求解上述问题：

1. `make_grid();` 创建网格，定义求解域 $\Omega$.
2. `setup_system();` 定义有限元空间，分配自由度，初始化矩阵和向量。
3. `assemble_system();` 组装系统矩阵 $\mathbf{A}$ 和右端向量 $\mathbf{F}$。
4. `solve();` 求解线性系统 $\mathbf{A} \mathbf{U} = \mathbf{F}$。
5. `output_results();` 输出结果，进行可视化。

# 敬请老师和同学们批评指正！

<hr/>

参考资料有：

[Wolfgang Bangerth's Video Lectures](https://www.math.colostate.edu/~bangerth/videos.html)

[The deal.II Library](https://dealii.org/current/doxygen/deal.II/)

如果您有任何意见或者建议，欢迎通过电子邮件联系我：

[yxyou@mail.bnu.edu.cn](mailto:yxyou@mail.bnu.edu.cn)

</div>
</div>

<p style="color: #ccccccff;">此页被有意地留作空白。</p>

## COPYRIGHT

**图片来源：**

除了在前页中标明的参考资料外，其他图片均来自于：

- [いらすとや](https://www.irasutoya.com/) ；
- [Wikipedia](https://www.wikipedia.org/) ；

</div>
<div class="col">

**其他说明：**

- 本站幻灯片及相关样式/模板部分改编自 TonyCrane 的 [slide-template](https://github.com/TonyCrane/slide-template) ；

</div>