Calculate flux over boundary

variational formulation
1

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 = fem.Expression(ufl.inner(ufl.grad(uh),n), V.element.interpolation_points())
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

2

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
3

Thanks for your reply

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], 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[1], 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)
a = dot(grad(u), grad(v)) * dx


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 = Expression(inner(grad(uh),n), V.element.interpolation_points())
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.

4

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)

5

thanks for the answer.
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[2]/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[0]/u_norm, nh[1]/u_norm)), V.element.interpolation_points())
    elif dim == 3:
        nh_normalized = dfx.fem.Expression(ufl.as_vector((nh[0]/u_norm, nh[1]/u_norm, nh[2]/u_norm)), V.element.interpolation_points())
    else:
        raise ValueError("Unsupported mesh dimension: {}".format(dim))

many thanks :slight_smile:

6

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

Any suggestions beside adding epsilon?

7

Hi p1rx,

Thanks for the input on the code not being applicable in other than 3 dimensions, I have now updated the gist so that the dimension of the mesh is taken into account when normalizing the facet vectors.

I’m not sure whether I’ve understood your problem correctly, but I think a solution could be to assemble the flux expression only on the boundary (i.e. only assemble the flux expression on parts of the mesh where you have computed facet vector approximations). Another way could be to take the vector you have now and extract the values at boundary dofs, such as to discard the non-boundary values of the vector.

Cheers, H