# Calculate flux over boundary

Hi,
i have a problem regarding solution postprocessing. I would like to calculate the flux field over some boundary → grad(u)*n. I tried interpolation

``````n = -ufl.FacetNormal(Mesh1_MESH)
flux_function = dolfinx.fem.Function(V)
flux_function.interpolate(flux)
``````

but the program crashes while defining the fem.expression->

``````VerificationError: CompileError: command '/bin/x86_64-linux-gnu-gcc' failed with exit code 1
``````

Probably this problem appears from the way the ufl.FacetNormal is created and the fact that it works in different manner than i initially thought. When trying to interpolate with other types of vectorfields there is no problem.

I would really apprecieate any suggestions. Thanks in advance

Note - my aim is not to calculate the flux integral over boundary but to find the flux/flux magnitude at each point of named boundary/interface

Please create a minimal working example – a small self-contained program that other people can copy-paste and run to see whether they can recreate your error. Also please specify which version of FEniCS(x) you’re using.

1 Like

Sure thing. My dolfinx version → 0.7.1 from docker image. As for the example: lets assume that the problem is a simple mass transfer with Dbc = 0 on top and = 1 on the bottom edge :

``````from mpi4py import MPI
from dolfinx import mesh
import numpy as np
from dolfinx import default_scalar_type
from dolfinx.fem import (Constant, Function, FunctionSpace, Expression,
assemble_scalar, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square
from dolfinx.plot import vtk_mesh

from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, FacetNormal, inner

import numpy as np
import pyvista

domain = mesh.create_unit_square(MPI.COMM_WORLD, 8, 8, mesh.CellType.quadrilateral)
V = FunctionSpace(domain, ("Lagrange", 1))

tdim = domain.topology.dim
fdim = tdim - 1
domain.topology.create_connectivity(fdim, tdim)
boundary_facets = mesh.exterior_facet_indices(domain.topology)

def boundary_Dtop(x):
return np.isclose(x, 1)
dofs_Dzero = locate_dofs_geometrical(V, boundary_Dtop)
u_bc_zero = Constant(domain, 0.0)
bc_top = dirichletbc(u_bc_zero, dofs_Dzero, V)

def boundary_Dbottom(x):
return np.isclose(x, 0)
dofs_Dzero = locate_dofs_geometrical(V, boundary_Dbottom)
u_bc_zero = Constant(domain, 1.0)
bc_bottom = dirichletbc(u_bc_zero, dofs_Dzero, V)
bc = [bc_top, bc_bottom]

u = TrialFunction(V)
v = TestFunction(V)

x = SpatialCoordinate(domain)
f = Constant(domain, default_scalar_type(0))
L = f * v * dx

problem = LinearProblem(a, L, bcs=bc, petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()

``````

I would like to calculate spatial distribution of the flux across the border.
The issue ive encountered is with this simple approach to calculate the grad(u) * n and interpolate it on some function space:

``````n = FacetNormal(domain)
flux_function = Function(V)
flux_function.interpolate(flux)

``````

as i did in the past f.eg. when i wanted to assess the vector norm. But here i get error as stated in the previous post.
The error probably comes from the fact that the normal facet is not a vector field as i thought, but is some abstract representation calculated during the integration.

Never the less maybe someone have an idea how to solve this error or calculate the spatial distribution of flux accross the boundary.

For now, expression does not support visualization over facets. You would have to project a function into the appropriate space for now (see for instance:How to compute tangent vectors on facets (2D surface parametrization) - #3 by hherlyng)

Code under:
Facet normal and facet tangent vector approximation in dolfinx
gets the job done + its pretty informative as for how the dolfinx works.
But one observation for the author. I had a 2D mesh, and calling the:

``````nh/u_norm
``````

raises an error:

``````ValueError: Index out of bounds.
``````

My advice is to substitute ori normalisation expression with:

``````# Get the dimension of the mesh (2 for 2D, 3 for 3D)
dim = V.mesh.geometry.dim
# Normalize the vectors based on the dimension
if dim == 2:
nh_normalized = dfx.fem.Expression(ufl.as_vector((nh/u_norm, nh/u_norm)), V.element.interpolation_points())
elif dim == 3:
nh_normalized = dfx.fem.Expression(ufl.as_vector((nh/u_norm, nh/u_norm, nh/u_norm)), V.element.interpolation_points())
else:
raise ValueError("Unsupported mesh dimension: {}".format(dim))
``````

many thanks so the proposed above normalisation doesnt reproduce the: Index out of bounds error, but vectors outside of the boundary are returned as NaNs → which in turn means that every operation with that vectorfield is useless. Even the vector norm is returned as NaN scalar field.
Obviously when omitting the normalisation stage the vector field is fully defined and the problem is not prevalent.
The Nans are probably related to division by 0 in vertices with 0 norm vectors