-
Notifications
You must be signed in to change notification settings - Fork 82
/
ex01.py
44 lines (31 loc) · 998 Bytes
/
ex01.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
from skfem import *
from skfem.helpers import dot, grad
# # enable additional mesh validity checks, sacrificing performance
# import logging
# logging.basicConfig(format='%(levelname)s %(asctime)s %(name)s %(message)s')
# logging.getLogger('skfem').setLevel(logging.DEBUG)
# create the mesh
m = MeshTri().refined(6)
# or, with your own points and cells:
# m = MeshTri(points, cells)
# or, load from file
# m = MeshTri.load("mesh.msh")
e = ElementTriP1()
basis = Basis(m, e)
@BilinearForm
def laplace(u, v, _):
return dot(grad(u), grad(v))
@LinearForm
def rhs(v, _):
return 1.0 * v
A = laplace.assemble(basis)
b = rhs.assemble(basis)
# enforce Dirichlet boundary conditions
A, b = enforce(A, b, D=m.boundary_nodes())
# solve -- can be anything that takes a sparse matrix and a right-hand side
x = solve(A, b)
def visualize():
from skfem.visuals.matplotlib import plot
return plot(m, x, shading='gouraud', colorbar=True)
if __name__ == "__main__":
visualize().show()