How I can plot eigenvectors as a function of two dimension (x,y))?

Linear Algebra
1 
from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
from mshr import *
#################################
V_con = 20
nm = 10 ** -9
a_d = 10 ** - 9

rc = 30 * nm / a_d
zc = 50 * nm/  a_d
H  = 10 * nm / a_d
r  = 8 * nm / a_d
n  = 15

domain_w = Rectangle(Point(0,-H), Point(r,2 * H))
domain_b = Rectangle(Point(0,-zc), Point(rc,zc))
domain_0 = domain_b - domain_w

    # mark domain
domain_b.set_subdomain(1, domain_w)
domain_b.set_subdomain(2, domain_0)
    #creat the mesh
mesh = generate_mesh(domain_b,n)
marker_subdomain = MeshFunction('size_t', mesh,2, mesh.domains())

V =  FunctionSpace(mesh, "CG", 1)

bc = DirichletBC(V, Constant(0), 'on_boundary')

dx = Measure("dx")[marker_subdomain]
u = TrialFunction(V)
v = TestFunction(V)
v_con = Constant(V_con)
r = Expression("x[0]" , degree=2)


a = r * ( inner(nabla_grad(u), nabla_grad(v))) * dx(1) + r * ( inner(nabla_grad(u), nabla_grad(v)) + v_con * u * v ) * dx(2)
b = r * u * v * dx(1) + r * u * v * dx(2)

   # Assemble matrices
H = PETScMatrix()
M = PETScMatrix()

assemble(a, tensor=H)
assemble(b, tensor=M)

    # Create eigensolver for generalized EVP
eigensolver = SLEPcEigenSolver(H, M)
eigensolver.parameters["spectrum"] = "smallest magnitude"
eigensolver.parameters["solver"] = "lapack"

    # Compute generalized eigenvalue decomposition
print("Computing eigenvalue decomposition. This may take a while.")
eigensolver.solve()
i = 0
E0 = []
phi = []
while i < eigensolver.get_number_converged():
    r, c, rx, cx = eigensolver.get_eigenpair(i)
    E0.append(r)
    phi.append(np.array(rx))                   #
    i = i + 1


plt.xlabel('x')
plt.ylabel('y')

plot(marker_subdomain)
plot(mesh)
plt.show()

Thank you in advance for your helping.

2

See here, for example.

1 Like
3

Thank you for you answer.
I think you didn’t understand what I want exactly ?
,I want to plot eigenvector in three dimensions
Similar to image that attached below.
GS2DAnd3DPlotsExample_05

4

Hello,
define before the loop phi = Function(V) and inside the loop do

phi.vector().set_local(rx)
plot(phi, mode="warp")
1 Like
5

That’s very kind of you. Thank you.

6

from dolfin import *
import matplotlib.pyplot as plt
import numpy as np
from mshr import *
#################################
V_con = 20
nm = 10 ** -9
a_d = 10 ** - 9

rc = 30 * nm / a_d
zc = 50 * nm/ a_d
H = 10 * nm / a_d
r = 8 * nm / a_d
n = 15

domain_w = Rectangle(Point(0,-H), Point(r,2 * H))
domain_b = Rectangle(Point(0,-zc), Point(rc,zc))
domain_0 = domain_b - domain_w

# mark domain

domain_b.set_subdomain(1, domain_w)
domain_b.set_subdomain(2, domain_0)
#creat the mesh
mesh = generate_mesh(domain_b,n)
marker_subdomain = MeshFunction(‘size_t’, mesh,2, mesh.domains())

V = FunctionSpace(mesh, “CG”, 1)

bc = DirichletBC(V, Constant(0), ‘on_boundary’)

dx = Measure(“dx”)[marker_subdomain]
u = TrialFunction(V)
v = TestFunction(V)
v_con = Constant(V_con)
r = Expression(“x[0]” , degree=2)

a = r * ( inner(nabla_grad(u), nabla_grad(v))) * dx(1) + r * ( inner(nabla_grad(u), nabla_grad(v)) + v_con * u * v ) * dx(2)
b = r * u * v * dx(1) + r * u * v * dx(2)

Assemble matrices

H = PETScMatrix()
M = PETScMatrix()

assemble(a, tensor=H)
assemble(b, tensor=M)

# Create eigensolver for generalized EVP

eigensolver = SLEPcEigenSolver(H, M)
eigensolver.parameters[“spectrum”] = “smallest magnitude”
eigensolver.parameters[“solver”] = “lapack”

# Compute generalized eigenvalue decomposition

print(“Computing eigenvalue decomposition. This may take a while.”)
eigensolver.solve()
i = 0
E0 = []
phi = Function(V)
r, c, rx, cx = eigensolver.get_eigenpair(0)
phi.vector().set_local(rx)
plot(phi, mode=“warp”)

plt.xlabel(‘x’)
plt.ylabel(‘y’)

plt.show()

Dear Sir/ Madam
Thank you for your answer. from your answer I rise a one question.
the plotted function is smaller versus figure, How I can change it ?
I made a lot of searches but unfortunately I didn’t find anything.

7

Please encapsulate your code in ``` for proper formatting, and reduce the size of your code by removing all lines of the code that is after the error message to make it more readable for others.
See Read before posting: How do I get my question answered? on how to post a minimal working example.

1 Like
8

Thank for your reply.
I don’t have a problem or error in my code.
I just want to make the figure very clearly.

Figure_1-4
Figure_1-41533×811 54.9 KB

9

Then why did you post all the code, it obscures the question. And whenever you post code it should be formatted. To run your code, anyone has to do a lot of manual formatting.

I ran you script using the latest development docker image, replacing plt.show() with plt.savefig("test.png") and got the following:
test
Anyhow, 3D visualization is often easiest in Paraview, where there are functionality such as warp as vector and warp as scalar.

1 Like
10

Thank you very much.

11

vtkplotter can visualize 2D and 3D solutions easily e.g.:

# ...
from vtkplotter.dolfin import plot

plot(phi, warpZfactor=150, cmap='gist_earth_r', lighting='plastic')

image
image606×532 76.8 KB

3 Likes
12

Thank you very much for your help in this matter.

1 Like