【我已经发了同名西瓜视频，本文仅是视频中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_EPSbc = 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)*dxL = f*v*dx - u0*dot(n,tau)*ds

# 求解w = Function(W)solve(a == L, w, bc)(sigma, u) = w.split()# 绘图plot(sigma)plot(u)

# 将解保存为VTK格式vtkfile = File("mixed_poisson/sigma.pvd")vtkfile << sigmavtkfile = File("mixed_poisson/u.pvd")vtkfile << u

# 绘制矢量场图(如果是矢量场的话)# 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