Hi,

this works exactly the same, independently of mechanics. If you take the action of the thermal bilinear form on a field with unit value on a given boundary, you will get the total flux on this boundary in a consistent manner.

1 Like

Thanks a lot Jeremy, I am very excited, I will post my (failed) attempt when I have the time. So far my whole personal project solely depends on evaluating these fluxes as accurately as possible, I was thinking about ways to bypass the problem but I did not find any remedy. I now have big hopes!

My current attempt looks like so:

```
import numpy as np
from dolfinx import log, fem
from dolfinx.fem import (Constant, Function, functionspace, FunctionSpace, assemble_scalar,
dirichletbc, form, locate_dofs_geometrical, VectorFunctionSpace,
locate_dofs_topological, Expression, create_matrix, assemble_matrix)
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
from mpi4py import MPI
from petsc4py import PETSc
from petsc4py.PETSc import ScalarType
from ufl import (FiniteElement, Measure, MixedElement, SpatialCoordinate,
TestFunction, TrialFunction, as_tensor, dot, dx, grad, inner,
inv, split, derivative, FacetNormal, div, curl, as_vector, action)
import gmsh
from dolfinx.io import gmshio, VTXWriter
mesh, cell_markers, facet_markers = gmshio.read_from_msh("meshes/mesh2.msh", MPI.COMM_WORLD, gdim=2)
def thermal_problem(ΔT=0):
current_curves=(1,2)
voltages_curves=(1,2)
inj_current_curve, out_current_curve = current_curves
reading_voltage_curve_0, reading_voltage_curve_1 = voltages_curves
# Define FE function space
deg = 1
el = FiniteElement("CG", mesh.ufl_cell(), deg)
V = FunctionSpace(mesh, el)
u = TestFunction(V)
temp = Function(V)
κ = 1.8
T_cold = 300
# Define the boundary conditions
left_facets = facet_markers.find(inj_current_curve)
right_facets = facet_markers.find(out_current_curve)
left_dofs_temp = locate_dofs_topological(V, mesh.topology.dim-1, left_facets)
right_dofs_temp = locate_dofs_topological(V, mesh.topology.dim-1, right_facets)
bc_temp_left = dirichletbc(ScalarType(T_cold), left_dofs_temp, V)
bc_temp_right = dirichletbc(ScalarType(T_cold + ΔT), right_dofs_temp, V)
bcs = [bc_temp_left, bc_temp_right]
x = SpatialCoordinate(mesh)
dx = Measure("dx", domain=mesh, subdomain_data=cell_markers)
ds = Measure("ds", domain=mesh, subdomain_data=facet_markers)
# Weak form.
weak_form = dot(-κ * grad(temp), grad(u)) * dx
print(f''' ------- Pre-processing --------
Length of the side where heat enters: {assemble_scalar(form(1 * ds(inj_current_curve, domain=mesh)))}
Length of the side where heat leaves the material: {assemble_scalar(form(1 * ds(out_current_curve, domain=mesh)))}
''')
# Solve the PDE.
problem = NonlinearProblem(weak_form, temp, bcs=bcs)
solver = NewtonSolver(MPI.COMM_WORLD, problem)
ksp = solver.krylov_solver
opts = PETSc.Options()
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.WARNING)
n, converged = solver.solve(temp)
assert (converged)
print(f'''------- Processing --------
Number of interations: {n:d}''')
avg_temp = assemble_scalar(form(temp * dx(domain=mesh))) / assemble_scalar(form(1 * dx(domain=mesh)))
print(f'avg_temp: {avg_temp}')
# Compute fluxes on boundaries
n = FacetNormal(mesh)
down_heat_flux = form(dot(-κ * grad(temp), n)*ds(inj_current_curve))
Q_out = assemble_scalar(down_heat_flux)
top_heat_flux = form(dot(-κ * grad(temp), n)*ds(out_current_curve))
Q_in = assemble_scalar(top_heat_flux)
print(f'''------- Post processing --------
Q_in: {Q_in}
Q_out: {Q_out}''')
Q = functionspace(mesh, ("DG", deg-1, (mesh.geometry.dim,)))
q = Function(Q)
flux_calculator = Expression(-κ * grad(temp), Q.element.interpolation_points())
q.interpolate(flux_calculator)
# Try an alternate way to compute the fluxes.
virtual_temp = Function(V)
virtual_heat_flux_form = form(action(weak_form, virtual_temp))
virtual_heat_flux = assemble_scalar(virtual_heat_flux_form)
print(f'Virtual heat flux computation: {virtual_heat_flux}')
return Q_in, Q_out
t_b = thermal_problem(ΔT=10)
print(t_b)
```

Producing the output:

```
------- Pre-processing --------
Length of the side where heat enters: 0.10000000000000002
Length of the side where heat leaves the material: 0.09999999999999998
------- Processing --------
Number of interations: 1
avg_temp: 305.4445644560497
------- Post processing --------
Q_in: -7.434551004310087
Q_out: 7.494950253817867
Virtual heat flux computation: 0.0
(-7.434551004310087, 7.494950253817867)
```

In other words, it computes as 0 heat flux instead of something like 7.49 or so.

My mesh is similar to the file.msh I posted in post #10 just above.

Any help is appreciated.