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).
