码迷,mamicode.com
首页 > 其他好文 > 详细

scikit-FEM-例1-求解Possion边值问题

时间:2018-08-08 22:05:13      阅读:353      评论:0      收藏:0      [点我收藏+]

标签: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()

画出解的图象:

技术分享图片

 

scikit-FEM-例1-求解Possion边值问题

标签:lin   质量   rip   ret   图片   main   des   nod   form   

原文地址:https://www.cnblogs.com/wangshixi12/p/9445858.html

(0)
(0)
   
举报
评论 一句话评论(0
登录后才能评论!
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码!