fealpy 介绍
本文程序部分使用了湘潭大学魏老师的有限元分析python程序包,下面给出地址。FEALPy 帮助与安装 · 魏华祎的个人网站weihuayi/fealpy
一般椭圆型偏微分方程的求解
对如下椭圆形偏微分方程构造一个有真解模型的例子,用
次元进行求解,画出误差阶,并输出误差表格.
求解区域
,
取真解为
. 则有
定义 PDE
import numpy as np
from fealpy.decorator import cartesian, barycentric
class elliptic_example1:
"""
The class which define the pde.
@author: chtld
@date: 2030
"""
def __init__(self):
pass
@property
def domain(self):
return np.array([0, 1, 0, 1])
@cartesian
def source(self, p):
x = p[..., 0]
y = p[..., 1]
pi = np.pi
val = (12 * pi ** 2 + 1 + x ** 2 + y ** 2) * np.cos(pi * x) * np.cos(pi * y) \
+ 2 * pi ** 2 * np.sin(pi * x) * np.sin(pi * y) \
- pi * np.sin(pi * x) * np.cos(pi * y) - pi * np.cos(pi * x) * np.sin(pi * y)
return val
@cartesian
def exact_solution(self, p):
x = p[..., 0]
y = p[..., 1]
pi = np.pi
val = np.cos(pi * x) * np.cos(pi * y)
return val
@cartesian
def gradient(self, p):
x = p[..., 0]
y = p[..., 1]
pi = np.pi
val = np.zeros(p.shape, dtype = np.float64)
val[..., 0] = -pi * np.sin(pi * x) * np.cos(pi * y)
val[..., 1] = -pi * np.cos(pi * x) * np.sin(pi * y)
return val
@cartesian
def flux(self, p):
return -self.gradient(p)
@cartesian
def dirichlet(self, p):
return self.exact_solution(p)
@cartesian
def diffusion_coefficient(self, p):
x = p[..., 0]
y = p[..., 1]
return np.array([[10.0 * x * y / x / y, -1.0 * x * y / x / y], [-1.0 * x * y / x / y, 2.0 * x * y / x / y]], dtype=np.float64)
@cartesian
def convection_coefficient(self, p):
x = p[..., 0]
y = p[..., 1]
return np.array([1.0 * x * y / x / y, 1.0 * x * y / x / y], dtype=np.float64)
@cartesian
def reaction_coefficient(self, p):
x = p[..., 0]
y = p[..., 1]
return 1 + x**2 + y**2
弱形式
代数系统
其中
一次性组装矩阵和向量
def assemble_matrix_and_vector(pde, mesh, femspace, q):
""" 组装矩阵和向量
input:
@pde: 偏微分方程对象
@mesh: 网格对象
@femspace: 有限元空间
@q: gauss积分公式
output:
matrix and vector
@author: chtld
@date:2030
"""
quadrature_rule = mesh.integrator(q, 'cell')
barycentric_points, weights = quadrature_rule.get_quadrature_points_and_weights()
cell_measure = mesh.entity_measure('cell')
gphi = femspace.grad_basis(barycentric_points)
cartesian_points = mesh.bc_to_point(barycentric_points)
diffusion_coef = pde.diffusion_coefficient(cartesian_points)
diffusion_matrix = np.einsum('i, rcij, ijnr, ijmc, j -> jmn', weights, diffusion_coef, gphi, gphi, cell_measure)
gdof = femspace.number_of_global_dofs()
cell2dof = femspace.cell_to_dof()
I = np.broadcast_to(cell2dof[:, :, None], shape = diffusion_matrix.shape)
J = np.broadcast_to(cell2dof[:, None, :], shape = diffusion_matrix.shape)
diffusion_matrix = csr_matrix((diffusion_matrix.flat, (I.flat, J.flat)), shape = (gdof, gdof))
phi = femspace.basis(barycentric_points)
convection_coef = pde.convection_coefficient(cartesian_points)
convection_matrix = np.einsum('i, rij, ijnr, ijm, j -> jmn', weights, convection_coef, gphi, phi, cell_measure)
convection_matrix = csr_matrix((convection_matrix.flat, (I.flat, J.flat)), shape = (gdof, gdof))
reaction_coef = pde.reaction_coefficient(cartesian_points)
reaction_matrix = np.einsum('i, ij, ijn, ijm, j -> jmn', weights, reaction_coef, phi, phi, cell_measure)
reaction_matrix = csr_matrix((reaction_matrix.flat, (I.flat, J.flat)), shape = (gdof, gdof))
righthandside = pde.source(cartesian_points)
righthand_vector = np.einsum('i, ij, ijm, j -> jm', weights, righthandside, phi, cell_measure)
righthand = np.zeros(gdof, dtype = np.float64)
np.add.at(righthand, cell2dof, righthand_vector)
return diffusion_matrix + convection_matrix + reaction_matrix, righthand
分别组装矩阵和向量
组装扩散项
def assemble_diffusion_matrix(pde, mesh, femspace, q):
""" 组装扩散矩阵
input:
@pde: 偏微分方程对象
@mesh: 网格对象
@femspace: 有限元空间
@q: gauss积分公式
output:
matrix
@author: chtld
@date:2030
"""
quadrature_rule = mesh.integrator(q, 'cell')
barycentric_points, weights = quadrature_rule.get_quadrature_points_and_weights()
cell_measure = mesh.entity_measure('cell')
gphi = femspace.grad_basis(barycentric_points)
cartesian_points = mesh.bc_to_point(barycentric_points)
diffusion_coef = pde.diffusion_coefficient(cartesian_points)
# 计算单元矩阵
diffusion_matrix = np.einsum('i, rcij, ijnr, ijmc, j -> jmn', weights, diffusion_coef, gphi, gphi, cell_measure)
gdof = femspace.number_of_global_dofs()
cell2dof = femspace.cell_to_dof()
I = np.broadcast_to(cell2dof[:, :, None], shape = diffusion_matrix.shape)
J = np.broadcast_to(cell2dof[:, None, :], shape = diffusion_matrix.shape)
# 组装进总体矩阵
diffusion_matrix = csr_matrix((diffusion_matrix.flat, (I.flat, J.flat)), shape = (gdof, gdof))
return diffusion_matrix
组装对流项
def assemble_convection_matrix(pde, mesh, femspace, q):
""" 组装对流矩阵
input:
@pde: 偏微分方程对象
@mesh: 网格对象
@femspace: 有限元空间
@q: gauss积分公式
output:
matrix
@author: chtld
@date:2030
"""
quadrature_rule = mesh.integrator(q, 'cell')
barycentric_points, weights = quadrature_rule.get_quadrature_points_and_weights()
cell_measure = mesh.entity_measure('cell')
gphi = femspace.grad_basis(barycentric_points)
cartesian_points = mesh.bc_to_point(barycentric_points)
phi = femspace.basis(barycentric_points)
convection_coef = pde.convection_coefficient(cartesian_points)
# 形成单元矩阵
convection_matrix = np.einsum('i, rij, ijnr, ijm, j -> jmn', weights, convection_coef, gphi, phi, cell_measure)
gdof = femspace.number_of_global_dofs()
cell2dof = femspace.cell_to_dof()
I = np.broadcast_to(cell2dof[:, :, None], shape = convection_matrix.shape)
J = np.broadcast_to(cell2dof[:, None, :], shape = convection_matrix.shape)
# 组装进入总体矩阵
convection_matrix = csr_matrix((convection_matrix.flat, (I.flat, J.flat)), shape = (gdof, gdof))
return convection_matrix
组装反应项
def assemble_reaction_matrix(pde, mesh, femspace, q):
""" 组装反应矩阵
input:
@pde: 偏微分方程对象
@mesh: 网格对象
@femspace: 有限元空间
@q: gauss积分公式
output:
matrix
@author: chtld
@date:2030
"""
quadrature_rule = mesh.integrator(q, 'cell')
barycentric_points, weights = quadrature_rule.get_quadrature_points_and_weights()
cell_measure = mesh.entity_measure('cell')
cartesian_points = mesh.bc_to_point(barycentric_points)
phi = femspace.basis(barycentric_points)
reaction_coef = pde.reaction_coefficient(cartesian_points)
# 计算单元矩阵
reaction_matrix = np.einsum('i, ij, ijn, ijm, j -> jmn', weights, reaction_coef, phi, phi, cell_measure)
gdof = femspace.number_of_global_dofs()
cell2dof = femspace.cell_to_dof()
I = np.broadcast_to(cell2dof[:, :, None], shape = reaction_matrix.shape)
J = np.broadcast_to(cell2dof[:, None, :], shape = reaction_matrix.shape)
# 组装进总体矩阵
reaction_matrix = csr_matrix((reaction_matrix.flat, (I.flat, J.flat)), shape = (gdof, gdof))
return reaction_matrix
组装右端源项
def assemble_source_vector(pde, mesh, femspace, q):
""" 组装右端源向量
input:
@pde: 偏微分方程对象
@mesh: 网格对象
@femspace: 有限元空间
@q: gauss积分公式
output:
vector
@author: chtld
@date:2030
"""
quadrature_rule = mesh.integrator(q, 'cell')
barycentric_points, weights = quadrature_rule.get_quadrature_points_and_weights()
cell_measure = mesh.entity_measure('cell')
cartesian_points = mesh.bc_to_point(barycentric_points)
gdof = femspace.number_of_global_dofs()
cell2dof = femspace.cell_to_dof()
phi = femspace.basis(barycentric_points)
righthandside = pde.source(cartesian_points)
righthand_vector = np.einsum('i, ij, ijm, j -> jm', weights, righthandside, phi, cell_measure)
righthand = np.zeros(gdof, dtype = np.float64)
np.add.at(righthand, cell2dof, righthand_vector)
return righthand
有限元求解器
import sys
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import spsolve
from fealpy.mesh import MeshFactory
from fealpy.functionspace import LagrangeFiniteElementSpace
from fealpy.boundarycondition import DirichletBC
from fealpy.tools.show import showmultirate, show_error_table
def elliptic_example1_solver(pde, n, order, refine):
mf = MeshFactory()
mesh = mf.boxmesh2d(pde.domain, nx = n, ny = n, meshtype = 'tri')
number_of_dofs = np.zeros(refine, dtype = mesh.itype)
error_matrix = np.zeros((2, refine), dtype = mesh.ftype)
error_type = ['$||u - u^{h}||_{0}$', '$||\\nabla u - \\nabla u^{h}||_{0}$']
for i in range(refine):
femspace = LagrangeFiniteElementSpace(mesh, p = order)
number_of_dofs[i] = femspace.number_of_global_dofs()
uh = femspace.function()
#A, F = assemble_matrix_and_vector(pde, mesh, femspace, order * 2)
A = assemble_diffusion_matrix(pde, mesh, femspace, order * 2)
C = assemble_convection_matrix(pde, mesh, femspace, order * 2)
R = assemble_reaction_matrix(pde, mesh, femspace, order * 2)
F = assemble_source_vector(pde, mesh, femspace, order * 2)
bc = DirichletBC(femspace, pde.dirichlet)
A, F = bc.apply(A + C + R, F, uh)
uh[:] = spsolve(A, F)
error_matrix[0, i] = femspace.integralalg.L2_error(pde.exact_solution, uh.value)
error_matrix[1, i] = femspace.integralalg.L2_error(pde.gradient, uh.grad_value)
if i < refine - 1:
mesh.uniform_refine()
fig = plt.figure()
axes = fig.gca(projection = '3d')
uh.add_plot(axes, cmap = 'rainbow')
showmultirate(plt, 0, number_of_dofs, error_matrix, error_type, propsize = 20)
show_error_table(number_of_dofs, error_type, error_matrix, f = 'e', pre = '4', sep = '&', out = sys.stdout, end = '\n')
plt.show()
线性元求解
pde = elliptic_example1()
n, order, refine = 10, 1, 5
elliptic_example1_solver(pde, n, order, refine)
run boxmesh2d with time: 0.001358699999999935
\begin{table}[!htdp]
\begin{tabular}[c]{|c|c|c|c|c|c|}\hline
Dof& 121& 441& 1681& 6561&25921
\\\hline
$||u - u^{h}||_{0}$&1.3167e-02&3.3584e-03&8.4405e-04&2.1130e-04&5.2842e-05
\\\hline
Order&--&1.97&1.99&2. &2.
\\\hline
$||\nabla u - \nabla u^{h}||_{0}$&3.4685e-01&1.7421e-01&8.7203e-02&4.3614e-02&2.1808e-02
\\\hline
Order&--&0.99&1. &1. &1.
\\\hline
\end{tabular}
\end{table}
二次元求解
pde = elliptic_example1()
n, order, refine = 10, 2, 5
elliptic_example1_solver(pde, n, order, refine)
run boxmesh2d with time: 0.0008777000000002033
\begin{table}[!htdp]
\begin{tabular}[c]{|c|c|c|c|c|c|}\hline
Dof& 441& 1681& 6561& 25921&103041
\\\hline
$||u - u^{h}||_{0}$&2.3859e-04&2.9421e-05&3.6621e-06&4.5725e-07&5.7140e-08
\\\hline
Order&--&3.02&3.01&3. &3.
\\\hline
$||\nabla u - \nabla u^{h}||_{0}$&2.1603e-02&5.4056e-03&1.3513e-03&3.3779e-04&8.4446e-05
\\\hline
Order&--&2.&2.&2.&2.
\\\hline
\end{tabular}
\end{table}
三次元求解
pde = elliptic_example1()
n, order, refine = 10, 3, 5
elliptic_example1_solver(pde, n, order, refine)
run boxmesh2d with time: 0.0005549000000044657
\begin{table}[!htdp]
\begin{tabular}[c]{|c|c|c|c|c|c|}\hline
Dof& 961& 3721& 14641& 58081&231361
\\\hline
$||u - u^{h}||_{0}$&7.1793e-06&4.2430e-07&2.5829e-08&1.5947e-09&9.9091e-11
\\\hline
Order&--&4.08&4.04&4.02&4.01
\\\hline
$||\nabla u - \nabla u^{h}||_{0}$&8.5307e-04&1.0568e-04&1.3148e-05&1.6396e-06&2.0471e-07
\\\hline
Order&--&3.01&3.01&3. &3.
\\\hline
\end{tabular}
\end{table}