fealpy 介绍

本文程序部分使用了湘潭大学魏老师的有限元分析python程序包，下面给出地址。FEALPy 帮助与安装 · 魏华祎的个人网站​www.weihuayi.cnweihuayi/fealpy​github.com

一般椭圆型偏微分方程的求解

对如下椭圆形偏微分方程构造一个有真解模型的例子，用

次元进行求解，画出误差阶，并输出误差表格.

求解区域

,

取真解为

. 则有

定义 PDE

import numpy as np

from fealpy.decorator import cartesian, barycentric

class elliptic_example1:

"""

The class which define the pde.

@author: chtld

@date: 20200730

"""

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:20200730

"""

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:20200730

"""

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:20200730

"""

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:20200730

"""

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:20200730

"""

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}