How to extract and plot the eigen-function

dolfinx
1

Hi,
I’m using SLEPC to solve an eigenvalue problem.
I want to plot the eigenfunctions with fenicsx.
I have tried from
[A simple eigenvalue solver in DOLFIN-X - #3 by bhaveshshrimali]

the part:

Blockquote
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
eigensolver.setOperators(A)
eigensolver.solve()
vr, vi = A.getVecs()
lmbda = eigensolver.getEigenpair(0, vr, vi)
u = Function(V)
u.vector.setArray(vr.array)
plt.figure(figsize=(8,8))
eig1 = plot(u, cmap=plt.cm.jet)
plot(mesh)
plt.colorbar(eig1)

but got the following error:
" TypeError: ‘module’ object is not callable"

2

The plotting functionality in dolfinx has been completely rewritten Since 2020. See for instance

https://jsdokken.com/dolfinx-tutorial/
Or

for examples

3

I am facing this error:
cells, types, x = plot.vtk_mesh(space)
AttributeError: module ‘dolfinx.plot’ has no attribute ‘vtk_mesh’

4

As you are probably using v0.6.x, see the appropriate tag: https://github.com/FEniCS/dolfinx/blob/v0.6.0/python/demo/demo_pyvista.py

5

I’m still stuck in plotting the solution. Can you guide me more on how to plot the eigenfunction it with fenicsx.

6

Without posting the error message and corresponding code that you are using, i cannot give you much further guidance.

I’ve already referenced several links that show how to plot functions in dolfinx.

Alternatively, you can save the function to a file (using XDMFFile, VTKFile, VTXWriter or FidesWriter) and use external tools such as paraview to visualize the solution

7

Hi again, Please find the code and the attached error:
(sorry don’t know how to quote it correctly, I use "Block quote )

Blockquote
from dolfinx.mesh import create_unit_square, create_box, CellType
from dolfinx import mesh, fem
from dolfinx.fem import locate_dofs_geometrical, Constant, form
import pyvista # visualizing the mesh using pyvista, an interface to the VTK toolkit.
print(pyvista.global_theme.jupyter_backend)
import dolfinx.plot as plot
from ufl import TrialFunction, TestFunction, dot, div, inner ,dx, SpatialCoordinate
import numpy as np
from petsc4py import PETSc
from slepc4py import SLEPc
import matplotlib.pyplot as plt
import math
from mpi4py import MPI
import sys, io, slepc4py, os.path
factor = math.pi*math.pi
slepc4py.init(sys.argv)
opts = PETSc.Options()
import dolfinx
print(f"DOLFINx version: {dolfinx.version} based on GIT commit: {dolfinx.git_commit_hash} of GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment")

Blockquote
#==================
import gmsh
import math
import sys
import meshio
#==================
def lshape_unstructQuad_gmsh2(element, num_elem):
print(’ The number of elements:‘, num_elem)
if (num_elem%2) != 0:
print(f"Insurt an even number of elements, the given number is {num_elem}")
exit()
num_of_points = num_elem + 1
gmsh.initialize()
lc = 1e-2
# boundaries: [-1,1]x[-1,1]
p1 = gmsh.model.geo.addPoint(-1 , -1 , 0, lc, 1)
p2 = gmsh.model.geo.addPoint(0 , -1 , 0, lc, 2)
p3 = gmsh.model.geo.addPoint(0 , 0 , 0, lc, 3)
p4 = gmsh.model.geo.addPoint(1 , 0 , 0, lc, 4)
p5 = gmsh.model.geo.addPoint(1 , 1 , 0, lc, 5)
p6 = gmsh.model.geo.addPoint(-1 , 1 , 0, lc, 6)
input(‘The boundaries of the domain are [-1,-1]x[1,1]:’)
l1 = gmsh.model.geo.addLine(p1, p2, 1)
l2 = gmsh.model.geo.addLine(p2, p3, 2)
l3 = gmsh.model.geo.addLine(p3, p4, 3)
l4 = gmsh.model.geo.addLine(p4, p5, 4)
l5 = gmsh.model.geo.addLine(p5, p6, 5)
l6 = gmsh.model.geo.addLine(p6, p1, 6)
# Adding points, lines, Create surface
gmsh.model.geo.addCurveLoop([l1, l2, l3, l4, l5, l6], 1)
gmsh.model.geo.addPlaneSurface([1], 1)
gmsh.model.addPhysicalGroup(0, [1, 2, 3, 4, 5, 6], 1)
gmsh.model.addPhysicalGroup(1, [1, 2, 3, 4, 5, 6], 2)
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.geo.mesh.setTransfiniteCurve(l1, int(num_of_points/2)+1)
gmsh.model.geo.mesh.setTransfiniteCurve(l2, int(num_of_points/2)+1)
gmsh.model.geo.mesh.setTransfiniteCurve(l3, int(num_of_points/2)+1)
gmsh.model.geo.mesh.setTransfiniteCurve(l4, int(num_of_points/2)+1)
gmsh.model.geo.mesh.setTransfiniteCurve(l5, num_of_points)
gmsh.model.geo.mesh.setTransfiniteCurve(l6, num_of_points)
# # To create quadrangles instead of triangles, one can use the `setRecombine’
# # constraint:
gmsh.model.geo.mesh.setRecombine(2, 1)
gmsh.model.geo.synchronize()
# Finally we apply an elliptic smoother to the grid to have a more regular
# mesh:
# gmsh.option.setNumber(“Mesh.Smoothing”, 100)
gmsh.model.mesh.generate(2)
#===============testing
# if ‘close’ not in sys.argv:
# gmsh.fltk.run() # Write mesh data:
# gmsh.write(“mymeshmsh”)
# gmsh.finalize()
#=============end testing
#-----------fenicsx
print(“Creating the mesh in fenicsx”)
from dolfinx.io import gmshio
from mpi4py import MPI
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD
domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model,
mesh_comm,
gmsh_model_rank,
gdim=2)
print(‘done with the mesh…’)
total_num_of_elems = num_elem*num_elem
print(‘The number of elments in the L-shape domain:’,total_num_of_elems)
# Dimension of the space
tdim = domain.topology.dim
# print(tdim)
topology, cell_types, geometry = plot.create_vtk_mesh(domain, tdim)
print(‘Constructing the mesh …’)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
# To view the mesh
plotter.show()
print(‘done…’)
return domain, tdim, total_num_of_elems, topology, cell_types, geometry

Blockquote
num_elemnts = 6
RT_order = 1 # order=1 is RT_0 lowest order RT for xfenics
element = ‘lshape_unstr_quad’
print(“number of elements %i”%num_elemnts)
‘’’ The mesh’‘’
domain, tdim, total_num_of_elems, ,,_ = lshape_unstructQuad_gmsh2(element, num_elemnts)
‘’‘The function spaces’‘’
print(‘Constructing the space …’)
space = fem.FunctionSpace(domain, (“RT”, RT_order))
print(‘done …\n’)
‘’‘Trial and test functions’‘’
print(‘The trial and test functions…’)
u = TrialFunction(space)
v = TestFunction(space)
print(‘done…\n’)
‘’‘The operators and the bilinear form’‘’
‘’‘Left hand side’‘’
print(‘Constructing the bilinear forms…’)
a = div(u)div(v)dx
stiff_bilinear_form = fem.form(a)
‘’‘Right hand side’‘’
b = inner(u,v)dx
mass_bilinear_form = fem.form(b)
print(‘done…\n’)
‘’‘----------------Setting the BC---------------------’‘’
print(‘Impossing the BC…’)
‘’’ Create facet to cell connectivity required to determine boundary facets’‘’
tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)
boundary_dofs = fem.locate_dofs_topological(space, fdim, boundary_facets)
ubc = fem.Function(space)
ubc.vector.set(0.0)
bc = fem.dirichletbc(ubc, boundary_dofs)
local_range = space.dofmap.index_map.local_range
dofs = np.arange(local_range)
print(‘INFO :The Dofs ‘, dofs)
‘’‘Assempleing the matrices’’’
print(‘Assembling the system…’)
A = fem.petsc.create_matrix(stiff_bilinear_form, )
A.setOption(PETSc.Mat.Option.SYMMETRIC, True)
A.setOption(PETSc.Mat.Option.SYMMETRY_ETERNAL, True)
A.setOption(PETSc.Mat.Option.IGNORE_ZERO_ENTRIES, True)
fem.petsc.assemble_matrix(A, stiff_bilinear_form,bcs=[bc])
A.assemble()
A.setOption(PETSc.Mat.Option.NEW_NONZERO_LOCATIONS, False)
‘’’ Mass matrix ‘’’
B = fem.petsc.create_matrix(mass_bilinear_form)
B.setOption(PETSc.Mat.Option.SYMMETRIC, True)
B.setOption(PETSc.Mat.Option.SYMMETRY_ETERNAL, True)
B.setOption(PETSc.Mat.Option.IGNORE_ZERO_ENTRIES, True)
fem.petsc.assemble_matrix(B, mass_bilinear_form,bcs=[bc])
B.assemble()
B.setOption(PETSc.Mat.Option.NEW_NONZERO_LOCATIONS, False)
print(‘done…\n’)
‘’‘Imposing the zero boundary conditions ‘’’
print(‘Imposing the zero Boundary conditions’)
B.zeroRowsLocal(bc.dof_indices()[0], 1.)
print(‘done…\n’)
print(‘Setting the solver…’)
shift = SLEPc.ST().create(MPI.COMM_WORLD)
shift.setType(‘sinvert’) # spectral transform
shift.setShift(1/factor) # spectral shift
eigensolver = SLEPc.EPS().create(MPI.COMM_WORLD)
numb_eigs = 1000
eigensolver.setDimensions(numb_eigs) # set number of eigenvalues to compute
eigensolver.setProblemType(eigensolver.ProblemType.GNHEP)
eigensolver.setWhichEigenpairs(eigensolver.Which.TARGET_MAGNITUDE) # For shift-and-invert
eigensolver.setST(shift)
eigensolver.setOperators(A,B)
eigensolver.setFromOptions() #any options specified at run time in the command line are
print(‘done…\n’)
print(‘Solving the problem…’)
eigensolver.solve()
print(‘done…\n’)
print("“)
print(” SLEPc Solution Results “)
print(”****“)
print()
num_of_converged_eig_val = eigensolver.getConverged()
vr, vi = A.createVecs()
real_eigs_sorted =
loop = 0
print( “Number of converged eigenpairs %d” % num_of_converged_eig_val )
print(‘’)
if num_of_converged_eig_val > 0:
for i in range (num_of_converged_eig_val):
l = eigensolver.getEigenpair(i ,vr, vi)
if element == ‘lshape_unstr_quad’:
if l.real > 1.4:
print(f"Mode {i} with value {l.real}”)
real_eigs_sorted += [l.real]
loop +=1
if loop == 20: #This is tocontrol the nukmber f eigvalues, not to spit them all
break

Blockquote
#To plot the first eigfunction
lmbda = eigensolver.getEigenpair(0, vr, vi)
x = SpatialCoordinate(domain)
topology1, cell_types1, x = plot.create_vtk_mesh(space)
grid1 = pyvista.UnstructuredGrid(topology1, cell_types1, x)
grid.point_data[“u”] = vr.x.array
warped = grid1.warp_by_scalar(“u”, factor=25)
plotter1 = pyvista.Plotter()
plotter1.add_mesh(warped, show_edges=True, show_scalar_bar=True, scalars=“u”)
plotter1.show()

The error:
raise RuntimeError(“Can only create meshes from continuous or discontinuous Lagrange spaces”)
RuntimeError: Can only create meshes from continuous or discontinuous Lagrange spaces

8

Please do not use blockquote.
Use 3x`, i.e.

```python
# Add code here
```

The error you are getting is because RT spaces do not have dof coordinates, (as their functionals are integrals: DefElement)

Thus you would have to interpolate your solution into an appropriate space (Say a vector space with DG 2 elements) to visualize the solution.

9

I managed to plot the solution But I think have a problem. How can I impose zero Neumann boundary conditions for my problem I need
u . n = 0.

10

This is not a Neumann condition.

This is a dirichlet condition that you would apply to a dof associated with a facet (DefElement)
if you weakly enforce this normal component to be zero, with nitsches method (as you are solving an eigenvalue problem where strong enforcement of DIrichlet bc doesn’t necessarily make sense).

See for instance:
https://bitbucket.org/nate-sime/dolfin_dg/src/master/

11

Thank you for the correction, yes you are right. I didn’t know how to apply it on the boundary.
Can you help me with this?
Moreover, if I want to use a quad element, does RT in fenicsx support this?

12

Yes, see: DefElement: Nédélec (first kind)
for definitions: (i.e. "RTCE" (quadrilateral, Lagrange))

13

Hi Dokken,
I’m still stuck in plotting the eigenfunction of slepsc with RT.
The code is above and the plotting function is:

def fig_out2(eig_vect, domain, space):
	u = Function(space)
	u.vector.setArray(eig_vect.array)

	gdim = domain.geometry.dim
	V0 = fem.FunctionSpace(domain, ("Discontinuous Lagrange", 2))
	u0 = fem.Function(V0, dtype=np.float64)
	u0.interpolate(u)

	import pyvista
	plotter1 = pyvista.Plotter()

	topology1, cell_types1, x1 = plot.create_vtk_mesh(V0)
	grid1 = pyvista.UnstructuredGrid(topology1, cell_types1, x1)
	grid1.point_data["u"] = u0.x1.array.reshape(x.shape[0], V0.dofmap.index_map_bs)
	glyphs = grid1.glyph(orient="u", factor=0.1)
	plotter1.add_mesh(glyphs)
	plotter1.show()

1-If I choose V0 such that:
V0 = fem.FunctionSpace(domain, (“Discontinuous Lagrange”, 2),(gdim,))

I get the following error:
assert mesh is None
AssertionError
2- But if I choose:
V0 = fem.FunctionSpace(domain, (“Discontinuous Lagrange”, 2))
RuntimeError: Interpolation: elements have different value dimensions

14

This should be a vector function space

15

Hi Dokken,
I still have an issue with plotting
This is what I did to plot:

def fig_out2(eig_vect, domain, space):

	#To plot the eigenfunctions, eig_vect from SLEPC solver
	u = Function(space)
	u.vector.setArray(eig_vect.array)

	#---------------------------------
	import pyvista
	plotter = pyvista.Plotter()
	pyvista_topology, pyvista_cell_types, x = plot.create_vtk_mesh(domain)
	grid = pyvista.UnstructuredGrid(pyvista_topology, 
									pyvista_cell_types, 
									x)
	
	plotter.add_mesh(grid, show_edges=True)

	# Exact visualization of the RT spaces requires a Lagrange or
	# discontinuous Lagrange finite element functions. Therefore, we
	# interpolate the RT function into a 2nd-order discontinuous
	# Lagrange space.
	gdim = domain.geometry.dim
	V0 = fem.VectorFunctionSpace(domain, ("Discontinuous Lagrange", 2))
	u0 = fem.Function(V0, dtype=np.float64)
	u0.interpolate(u)

	# Create a second grid, whose geometry and topology are based on the
	# output function space
	topology, cell_types, x = plot.create_vtk_mesh(V0)
	grid = pyvista.UnstructuredGrid(topology, cell_types, x)
	
	# Create point cloud of vertices, and add the vertex values to the cloud
	grid.point_data["u"] = u0.x.array.reshape(x.shape[0], V0.dofmap.index_map_bs)
	glyphs = grid.glyph(orient="u", factor=0.1)

On the last step, I get the following error

ValueError: Data field (u) with type (FieldAssociation.POINT) could not be set as the active vectors

I have a problem with the dimensions for each space, RT, and vector DG. They don’t seem compatible.
I’m using the source
Visualization with PyVista — DOLFINx 0.8.0.0 documentation

16

Please read: Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial

17

Thank you very much.
I managed to plot it with the last documentation.