-
Notifications
You must be signed in to change notification settings - Fork 4
/
flux_test_2D.py
98 lines (70 loc) · 1.98 KB
/
flux_test_2D.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
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
from pylab import *
from fenics import *
from matrices import plot_matrix
n = 2
mesh = UnitSquareMesh(n,n)
#mesh = Mesh('meshes/unit_square_mesh.xml')
#mesh = Mesh('meshes/triangle.xml')
# refine mesh :
origin = Point(0.0,0.5,0.0)
#for i in range(1,8):
# cell_markers = CellFunction("bool", mesh)
# cell_markers.set_all(False)
# for cell in cells(mesh):
# p = cell.midpoint()
# if p.distance(origin) < 1.0/i:
# cell_markers[cell] = True
# mesh = refine(mesh, cell_markers)
## refine mesh :
#for i in range(1,1):
# cell_markers = CellFunction("bool", mesh)
# cell_markers.set_all(False)
# for cell in cells(mesh):
# p = cell.midpoint()
# if p.y() <= 1.0/i:
# cell_markers[cell] = True
# mesh = refine(mesh, cell_markers)
Q = FunctionSpace(mesh, 'CG', 1)
N = FacetNormal(mesh)
A = FacetArea(mesh)
w = TrialFunction(Q)
v = TestFunction(Q)
u = Function(Q, name='u')
f = Constant(1.0)
a = inner(grad(w), grad(v)) * dx
l = f * v * dx
def left(x, on_boundary):
return x[0] == 0 and on_boundary
def right(x, on_boundary):
return x[0] == 1 and on_boundary
def top(x, on_boundary):
return x[1] == 1 and on_boundary
def bottom(x, on_boundary):
return x[1] == 0 and on_boundary
bcl = DirichletBC(Q, 0.0, left)
bcr = DirichletBC(Q, 0.0, right)
bct = DirichletBC(Q, 0.0, top)
bcb = DirichletBC(Q, 0.0, bottom)
bc = [bcl]
solve(a == l, u, bc)
File('output/u.pvd') << u
uv = u.vector().array()
s = assemble(dot(grad(u), N) * v * ds)
M = assemble(w * v * dx)
b = assemble(l)
K = assemble(a)
t = (np.dot(K.array(),uv) - b.array())
q = Function(Q, name='q')
q.vector().set_local(t)
q.vector().apply('insert')
File('output/q.pvd') << q
fx = Function(Q, name='fx')
solve(M, fx.vector(), s)
File('output/fx.pvd') << fx
fig = figure(figsize=(10,5))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
plot_matrix(M.array(), ax1, r'mass matrix $M$', continuous=True)
plot_matrix(K.array(), ax2, r'stiffness matrix $K$', continuous=False)
tight_layout()
show()