waitingwinter的博客

Fealpy 是湘潭大学魏华祎老师开发的用来求解有限元的库，其中设置了多种网格（三角网格，矩形网格，多边形网格），多种有限元空间（Lagrange有限元空间，缩放多项式空间等），自适应方法等。本系列是笔者听完魏老师在天元数学东北中心开设的暑期短课，以及自身对有限元程序的需求之下写的一些回忆总结性片段，一来用于记录，二来方便与各位朋友交流。相关的参考链接有：Fealpy的安装教程(魏老师主页)Fealpy2020暑期课程...

Fealpy 库中实现了一个类 mesh 用于网格生成, 我们先介绍规则的网格（默认大家已经完成 fealpy 的安装， 见前言）首先引入一些库：import numpy as npimport matplotlib.pyplot as pltfrom fealpy.mesh import QuadrangleMeshfrom fealpy.mesh import HalfEdgeMesh2d, PolygonMeshfrom fealpy.mesh import TriangleMeshfr

我们先给出三角形单元的中心坐标公式对应的代码为：import numpy as npfrom fealpy.mesh import MeshFactorymf = MeshFactory()box = [0, 1, 0, 1]mesh = mf.boxmesh2d(box, nx=1, ny=1, meshtype='tri')NC = mesh.number_of_cells()node = mesh.entity('node')cell = mesh.entity('cell')

这一讲，我们来介绍 Fealpy 中的有限元空间。 Fealpy 中实现了几种常用的有限元空间， 比如拉格朗日有限元空间， 缩放多项式空间等。首先我们从 Functionspace 中引入 LagrangeFiniteElementSpace, 当然也要引入之前异界介绍的网格工厂（Meshfactory）:import numpy as npfrom fealpy.mesh import MeshFactoryfrom fealpy.functionspace import LagrangeFini

在介绍具体的细节之前， 我们先来介绍 numpy 的两个函数 以及稀疏矩阵的存储np.add.atindices = np.array([0, 0, 1, 5, 3], dtype=np.int_)a = np.zeros(6, dtype=np.float64)data = np.array([4.5, 3.5, 2.0, 1.0, 1.5])np.add.at(a, indices, data)print('a:' , a)add at 的命令是将data放在指定indices的位置上。

这一节， 我们来介绍一下有限元相关的知识。我们以 Poisson 方程为例：注意我们的测试函数 uuu 在 ΓD\Gamma_DΓD​ 上是 0。于是我们可以给出离散的代数系统接下来我们来分析边界条件， 若只有 Dirichlet 边界条件， 则有离散代数系统中的形式知 R=0,bN=bR=0,Au=b.R=0, b_N=b_R=0, Au=b.R=0,bN​=bR​=0,Au=b. 将 uuu 进行向量分解u=uI(内部自由度)+uD(边界自由度)，u=u_I(内部自由度)+u_D(边界自

给定一个单纯形网格 T\mathcal{T}T, 其有 NNNNNN 个节点， NCNCNC 个单元。 定义在 T\mathcal{T}T 的分片 ppp 次连续有限元空间 VhV_hVh​ 有 gdofgdofgdof 个基函数， 其组成的函数的行向量为：ϕ=[ϕ0,ϕ1,⋯ ,ϕgdof−1],\phi=[\phi_0,\phi_1,\cdots,\phi_{gdof-1}],ϕ=[ϕ0​,ϕ1​,⋯,ϕgdof−1​],限制在每个网格单元 τ\tauτ 上， 共有 ldofldofldof 个基函

这一节， 我们来介绍一下 Fealpy 中的PDE接口。 打开 fealpy/fealpy/pde, 你可以看到大量的pde接口：fealpy/example 中有很多求解的例子， 基本涵盖了 pde solver 等的东西，如下图：因为作者的文档写的很清楚了， 故我这里只举一个例子来简单介绍一下：混合边界条件的poisson 方程求解混合边界条件的方程如下：# 构造一个二维有真解的例子，用 p 次有限元来求解， 画出误差阶， 输出误差表格# -\Delta u + 3 u = f, in