/
symmetry.py
54 lines (40 loc) · 1.93 KB
/
symmetry.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
from dolfin import *
def eigenvalues(V, bcs, tol, symmetry):
u = TrialFunction(V)
v = TestFunction(V)
dummy = inner(Constant((1, 0)), v)*dx
# Assemble matrix for LHS
a = inner(curl(u), curl(v))*dx
asm = SystemAssembler(a, dummy, bcs)
A = PETScMatrix(); asm.assemble(A)
# Assemble matrix for RHS
b = inner(u, v)*dx
asm = SystemAssembler(b, dummy, bcs)
vec = PETScVector(); asm.assemble(vec)
B = PETScMatrix(); asm.assemble(B)
#assemble_system(a, dummy, bcs, A_tensor=A, b_tensor=vec)
#assemble_system(b, dummy, bcs, A_tensor=B, b_tensor=vec)
if symmetry: [bc.zero_columns(B, vec) for bc in bcs]
[bc.zero(B) for bc in bcs]
return bool(True*A.mat().isSymmetric(tol)), bool(True*B.mat().isSymmetric(tol))
def compare_eigenvalues(mesh, tol, symmetry):
nedelec_V = FunctionSpace(mesh, "N1curl", 1)
nedelec_bcs = [DirichletBC(nedelec_V, Constant((0.0, 0.0)), DomainBoundary())]
nedelec_eig = eigenvalues(nedelec_V, nedelec_bcs, tol, symmetry)
lagrange_V = VectorFunctionSpace(mesh, "Lagrange", 1)
lagrange_bcs = [DirichletBC(lagrange_V.sub(1), 0, "near(x[0], 0) || near(x[0], pi)"),
DirichletBC(lagrange_V.sub(0), 0, "near(x[1], 0) || near(x[1], pi)")]
lagrange_eig = eigenvalues(lagrange_V, lagrange_bcs, tol, symmetry)
ans = ', '.join(['Ned symmetry A:%s, B:%s' % nedelec_eig,
'Lag symmetry A:%s, B:%s' % lagrange_eig])
return ans
# -----------------------------------------------------------------------------
symmetry = False
mesh = RectangleMesh(Point(0, 0), Point(pi, pi), 40, 40)
print("\ndiagonal mesh")
for tol in [10.**-i for i in range(12)]:
print '%.2E'% tol, compare_eigenvalues(mesh, tol, symmetry)
mesh = RectangleMesh(Point(0, 0), Point(pi, pi), 40, 40, "crossed")
print("\ncrossed mesh")
for tol in [10.**-i for i in range(12)]:
print '%.2E' % tol, compare_eigenvalues(mesh, tol, symmetry)