FEniCS学习笔记03：求解泊松方程数学原理及程序实现

Zhang Yiyi

已于 2024-03-21 18:28:42 修改

阅读量653

点赞数 8

文章标签：学习笔记 python 开源软件

于 2024-03-13 20:32:41 首次发布

本文链接：https://blog.csdn.net/huahua_xiao/article/details/136647525

版权

内容来自官网tutorial：https://jsdokken.com/dolfinx-tutorial/chapter1/fundamentals_code.html#generating-simple-meshes

第一次运行可能会报错，解决方法见：https://blog.csdn.net/huahua_xiao/article/details/136690906?spm=1001.2014.3001.5502

1.有限元求解泊松方程数学原理

泊松方程基本形式为：

$-\bigtriangledown^2u(x)=f(x), x\in \Omega$

边界条件：

用有限元方法求解上式：将求解域离散成单元，寻找一近似解u严格满足边界条件，近似满足域内方程，再使用加权残量方法使误差在单元上最小。转换成数学表达就是对泊松方程两边同时乘以检验函数（test function）v，然后分别积分，保证等式依然成立：

对等式左边的二阶微分项进行分部积分，并引入检验函数v在边界上为0的设置条件，化简后可以得到泊松方程的弱形式（由于对方程做了降阶处理，导数连续性要求下降，所以方程限制更弱了）：

如此，我们就可以在规则的离散单元上（如三角形、四边形、四面体、六面体）用多项式近似解去拟合原本形式复杂的真解。

2.程序实现

step1:生成8*8的单位网格[0,1]*[1,0]

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)

注意：其中MPI.COMM_WORLD是并行计算时会用到的参数，表示创建一个网格，其数据分布在我们想要使用的处理器数量上。例如想在2个处理器上并行计算，可使用代码“ mpirun -n 2 python3 t1.py”来运行该程序。另外，可以使用MPI.COMM_SELF表示在每个处理器上分别创建网格，适用于不同参数设置的小型问题计算。

step2:# 定义函数空间（指定网格类型）

V = FunctionSpace(domain, ("Lagrange", 1))

表示拉格朗日1型网格，即三角形网格，节点在三个顶点上。

step3：设置dirichlet边界条件

tdim = domain.topology.dim
# 这行代码将变量dim设置为域（domain）的拓扑维度。拓扑维度表示域的几何维度，例如，一个二维域的拓扑维度为2，一个三维域的拓扑维度为3
fdim = tdim - 1
# 这行代码将变量fdim设置为域的拓扑维度减去1。这通常在有限元分析中用于计算边界面（facet）的维度
domain.topology.create_connectivity(fdim, tdim)
# 这行代码创建了域上的连接信息，以便在后续步骤中可以找到边界面和顶点之间的关联关系。fdim是边界面的维度，tdim是域的拓扑维度
boundary_facets = mesh.exterior_facet_indices(domain.topology)
# 这行代码找到域的外部边界面的索引。mesh是定义域的网格对象，domain.topology是先前创建的域的拓扑对象

boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
# 这行代码使用fem.locate_dofs_topological函数找到边界面上的自由度。V是定义有限元空间的对象，fdim是边界面的维度，boundary_facets是边界面的索引。
bc = fem.dirichletbc(uD, boundary_dofs)
# 这行代码创建一个边界条件对象bc，指定边界上的自由度以及应用于这些自由度上的边界值uD。fem.dirichletbc函数用于创建强制边界条件

本节代码注释来源于ChatGPT3.5

step4：设置试函数和检验函数、源项、变分形式

# 设置试函数和检验函数
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# 设置源项
f = fem.Constant(domain, default_scalar_type(-6))
# 设置变分形式
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx

其中，ulf.dot函数表示向量点积。若是计算矩阵内积，需使用ulf.inner。对于向量来说，两函数等价。

step5：计算方程结果和数值解误差

数值误差计算了两类。L2_error可以理解为误差（uh-uex）向量的2范数。error_L2是考虑到如果采用了多个处理器计算结果，那么需要用MPI.SUM将它们收集到一个处理器上。如果只使用了一个处理器计算，那么L2_error和error_L2等价。error_max可以理解为误差向量的无穷范数。

# 求解线性系统
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# 计算误差
V2 = fem.FunctionSpace(domain, ("Lagrange", 2))
uex = fem.Function(V2)
uex.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
L2_error = fem.form(ufl.inner(uh - uex, uh - uex) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = numpy.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))
error_max = numpy.max(numpy.abs(uD.x.array-uh.x.array))
# Only print the error on one process
if domain.comm.rank == 0:
    print(f"Error_L2 : {error_L2:.2e}")
    print(f"Error_max : {error_max:.2e}")

step6：结果后处理

官网上pyvista.start_xvfb()这行代码可能会报错，提示你安装xvfb，有图形界面的linux系统可以跳过这行代码，不影响可视化结果。

# 在网格上绘制计算结果
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()
# 绘制3D结果
warped = u_grid.warp_by_scalar()
plotter2 = pyvista.Plotter()
plotter2.add_mesh(warped, show_edges=True, show_scalar_bar=True)
if not pyvista.OFF_SCREEN:
    plotter2.show()
# 导出结果文件results
results_folder = Path("results")
results_folder.mkdir(exist_ok=True, parents=True)
filename = results_folder / "fundamentals"
with io.VTXWriter(domain.comm, filename.with_suffix(".bp"), [uh]) as vtx:
    vtx.write(0.0)
with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)

3.完整代码

# 引入必要模块
from mpi4py import MPI
from dolfinx import mesh
from dolfinx.fem import FunctionSpace
from dolfinx import fem
import numpy
import ufl
from dolfinx import default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import pyvista
from dolfinx import plot
from dolfinx import io
from pathlib import Path

# 生成8*8的单位网格[0,1]*[1,0]
domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
# 定义函数空间
V = FunctionSpace(domain, ("Lagrange", 1))
# 设置边界条件 uD=1+x^2+2y^2
uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
# 设置dirichlet边界条件
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)

boundary_dofs = fem.locate_dofs_topological(V, fdim, boundary_facets)
bc = fem.dirichletbc(uD, boundary_dofs)
# 设置试函数和检验函数
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
# 设置源项
f = fem.Constant(domain, default_scalar_type(-6))
# 设置变分形式
a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
L = f * v * ufl.dx
# 求解线性系统
problem = LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
# 计算误差
V2 = fem.FunctionSpace(domain, ("Lagrange", 2))
uex = fem.Function(V2)
uex.interpolate(lambda x: 1 + x[0]**2 + 2 * x[1]**2)
L2_error = fem.form(ufl.inner(uh - uex, uh - uex) * ufl.dx)
error_local = fem.assemble_scalar(L2_error)
error_L2 = numpy.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))
error_max = numpy.max(numpy.abs(uD.x.array-uh.x.array))
# Only print the error on one process
if domain.comm.rank == 0:
    print(f"Error_L2 : {error_L2:.2e}")
    print(f"Error_max : {error_max:.2e}")
# 画网格
print(pyvista.global_theme.jupyter_backend)
topology, cell_types, geometry = plot.vtk_mesh(domain, tdim)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("fundamentals_mesh.png")
# 在网格上绘制计算结果
u_topology, u_cell_types, u_geometry = plot.vtk_mesh(V)
u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["u"] = uh.x.array.real
u_grid.set_active_scalars("u")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()
# 绘制3D结果
warped = u_grid.warp_by_scalar()
plotter2 = pyvista.Plotter()
plotter2.add_mesh(warped, show_edges=True, show_scalar_bar=True)
if not pyvista.OFF_SCREEN:
    plotter2.show()
# 导出结果文件results
results_folder = Path("results")
results_folder.mkdir(exist_ok=True, parents=True)
filename = results_folder / "fundamentals"
with io.VTXWriter(domain.comm, filename.with_suffix(".bp"), [uh]) as vtx:
    vtx.write(0.0)
with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(uh)