泊松方程试探法，混合形式泊松方程之FEniCS求解
泊松方程试探法，混合形式泊松方程之FEniCS求解稳态情况下，根据能量守恒，通过任意闭合曲面流出的能量（这里指热量） 必须等于 曲面内热源发出的能量：而泊松方程实际就是傅里叶定律＋能量守恒的结果。如图：现在我们要将其改成一阶方程组，引入一个叫热通量（或简称 通量）的辅助矢量，不难得到：这个通量是有物理意义的，那就是所谓的傅里叶定律：热通量 与温度 的负梯度成比例。这里我们取而已。
【我已经发了同名西瓜视频，本文仅是视频中Markdown文件内容，不包括视频中讲解】
如果你希望和我一样Win10上安装FEniCS，建议采用Win10 WSL2 Ubuntu。
FEniCS讲座系列第1讲：有限元法求解偏微分方程之FEniCS入门讲解【上一讲】
第２讲要点
- 将高阶偏微分方程 变换成 1阶偏微分方程组
- 1阶偏微分方程组 变换成 1个变分方程
- 自然边界 和 基本边界
- 如何定义混合空间
- 如何自定义表达式
- FEniCS快捷可视化
以泊松方程为例，这是一个二阶方程：
现在我们要将其改成一阶方程组，引入一个叫热通量（或简称 通量）的辅助矢量，不难得到：
这个通量是有物理意义的，那就是所谓的傅里叶定律：热通量 与温度 的负梯度成比例。
这里我们取而已。
而泊松方程实际就是傅里叶定律＋能量守恒的结果。如图：
稳态情况下，根据能量守恒，通过任意闭合曲面流出的能量（这里指热量） 必须等于 曲面内热源发出的能量：
进而
由于闭合区域的任意选择性，所以在整个，下式成立
对应的边界条件则需要改成：
转换成变分形式第一个方程是标量方程，选择标量测试函数；第二个方程是矢量方程，选择矢量测试函数。
方程１× 方程２· ，然后对全域积分，得到：
参考这个关系，可以具体地定义测试函数和试探函数。
将所有测试函数元组的集合记作；类似地，将所有试探函数元组的集合记作：
请注意，这里和上一讲的不同之处在于， 上的第一类边界条件 进入了变分公式（现在是自然边界条件），而第二类边界条件条件则在上进入试探空间的定义（现在是基本边界条件）。
将上面的结果整理得到最终的变分方程：
改成标准变分问题表述：寻求 满足
求解这个混合形式泊松方程- (一个单位矩阵)
- (Dirichlet边界)
- (Neumann边界)
- (法向导数)
- (源项)
from dolfin import *
# 创建单位矩形网格
mesh = UnitSquareMesh(32 32)
# 定义有限元空间和混合空间
BDM = FiniteElement("BDM" mesh.ufl_cell() 1)
DG = FiniteElement("DG" mesh.ufl_cell() 0)
W = FunctionSpace(mesh BDM * DG)第二步：定义已提供的表达式
# 第一类界值
u0 = Constant(0.0)
# 定义源函数
f = Expression("10*exp(-(pow(x[0] - 0.5 2) pow(x[1] - 0.5 2)) / 0.02)" degree=2)
# 定义满足G·n = -g的函数G(实际对应公式中的σ)
# g = sin(5x)
class BoundarySource(UserExpression):
def __init__(self mesh **kwargs):
self.mesh = mesh
super().__init__(**kwargs)
def eval_cell(self values x ufc_cell):
cell = Cell(self.mesh ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = sin(5*x[0])
values[0] = -g*n[0]
values[1] = -g*n[1]
def value_shape(self):
return (2 )
G = BoundarySource(mesh degree=2)第三步：创建和应用基本边界条件
# 定义基本边界条件
def boundary(x):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
bc = DirichletBC(W.sub(0) G boundary)第四步：定义变分问题
# 定义试探函数和测试函数
(sigma u) = TrialFunctions(W)
(tau v) = TestFunctions(W)
# 定义变分形式
n = FacetNormal(mesh)
a = (div(sigma)*v dot(sigma tau) - div(tau)*u)*dx
L = f*v*dx - u0*dot(n tau)*ds第五步：求解并绘图
# 求解
w = Function(W)
solve(a == L w bc)
(sigma u) = w.split()
# 绘图
plot(sigma)
plot(u)ParaView可视化
# 将解保存为VTK格式
vtkFile = File("mixed_poisson/sigma.pvd")
vtkfile << sigma
vtkfile = File("mixed_poisson/u.pvd")
vtkfile << uFEniCS快捷可视化
# 绘制矢量场图(如果是矢量场的话)
# plot(sigma mode="glyphs")
plot(sigma)
# 绘制位移图(如果是矢量场的话)
plot(sigma mode="displacement")
# 等高色图(如果是标量场的话)
# plot(u mode = "contourf" levels = 40)
# plot(u mode = "color" shading = "gouraud")
plot(u)
# 标量场3d图（曲面）
plot(u mode="warp" linewidths=0)
# 等高线(如果是标量场的话)
plot(u mode="contour")
END
