标签:lin 质量 rip ret 图片 main des nod form
""" Author: kinnala Solve the problem -?2u = 1 with zero boundary conditions on a unit square. """
from skfem import *
调入 skfem 宏包
m = MeshTri()
m.refine(4)
三角剖分网格(MeshTri),加密(refine) 4 次
e = ElementTriP1()
basis = InteriorBasis(m, e)
ElememtTriP1: 三角形线性有限元 P1
InteriorBasis: 节点基函数
@bilinear_form def laplace(u, du, v, dv, w): return du[0]*dv[0] + du[1]*dv[1]
调用双线性形式模块 @bilinear_form
定义 laplace 函数
@linear_form def load(v, dv, w): return 1.0*v
调用线性泛函模块@linear_form
定义 load 函数
A = asm(laplace, basis)
b = asm(load, basis)
组装刚度矩阵 A
组装质量向量 b
I = m.interior_nodes()
取出内部网格 I
x = 0*b
x[I] = solve(*condense(A, b, I=I))
求解方程 Ax=b
if __name__ == "__main__": m.plot3(x) m.show()
画出解的图象:
标签:lin 质量 rip ret 图片 main des nod form
原文地址:https://www.cnblogs.com/wangshixi12/p/9445858.html