I am trying to convert some old (FEniCS) working code to the newer FEniCSx, but I am facing a problem. The temperature distribution I get is insensitive to the electric current, but this should not be the case, since I consider Joule heating.

In terms of code, I employ mixed elements and solve the heat equation (to get the temperature distribution), and div(J)=0, where J = sigma E. I expect a parabolic solution for the temperature distribution along z. As boundary conditions, I keep fixed the temperatures of the left and right ends of a wire (parallelepipede).

My full code starts with the file.msh:

```
$MeshFormat
4.1 0 8
$EndMeshFormat
$PhysicalNames
3
2 14 "left_end"
2 15 "right_end"
3 13 "the_wire"
$EndPhysicalNames
$Entities
8 12 6 1
1 0 0 1 0
2 0 0 0 0
3 0 0.01 1 0
4 0 0.01 0 0
5 0.01 0 1 0
6 0.01 0 0 0
7 0.01 0.01 1 0
8 0.01 0.01 0 0
1 -1e-07 -1e-07 -9.999999994736442e-08 1e-07 1e-07 1.0000001 0 2 2 -1
2 -1e-07 -1.000000000002735e-07 0.9999999000000001 1e-07 0.0100001 1.0000001 0 2 1 -3
3 -1e-07 0.009999900000000001 -9.999999994736442e-08 1e-07 0.0100001 1.0000001 0 2 4 -3
4 -1e-07 -1.000000000002735e-07 -1e-07 1e-07 0.0100001 1e-07 0 2 2 -4
5 0.009999900000000001 -1e-07 -9.999999994736442e-08 0.0100001 1e-07 1.0000001 0 2 6 -5
6 0.009999900000000001 -1.000000000002735e-07 0.9999999000000001 0.0100001 0.0100001 1.0000001 0 2 5 -7
7 0.009999900000000001 0.009999900000000001 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 0 2 8 -7
8 0.009999900000000001 -1.000000000002735e-07 -1e-07 0.0100001 0.0100001 1e-07 0 2 6 -8
9 -1.000000000002735e-07 -1e-07 -1e-07 0.0100001 1e-07 1e-07 0 2 2 -6
10 -1.000000000002735e-07 -1e-07 0.9999999000000001 0.0100001 1e-07 1.0000001 0 2 1 -5
11 -1.000000000002735e-07 0.009999900000000001 -1e-07 0.0100001 0.0100001 1e-07 0 2 4 -8
12 -1.000000000002735e-07 0.009999900000000001 0.9999999000000001 0.0100001 0.0100001 1.0000001 0 2 3 -7
1 -1e-07 -1.000000000002735e-07 -9.999999994736442e-08 1e-07 0.0100001 1.0000001 0 4 -1 4 3 -2
2 0.009999900000000001 -1.000000000002735e-07 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 0 4 -5 8 7 -6
3 -1.000000000002735e-07 -1e-07 -9.999999994736442e-08 0.0100001 1e-07 1.0000001 0 4 -9 1 10 -5
4 -1.000000000002735e-07 0.009999900000000001 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 0 4 -11 3 12 -7
5 -1.000000000002735e-07 -1.000000000002735e-07 -1e-07 0.0100001 0.0100001 1e-07 1 14 4 -4 9 8 -11
6 -1.000000000002735e-07 -1.000000000002735e-07 0.9999999000000001 0.0100001 0.0100001 1.0000001 1 15 4 -2 10 6 -12
1 -1.000000000002735e-07 -1.000000000002735e-07 -9.999999994736442e-08 0.0100001 0.0100001 1.0000001 1 13 6 -1 2 -3 4 -5 6
$EndEntities
$Nodes
15 44 1 44
0 1 0 1
1
0 0 1
0 2 0 1
2
0 0 0
0 3 0 1
3
0 0.01 1
0 4 0 1
4
0 0.01 0
0 5 0 1
5
0.01 0 1
0 6 0 1
6
0.01 0 0
0 7 0 1
7
0.01 0.01 1
0 8 0 1
8
0.01 0.01 0
1 1 0 9
9
10
11
12
13
14
15
16
17
0 0 0.09999999999999987
0 0 0.1999999999999997
0 0 0.2999999999999997
0 0 0.3999999999999997
0 0 0.4999999999999998
0 0 0.5999999999999998
0 0 0.6999999999999998
0 0 0.7999999999999999
0 0 0.8999999999999999
1 3 0 9
18
19
20
21
22
23
24
25
26
0 0.01 0.09999999999999987
0 0.01 0.1999999999999997
0 0.01 0.2999999999999997
0 0.01 0.3999999999999997
0 0.01 0.4999999999999998
0 0.01 0.5999999999999998
0 0.01 0.6999999999999998
0 0.01 0.7999999999999999
0 0.01 0.8999999999999999
1 5 0 9
27
28
29
30
31
32
33
34
35
0.01 0 0.09999999999999987
0.01 0 0.1999999999999997
0.01 0 0.2999999999999997
0.01 0 0.3999999999999997
0.01 0 0.4999999999999998
0.01 0 0.5999999999999998
0.01 0 0.6999999999999998
0.01 0 0.7999999999999999
0.01 0 0.8999999999999999
1 7 0 9
36
37
38
39
40
41
42
43
44
0.01 0.01 0.09999999999999987
0.01 0.01 0.1999999999999997
0.01 0.01 0.2999999999999997
0.01 0.01 0.3999999999999997
0.01 0.01 0.4999999999999998
0.01 0.01 0.5999999999999998
0.01 0.01 0.6999999999999998
0.01 0.01 0.7999999999999999
0.01 0.01 0.8999999999999999
2 5 0 0
2 6 0 0
3 1 0 0
$EndNodes
$Elements
3 64 1 64
2 5 2 2
1 2 4 6
2 6 4 8
2 6 2 2
3 1 5 3
4 5 7 3
3 1 4 60
5 9 2 4 6
6 4 9 6 27
7 27 6 4 8
8 4 9 27 18
9 8 4 27 18
10 8 27 36 18
11 10 9 18 27
12 18 10 27 28
13 28 27 18 36
14 18 10 28 19
15 36 18 28 19
16 36 28 37 19
17 11 10 19 28
18 19 11 28 29
19 29 28 19 37
20 19 11 29 20
21 37 19 29 20
22 37 29 38 20
23 12 11 20 29
24 20 12 29 30
25 30 29 20 38
26 20 12 30 21
27 38 20 30 21
28 38 30 39 21
29 13 12 21 30
30 21 13 30 31
31 31 30 21 39
32 21 13 31 22
33 39 21 31 22
34 39 31 40 22
35 14 13 22 31
36 22 14 31 32
37 32 31 22 40
38 22 14 32 23
39 40 22 32 23
40 40 32 41 23
41 15 14 23 32
42 23 15 32 33
43 33 32 23 41
44 23 15 33 24
45 41 23 33 24
46 41 33 42 24
47 16 15 24 33
48 24 16 33 34
49 34 33 24 42
50 24 16 34 25
51 42 24 34 25
52 42 34 43 25
53 17 16 25 34
54 25 17 34 35
55 35 34 25 43
56 25 17 35 26
57 43 25 35 26
58 43 35 44 26
59 1 17 26 35
60 26 1 35 5
61 5 35 26 44
62 26 1 5 3
63 44 26 5 3
64 44 5 7 3
$EndElements
```

Then my FEniCSx code:

```
from mpi4py import MPI
import gmsh
from dolfinx.io import gmshio
mesh, cell_markers, facet_markers = gmshio.read_from_msh("long_wire.msh", MPI.COMM_WORLD, gdim=3)
proc = MPI.COMM_WORLD.rank
print(proc)
import meshio
def create_mesh(mesh, cell_type, prune_z=False):
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("long_wire.msh")
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "triangle", prune_z=False)
tetrahedron_mesh = create_mesh(msh, "tetra", prune_z=False)
meshio.write("mt.xdmf", triangle_mesh)
meshio.write("mesh.xdmf", tetrahedron_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, MixedElement, FiniteElement, split)
import numpy as np
# Define ME function space
el = FiniteElement("CG", mesh.ufl_cell(), 1)
mel = MixedElement([el, el])
ME = FunctionSpace(mesh, mel)
u, v = TestFunction(ME)
TempVolt = Function(ME)
temp, volt = split(TempVolt)
rho = 1e-9
sigma = 1.0 / rho
kappa = 144
# Reminder of the mesh:
#Physical Volume("the_wire", 13) = {1};
#Physical Surface("left_end", 14) = {5};
#Physical Surface("right_end", 15) = {6};
inj_current_surface = 14
vanish_voltage_surface = 15 # 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 = 14
reading_voltage_surface_1 = 15
# Define the boundary conditions
from petsc4py.PETSc import ScalarType
left_facets = ft.find(inj_current_surface)
right_facets = ft.find(vanish_voltage_surface)
left_dofs = locate_dofs_topological(ME.sub(1), mesh.topology.dim-1, left_facets)
bc_volt = dirichletbc(ScalarType(V_current_contact_out), left_dofs, ME.sub(1))
left_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, left_facets)
right_dofs_temp = locate_dofs_topological(ME.sub(0), mesh.topology.dim-1, right_facets)
bc_temp_left = dirichletbc(ScalarType(300.0), left_dofs_temp, ME.sub(0))
bc_temp_right = dirichletbc(ScalarType(320.0), right_dofs_temp, ME.sub(0))
bcs = [bc_temp_left, bc_temp_right, bc_volt]
from dolfinx.fem import assemble_scalar, form
from ufl import Measure
x = SpatialCoordinate(mesh)
print(x)
#print(x.shape)
dx = Measure("dx", domain=mesh,subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)
the_current = 20.0 # Current, in amperes.
J = the_current / assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh)))
print('Current density:', J)
print('Surface area where current is injected', assemble_scalar(form(1 * ds(inj_current_surface, domain=mesh))))
print('Surface area where current leaves the wire', assemble_scalar(form(1 * ds(vanish_voltage_surface, domain=mesh))))
from ufl import dot
from petsc4py import PETSc
def temp_init(x):
values = np.full(x.shape[1], 300, dtype=PETSc.ScalarType)
return values
def volt_init(x):
values = np.full(x.shape[1], 1.0e-3, dtype=PETSc.ScalarType)
return values
# Initialize values for the temperature and voltage.
TempVolt.sub(0).interpolate(temp_init)
TempVolt.sub(1).interpolate(volt_init)
# Weak form.
J_vector = -sigma * grad(volt)
Fourier_term = dot(kappa * grad(temp), grad(u)) * dx
Joule_term = dot(rho * J_vector, J_vector) * u * dx
F_T = Fourier_term + Joule_term
F_V = dot(grad(v), sigma * grad(volt))*dx + v * J * ds(vanish_voltage_surface)
weak_form = F_T + F_V
# Solve the PDE.
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
problem = NonlinearProblem(weak_form, TempVolt, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-4
solver.report = True
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(TempVolt)
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")
with io.XDMFFile(mesh.comm, "solution/voltage.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(TempVolt.sub(1))
with io.XDMFFile(mesh.comm, "solution/temperature.xdmf", "w") as xdmf:
xdmf.write_mesh(mesh)
xdmf.write_function(TempVolt.sub(0))
```

I get a linear temperature distribution, from 300 K to 320 K, which indicates that the boundary conditions are correctly implemented, but not the influence of J on the temperature. I do expect a parabolic dependence, not a linear one. I do not see what is wrong, it is as if â€śtempâ€ť was not influenced by â€śJ_vectorâ€ť in my code. Any idea what could be wrong?