How to visualize/save a vector field in a particular case?

I am unable to visualize a vector field with pyvista, and to save it in a file (yes, I have already looked at the magnetostatics example in the tutorial, but it seems a very different case so that I cannot apply it, as is.). I am using a mesh constructed with Gmsh (I still haven’t figured out how to use the Python library for gmsh, in particular how can people figure out the numbers of the curves/lines if they do not run Gmsh GUI, but that’s for another topic I guess).

My fenicsx code solves Ohm’s law in a 2D mesh, given an input current. I.e. it computes the voltage reading between 2 surfaces. It also computes the electrostatics potential in the whole geometry, and then computes the electric field. And that’s the field I would like to visualize.

Here is my code, first the file.msh:

$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
4
1 13 "left_mat_A"
1 14 "right_mat_B"
2 11 "material A"
2 12 "material B"
$EndPhysicalNames
$Entities
6 7 2 0
2 1 0 0 0 
3 1 2 0 0 
5 -1 0 0 0 
8 -1 2 0 0 
9 0 0 0 0 
10 0 2 0 0 
9 -1.0000001 -9.999999994736442e-08 -1e-07 -0.9999999000000001 2.0000001 1e-07 1 14 2 8 -5 
10 0.9999999000000001 -9.999999994736442e-08 -1e-07 1.0000001 2.0000001 1e-07 1 13 2 3 -2 
14 -1e-07 -9.999999994736442e-08 -1e-07 1e-07 2.0000001 1e-07 0 2 10 -9 
15 -1.0000001 -1e-07 -1e-07 9.999999994736442e-08 1e-07 1e-07 0 2 5 -9 
16 -1.0000001 1.9999999 -1e-07 9.999999994736442e-08 2.0000001 1e-07 0 2 10 -8 
17 -9.999999994736442e-08 -1e-07 -1e-07 1.0000001 1e-07 1e-07 0 2 9 -2 
18 -9.999999994736442e-08 1.9999999 -1e-07 1.0000001 2.0000001 1e-07 0 2 3 -10 
1 -9.999999994736442e-08 -9.999999994736442e-08 -1e-07 1.0000001 2.0000001 1e-07 1 12 4 17 -10 18 14 
2 -1.0000001 -9.999999994736442e-08 -1e-07 9.999999994736442e-08 2.0000001 1e-07 1 11 4 15 -14 16 9 
$EndEntities
$Nodes
15 81 1 81
0 2 0 1
1
1 0 0
0 3 0 1
2
1 2 0
0 5 0 1
3
-1 0 0
0 8 0 1
4
-1 2 0
0 9 0 1
5
0 0 0
0 10 0 1
6
0 2 0
1 9 0 7
7
8
9
10
11
12
13
-1 1.75 0
-1 1.5 0
-1 1.25 0
-1 0.9999999999999998 0
-1 0.75 0
-1 0.5 0
-1 0.25 0
1 10 0 7
14
15
16
17
18
19
20
1 1.75 0
1 1.5 0
1 1.25 0
1 0.9999999999999998 0
1 0.75 0
1 0.5 0
1 0.25 0
1 14 0 7
21
22
23
24
25
26
27
0 1.75 0
0 1.5 0
0 1.25 0
0 0.9999999999999998 0
0 0.75 0
0 0.5 0
0 0.25 0
1 15 0 3
28
29
30
-0.75 0 0
-0.4999999999999999 0 0
-0.25 0 0
1 16 0 3
31
32
33
-0.2500000000000001 2 0
-0.5000000000000001 2 0
-0.75 2 0
1 17 0 3
34
35
36
0.2500000000000001 0 0
0.5000000000000001 0 0
0.75 0 0
1 18 0 3
37
38
39
0.75 2 0
0.4999999999999999 2 0
0.25 2 0
2 1 0 21
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
0.2500000000000001 0.25 0
0.2500000000000001 0.5 0
0.25 0.75 0
0.25 0.9999999999999998 0
0.25 1.25 0
0.25 1.5 0
0.25 1.75 0
0.5000000000000001 0.25 0
0.5 0.5 0
0.5 0.75 0
0.5 0.9999999999999998 0
0.5 1.25 0
0.4999999999999999 1.5 0
0.4999999999999999 1.75 0
0.75 0.25 0
0.75 0.5 0
0.75 0.75 0
0.75 0.9999999999999998 0
0.75 1.25 0
0.75 1.5 0
0.7499999999999999 1.75 0
2 2 0 21
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
-0.75 0.25 0
-0.75 0.5 0
-0.75 0.75 0
-0.75 0.9999999999999998 0
-0.75 1.25 0
-0.75 1.5 0
-0.75 1.75 0
-0.4999999999999999 0.25 0
-0.4999999999999999 0.5 0
-0.5 0.75 0
-0.5 0.9999999999999998 0
-0.5 1.25 0
-0.5 1.5 0
-0.5000000000000001 1.75 0
-0.25 0.25 0
-0.25 0.5 0
-0.25 0.75 0
-0.2500000000000001 0.9999999999999998 0
-0.25 1.25 0
-0.2500000000000001 1.5 0
-0.2500000000000001 1.75 0
$EndNodes
$Elements
4 144 1 144
1 9 1 8
1 4 7 
2 7 8 
3 8 9 
4 9 10 
5 10 11 
6 11 12 
7 12 13 
8 13 3 
1 10 1 8
9 2 14 
10 14 15 
11 15 16 
12 16 17 
13 17 18 
14 18 19 
15 19 20 
16 20 1 
2 1 2 64
17 5 34 27 
18 27 34 40 
19 27 40 26 
20 26 40 41 
21 26 41 25 
22 25 41 42 
23 25 42 24 
24 24 42 43 
25 24 43 23 
26 23 43 44 
27 23 44 22 
28 22 44 45 
29 22 45 21 
30 21 45 46 
31 21 46 6 
32 6 46 39 
33 34 35 40 
34 40 35 47 
35 40 47 41 
36 41 47 48 
37 41 48 42 
38 42 48 49 
39 42 49 43 
40 43 49 50 
41 43 50 44 
42 44 50 51 
43 44 51 45 
44 45 51 52 
45 45 52 46 
46 46 52 53 
47 46 53 39 
48 39 53 38 
49 35 36 47 
50 47 36 54 
51 47 54 48 
52 48 54 55 
53 48 55 49 
54 49 55 56 
55 49 56 50 
56 50 56 57 
57 50 57 51 
58 51 57 58 
59 51 58 52 
60 52 58 59 
61 52 59 53 
62 53 59 60 
63 53 60 38 
64 38 60 37 
65 36 1 54 
66 54 1 20 
67 54 20 55 
68 55 20 19 
69 55 19 56 
70 56 19 18 
71 56 18 57 
72 57 18 17 
73 57 17 58 
74 58 17 16 
75 58 16 59 
76 59 16 15 
77 59 15 60 
78 60 15 14 
79 60 14 37 
80 37 14 2 
2 2 2 64
81 3 28 13 
82 13 28 61 
83 13 61 12 
84 12 61 62 
85 12 62 11 
86 11 62 63 
87 11 63 10 
88 10 63 64 
89 10 64 9 
90 9 64 65 
91 9 65 8 
92 8 65 66 
93 8 66 7 
94 7 66 67 
95 7 67 4 
96 4 67 33 
97 28 29 61 
98 61 29 68 
99 61 68 62 
100 62 68 69 
101 62 69 63 
102 63 69 70 
103 63 70 64 
104 64 70 71 
105 64 71 65 
106 65 71 72 
107 65 72 66 
108 66 72 73 
109 66 73 67 
110 67 73 74 
111 67 74 33 
112 33 74 32 
113 29 30 68 
114 68 30 75 
115 68 75 69 
116 69 75 76 
117 69 76 70 
118 70 76 77 
119 70 77 71 
120 71 77 78 
121 71 78 72 
122 72 78 79 
123 72 79 73 
124 73 79 80 
125 73 80 74 
126 74 80 81 
127 74 81 32 
128 32 81 31 
129 30 5 75 
130 75 5 27 
131 75 27 76 
132 76 27 26 
133 76 26 77 
134 77 26 25 
135 77 25 78 
136 78 25 24 
137 78 24 79 
138 79 24 23 
139 79 23 80 
140 80 23 22 
141 80 22 81 
142 81 22 21 
143 81 21 31 
144 31 21 6 
$EndElements

Now the FEniCSx code:

from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("rectangles.msh", MPI.COMM_WORLD, gdim=2)
proc = MPI.COMM_WORLD.rank 
import meshio
def create_mesh(mesh, cell_type, prune_z=True):
    cells = mesh.get_cells_type(cell_type)
    cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
    points = mesh.points[:,:2] if prune_z else mesh.points
    out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read":[cell_data]})
    return out_mesh

if proc == 0:
    # Read in mesh
    msh = meshio.read("rectangles.msh")
   
    # Create and save one file for the mesh, and one file for the facets 
    triangle_mesh = create_mesh(msh, "triangle", prune_z=True)
    line_mesh = create_mesh(msh, "line", prune_z=True)
    meshio.write("mesh.xdmf", triangle_mesh)
    meshio.write("mt.xdmf", line_mesh)
# Now read the mesh files produced above.
from dolfinx.io import XDMFFile
with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "r") as xdmf:
    mesh = xdmf.read_mesh(name="Grid")
    ct = xdmf.read_meshtags(mesh, name="Grid")
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
with XDMFFile(MPI.COMM_WORLD, "mt.xdmf", "r") as xdmf:
    ft = xdmf.read_meshtags(mesh, name="Grid")
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, assemble_scalar,
                         form, locate_dofs_geometrical, locate_dofs_topological)
from ufl import (SpatialCoordinate, TestFunction, TrialFunction,
                 dx, grad, inner, as_tensor, inv)
import numpy as np

# Define function space
V = FunctionSpace(mesh, ("CG", 1))
volt = Function(V) # TrialFunction(V)
u = TestFunction(V)

rho_tensor = as_tensor([[1, 0], [0, 1]])
sigma_tensor = inv(rho_tensor)
# Reminder of the mesh:
# Physical Surface("material A", 11) = {2};
# Physical Surface("material B", 12) = {1};
# Physical Curve("left_mat_A", 13) = {8};
# Physical Curve("right_mat_B", 14) = {2};
inj_current_curve = 13
vanish_voltage_curve = 14 # This corresponds to the curve the current leaves the material.
V_current_contact_out = 0.0 # Voltage value of the curve where the current leaves the material.
reading_voltage_surface_0 = 13
reading_voltage_surface_1 = 14
# We define the boundary conditions
from petsc4py.PETSc import ScalarType
u_bc = Function(V)
left_facets = ft.find(inj_current_curve)
left_dofs = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
bcs = [dirichletbc(ScalarType(V_current_contact_out), left_dofs, V)]
from dolfinx.fem import assemble_scalar, form
from ufl import Measure
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
the_current = 2.0 # Current, in amperes.
J = the_current / assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))
print(J)
print('Length of curve where current is injected', assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh))))
print('Length of curve where current leaves the wire', assemble_scalar(form(1 * ds(vanish_voltage_curve, domain=mesh))))
from ufl import dot
# Weak form.
F = dot(grad(u), sigma_tensor * grad(volt))*dx + u * J * ds(vanish_voltage_curve)

# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
#solve(F == 0, volt, bcs, solver_parameters={"newton_solver":{"relative_tolerance": 1e-11, "absolute_tolerance": 2e-9, "maximum_iterations":1000}})
problem = NonlinearProblem(F, volt, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-4
solver.report = True

from petsc4py import PETSc
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "cg"
opts[f"{option_prefix}pc_type"] = "gamg"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
opts[f"{option_prefix}ksp_max_it"]= 10000
ksp.setFromOptions()
from dolfinx import log
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(volt)
assert(converged)
print(f"Number of interations: {n:d}")
lu_solver = ksp
viewer = PETSc.Viewer().createASCII("lu_output.txt")
lu_solver.view(viewer)
solver_output = open("lu_output.txt", "r")
for line in solver_output.readlines():
    print(line)
avg_voltage_curve_integral_0 = assemble_scalar(form(volt * ds(reading_voltage_surface_0, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_surface_0, domain=mesh)))
avg_voltage_curve_integral_1 = assemble_scalar(form(volt * ds(reading_voltage_surface_1, domain=mesh))) / assemble_scalar(form(1 * ds(reading_voltage_surface_1, domain=mesh)))
voltage_diff = avg_voltage_curve_integral_1 - avg_voltage_curve_integral_0
print(voltage_diff,'voltage measured (V)')

import pyvista
pyvista.set_jupyter_backend("pythreejs")


from dolfinx import plot
u_topology, u_cell_types, u_geometry = plot.create_vtk_mesh(V)

u_grid = pyvista.UnstructuredGrid(u_topology, u_cell_types, u_geometry)
u_grid.point_data["volt"] = volt.x.array.real
u_grid.set_active_scalars("volt")
u_plotter = pyvista.Plotter()
u_plotter.add_mesh(u_grid, show_edges=True)
u_plotter.view_xy()
if not pyvista.OFF_SCREEN:
    u_plotter.show()
# Let's visualize the electric field
E_field = -grad(volt)

I wish to visualize E_field in all the geometry, and save the corresponding file(s) (with vtk I guess, but I am unable to do it).

You need to interpolate E_field into a suitable vector space, which is similar to Electromagnetics example — FEniCSx tutorial
You just change the following


W = VectorFunctionSpace(mesh, ("DG", 0))
B = Function(W)
B_expr(E_field, W.element.interpolation_points())
B.interpolate(B_expr)

And follow what is in the tutorial

1 Like

Thank you very much dokken, this seems to have worked, i.e. I can now visualize the electric field over the mesh, this is really beautiful.

I am still having a hard time to save the E_field (and volt for that matter), into some vtk file, for post processing with Paraview. I don’t think this part is covered in that tutorial’s page. It is covered elsewhere in the tutorial, here is my attempt to adapt it to my code:

from dolfinx import io
with io.VTKFile(mesh.comm, "output.pvd", "w") as vtk:
    vtk.write([E_field._cpp_object])
with io.XDMFFile(mesh.comm, "output.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    xdmf.write_function(E_field)

However this yields:

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In [30], line 4
      2 from dolfinx import io
      3 with io.VTKFile(mesh.comm, "output.pvd", "w") as vtk:
----> 4     vtk.write([E_field._cpp_object])
      5 with io.XDMFFile(mesh.comm, "output.xdmf", "w") as xdmf:
      6     xdmf.write_mesh(domain)

RuntimeError: Unknown triangle layout. Number of nodes: 1

You Need to use XDMFFile to visualize DG 0 quantities.

For higher order Lagrange and DG spaces (order>1) you should use VTXWriter.

For first order (continuous) Lagrange spaces you can use any of the formats

1 Like

Thanks dokken, so I just commented out the lines

with io.VTKFile(mesh.comm, "output.pvd", "w") as vtk:
    vtk.write([E_field._cpp_object])

(and replaced “domain” by “mesh”) and this seems to have worked, I can use Paraview to visualize the E field. It looks a very bit different than in Python: the density of the arrows isn’t uniform in Paraview, while it is with Python. I wonder if something is wrong.


This depends on how you visualize glyphs in Paraview. You have not added any information about how you visualize it (what filters you use, what options you have set in the filters), and thus it is really hard to give any guidance.

I just loaded the file Efield.xdmf, and then set the default glyph filter.

Edit: Bingo… I found out an option to “show all points”, when I do so, I get exactly the same visual than with the Python code. Thanks!