# How to extract and plot the eigen-function

**URL:** <https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272>\
**Category:** dolfinx\
**Created:** [September 23, 2023, 11:44pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272 "2023-09-23T23:44:13Z")\
**Posts on this page:** 17\
**Page:** 1

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [September 23, 2023, 11:44pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/1 "2023-09-23T23:44:13Z")

</div>

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](https://fenicsproject.discourse.group/t/a-simple-eigenvalue-solver-in-dolfin-x/3790/3)]

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"

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 24, 2023, 6:09am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/2 "2023-09-24T06:09:14Z")

</div>

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

[https://jsdokken.com/dolfinx-tutorial/](https://jsdokken.com/dolfinx-tutorial/)  
Or

> <https://github.com/FEniCS/dolfinx/blob/main/python/demo/demo_pyvista.py>

for examples

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [September 24, 2023, 6:47am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/3 "2023-09-24T06:47:05Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 24, 2023, 7:05am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/4 "2023-09-24T07:05:01Z")

</div>

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](https://github.com/FEniCS/dolfinx/blob/v0.6.0/python/demo/demo_pyvista.py)

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [September 24, 2023, 8:04am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/5 "2023-09-24T08:04:20Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 24, 2023, 8:06am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/6 "2023-09-24T08:06:53Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [September 25, 2023, 9:11am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/7 "2023-09-25T09:11:42Z")

</div>

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](https://github.com/FEniCS/dolfinx/)")

> 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](http://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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 25, 2023, 10:46am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/8 "2023-09-25T10:46:11Z")

</div>

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

````python
```python
# Add code here
```

````

The error you are getting is because `RT` spaces do not have dof coordinates, (as their functionals are integrals: [DefElement](https://defelement.com/elements/examples/triangle-raviart-thomas-lagrange-1.html))

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

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [September 27, 2023, 8:36am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/9 "2023-09-27T08:36:43Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 27, 2023, 11:59am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/10 "2023-09-27T11:59:00Z")

</div>

> [@Lio97](#):
>
> How can I impose zero Neumann boundary conditions for my problem I need  
> u . n = 0.

This is not a Neumann condition.

This is a dirichlet condition that you would apply to a dof associated with a facet ([DefElement](https://defelement.com/elements/examples/triangle-raviart-thomas-lagrange-1.html))  
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/](https://bitbucket.org/nate-sime/dolfin_dg/src/master/)

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [September 28, 2023, 8:26am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/11 "2023-09-28T08:26:08Z")

</div>

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?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 28, 2023, 9:11am UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/12 "2023-09-28T09:11:06Z")

</div>

Yes, see: [DefElement: Nédélec (first kind)](https://defelement.com/elements/nedelec1.html)  
for definitions: (i.e. `"RTCE"` (quadrilateral, Lagrange))

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [October 12, 2023, 1:08pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/13 "2023-10-12T13:08:08Z")

</div>

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

```auto
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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 12, 2023, 2:29pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/14 "2023-10-12T14:29:40Z")

</div>

> [@Lio97](#):
>
> `V0 = fem.FunctionSpace(domain, ("Discontinuous Lagrange", 2))`

This should be a vector function space

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [October 16, 2023, 5:24pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/15 "2023-10-16T17:24:11Z")

</div>

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

```auto
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](https://docs.fenicsproject.org/dolfinx/main/python/demos/demo_pyvista.html)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 16, 2023, 5:38pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/16 "2023-10-16T17:38:53Z")

</div>

> [@Lio97](#):
>
> ValueError: Data field (u) with type (FieldAssociation.POINT) could not be set as the active

Please read: [Test problem 1: Channel flow (Poiseuille flow) — FEniCSx tutorial](https://jsdokken.com/dolfinx-tutorial/chapter2/ns_code1.html#visualization-of-vectors)

---

<div class="post-metadata">

**Author:** ![Lio97](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/lio97/32/5147_2.png) [@Lio97](https://fenicsproject.discourse.group/u/Lio97)\
**Post date:** [October 16, 2023, 6:24pm UTC](https://fenicsproject.discourse.group/t/how-to-extract-and-plot-the-eigen-function/12272/17 "2023-10-16T18:24:27Z")

</div>

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