(六)、Fealpy 组装刚度（质量）矩阵和载荷向量

给定一个单纯形网格 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 个基函数：
ϕτ=[ϕ0,ϕ1,⋯,ϕldof−1].\phi_{\tau}=[\phi_0,\phi_1,\cdots, \phi_{ldof-1}].ϕτ=[ϕ0,ϕ1,⋯,ϕldof−1].
此时 stiff matrix 为
A=∫Ω(∇ϕ)T∇ϕdx.A=\int_{\Omega} (\nabla \phi)^T\nabla\phi\,\mathrm{d}\bm{x}.A=∫Ω(∇ϕ)T∇ϕdx.
注意 AAA 是一个稀疏矩阵。 Fealpy 的想法是先在局部组装刚度矩阵，再以某种方式拼接在一起得到全局的刚度矩阵。实际计算中， Fealpy 采用数值积分的方式来实现：

# 刚度矩阵的计算
# 假设已经定义好 mesh 和 space
qf = mesh.integrator(q, 'cell')  # 获得第q个积分公式
bcs, ws = qf.get_quadrature_points_and_weights()  # 拿到相应的重心坐标和权重
cellmeasure = mesh.entity_measure('cell')  # 得到单元的面积
gphi = space.grad_basis(bcs)  # 计算基函数在重心坐标处的梯度值(NQ, NC, ldof, GD)
# NC: 单元个数 ldof:局部基函数个数 # GD: 几何维数 NQ:积分点个数
# 组装刚度矩阵 A.shape==(NC, ldof, ldof) <= (ldof, gd) * (gd, ldof)
A = np.einsum('i, ijkl, ijml, j -> jkl', ws, gphi, gphi, cellmeasure)

我们来说明一下上述代码中的最后一行：‘i’ 代表权重，‘l’代表gd, 由einsum知，这是对权重和gd求和，正好是我们要得到的式子。
有了单元刚度矩阵，我们便可以计算总体刚度矩阵。

import scipy.sparse import csr_matrixgdof = space.number_of_global_dof()  # 全局自由度的个数
# (NC, ldof) cell2dof[i,j] 存储第i个单元上局部第j个自由度的全局编号
cell2dof = space.cell_to_dof()# (NC, ldof) -> (Nc, ldof,1) -> (NC, ldof, ldof)
I = np.broadcast_to(cell2dof[:, :, None], shape=A.shape)
# (NC, ldof) -> (NC, 1, ldof) -> (NC, ldof, ldof)
J = np.broadcast_to(cell2dof[:, None, :], shape=A.shape)A = csr_matrix((A.flat, (I.flat, J.flat)), shape=(gdodf, gdof))

类似地，单元载荷向量的组装便很简单了：

实现的代码如下：

import numpy as np
from fealpy.mesh import MeshFactory
from fealpy.functionspace import LagrangeFiniteElementSpace
from matplotlib import pyplot as plt
from scipy.sparse import csr_matrix
mf = MeshFactory()
mesh = mf.boxmesh2d([0, 1, 0, 1], nx=1, ny=1, meshtype='tri')fig = plt.figure()
axes = fig.gca()
mesh.add_plot(axes)space = LagrangeFiniteElementSpace(mesh, p=1)
gdof = space.number_of_global_dofs()
print('gdof:', gdof)cell2dof = space.cell_to_dof()
print('cell2dof:', cell2dof)
cell = mesh.entity('cell')
print('cell:', cell)mesh.find_node(axes, showindex=True, fontsize=24)
mesh.find_cell(axes, showindex=True, fontsize=24)cellmeasure = mesh.entity_measure('cell')
print('cellmeasure:', cellmeasure)qf = mesh.integrator(1, 'cell')  # 积分公式
bcs, ws = qf.get_quadrature_points_and_weights()
print('bcs:', bcs, 'ws:', ws)  # 重心坐标与权重
gphi = space.grad_basis(bcs)  # [1, 2, 3, 2]  q * NC * gphi * dim 积分点*单元*基函数*导数分量(维数)A = np.einsum('i, ijkl, ijml, j-> jkm', ws, gphi, gphi, cellmeasure)  # (2, 3, 3) 单元 基函数 基函数I = np.broadcast_to(cell2dof[:, :, None], shape=A.shape)
J = np.broadcast_to(cell2dof[:, None, :], shape=A.shape)
# A, I, J  # (2, 3, 3)
for i, j, val in zip(I.flat, J.flat, A.flat):print('(', i, ',', j, '):', val)
A = csr_matrix((A.flat, (I.flat, J.flat)), shape=(gdof, gdof))  # (gdof, gdof)
A.toarray()# mass matrix
phi = space.basis(bcs)  # (NQ, 1, ) 1:NC 节省内存
M = np.einsum('i, ijk, ijm, j->jkm', ws, phi, phi, cellmeasure)
M = csr_matrix((M.flat, (I.flat, J.flat)), shape=(gdof, gdof))  # (gdof, gdof)
M.toarray()# 对流扩散方程
# 对流项组装
b = np.array([1, 1], dtype=np.float64)
B = np.einsum('i, q, ijkq, ijm, j-> jkm', ws, b, gphi, phi, cellmeasure)
B = csr_matrix((B.flat, (I.flat, J.flat)), shape=(gdof, gdof))  # (gdof, gdof)
B.toarray()
# 扩散项组装
D = np.array([[10.0, -1.0], [-1.0, 2.0]], dtype=np.float64)
B = np.einsum('i,sq, ijkq, ijls,j->jkl', ws, D, gphi, gphi, cellmeasure)
B = csr_matrix((B.flat, (I.flat, J.flat)), shape=(gdof, gdof))  # (gdof, gdof)
print('B', B.toarray())B = space.stiff_matrix(c=D)
print(B.toarray())
# 组装向量
from fealpy.decorator import cartesian
@cartesian
def f(p):x = p[..., 0]y = p[..., 1]return np.exp(x**2 + y**2)
# print(f.coordtype)
qf= mesh.integrator(3, 'cell')
bcs, ws = qf.get_quadrature_point_and_weight()
phi = space.basis(bcs)  # (NQ, 1, 3) 基函数
ps = mesh.bc_to_point(bcs)  # (6, 2, 2) 积分点 单元 分量
val = f(ps)
bb = np.einsum('i, ij, ijk, j -> jk', ws, val, phi, cellmeasure)
F = np.zeros(gdof, dtype=np.float64)
np.add.at(F, cell2dof, bb)plt.show()

笔者感悟：我在魏老师上课代码中加了一个扩散项的组装，经过与fealpy中space.stiff_matrix(c=diffusion_coefficient) 比对，结果是对的。我花了三小节的时间来整理、处理矩阵的组装，是因为我觉得这是fealpy 的灵魂所在：包括cell_to_dof， csr_matrix, np.einsum 的引入等等。学习felapy, 我们要学到其背后的原理与编程的想法，这对我我们解自己专业的问题至关重要。我这么说并不意味着用fealpy解PDE都要自己组装。事实上，魏老师写了很多PDE的接口，下一节我们就来介绍这些接口，和展示一些数值实验。

(六)、Fealpy 组装刚度（质量）矩阵和载荷向量相关推荐

matlab 单元质量矩阵,这样组装平面桁架的质量矩阵对不对呢？
%单元刚度矩阵 k1=[0.7375000E+08 0.0000000E+00 -0.7375000E+08 0.0000000E+00; 0.0000000E+00 0.0000000E+00 ...
ANSYS中关于质量矩阵刚度矩阵的提取【1】
ANSYS就想对模态结果进行二次处理,求的所需要的刚度矩阵及质量矩阵,方便对结果进行进一步的分析研究关于SIMWE仿真科技论坛.一个是 ANSYS中提取模态分析结果中得到整体结构的质量矩阵和刚度矩阵 ...
[振动力学] 使用能量法求质量矩阵的时候需要注意刚体运动分解
使用能量法求质量矩阵的时候,对于杆的动能,我有两种写法 ①求出的质量矩阵没有耦合项,而②有耦合项显然正确答案只有一个,那就是② 在①中我是把运动分解的基点取在了杆和圆柱体的连接点处,实际上并不可以把 ...
六西格玛，为质量人的职业发展保驾护航
常常会有做质量的人来问优思学院:在质量界混,应学什么好? 大部分时候,我都会提议他們:有空就要学六西格玛吧! 他会回应:六西格玛感觉太难懂了,学了也没什么用啊,哪有那么多公司搞这个啊? 的确讲得好像有 ...
用矩阵来运算向量与点的平移
用4*4的矩阵来描述向量与点: 1.为什么要用4*4的矩阵,而不是3*3的矩阵呢? 因为在3D世界中,描述一个点至少需要3个维度,如果使用3个维度来描述向量或者点, 那么点与向量就没法区别对待,但是点 ...
numpy向量转换为矩阵_Numpy之将矩阵拉成向量的实例
Numpy之将矩阵拉成向量的实例废话不多说,直接上代码吧! # 矩阵操作 # 将矩阵拉成向量 import numpy as np x = np.arange(10).reshape(2,5) pr ...
【机器学习线性代数】02 初识矩阵：让向量动起来
1.矩阵?一排向量,一堆数介绍完了向量,这一节我们开始介绍矩阵.对于矩阵而言,最直观的描述就是一个m×nm\times nm×n大小的数字方阵,他可以看作是nnn个mmm维列向量从左到右并排摆放,也 ...
matlab中矩阵怎么敲_Ansys刚度(质量、阻尼)矩阵的提取(part 1)
在做数值计算分析中要验证正确性,刚度.质量和阻尼矩阵是有限元分析最终最重要的3个矩阵.本文只针对刚度矩阵,不过我认为质量和阻尼矩阵应该是同理. 王新敏老师的<ansys工程结构数值分析>中 ...
高通量数据中批次效应的鉴定和处理（六）- 直接校正表达矩阵
生物信息学习的正确姿势 NGS系列文章包括NGS基础.转录组分析 (Nature重磅综述|关于RNA-seq你想知道的全在这).ChIP-seq分析 (ChIP-seq基本分析流程).单细胞测序分析 ...

(六)、Fealpy 组装刚度（质量）矩阵和载荷向量

(六)、Fealpy 组装刚度（质量）矩阵和载荷向量相关推荐

最新文章

热门文章