How to compute flux at particular boundary facets?

General
1

Hello, I’m having issues converting legacy fenics code to fenicsx. In particular, I would like to compute the flux along particular boundary facets, for example:

import numpy as np
from dolfinx.fem import (Constant, Function, FunctionSpace, 
                         assemble_scalar, assemble_vector, assemble_matrix, dirichletbc, form, locate_dofs_geometrical)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary
from mpi4py import MPI
from petsc4py.PETSc import ScalarType
from ufl import SpatialCoordinate, TestFunction, TrialFunction, dot, ds, dx, grad, FacetNormal
import ufl

mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)
V = FunctionSpace(mesh, ("CG", 1))
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) * dx

def u_exact(x):
    return 1 + x[0]**2 + 2*x[1]**2

def boundary_D(x):
    return np.logical_or(np.isclose(x[0], 0), np.isclose(x[0],1))

dofs_D = locate_dofs_geometrical(V, boundary_D)
u_bc = Function(V)
u_bc.interpolate(u_exact)
bc = dirichletbc(u_bc, dofs_D)

x = SpatialCoordinate(mesh)
g = -4 * x[1]
f = Constant(mesh, ScalarType(-6))
L = f * v * dx - g * v * ds

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

def neumann_boundary(x):
    return np.isclose(x[1], 1)

facets = locate_entities_boundary(mesh, dim=(mesh.topology.dim - 1),
                                       marker=neumann_boundary)
n = FacetNormal(mesh)
ds2 = ufl.Measure("ds", domain=mesh, subdomain_data=facets)
flux = form(dot(grad(uh), n) * ds2)
assemble_scalar(flux)

However, when I run this, I get the following error:

TypeError: init(): incompatible constructor arguments. The following argument types are supported:
1. dolfinx.cpp.fem.Form_float64(spaces: List[dolfinx::fem::FunctionSpace], integrals: Dict[dolfinx::fem::IntegralType, Tuple[List[Tuple[int, object]], dolfinx.cpp.mesh.MeshTags_int32]], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], need_permutation_data: bool, mesh: dolfinx.cpp.mesh.Mesh = None)
2. dolfinx.cpp.fem.Form_float64(form: int, spaces: List[dolfinx::fem::FunctionSpace], coefficients: List[dolfinx.cpp.fem.Function_float64], constants: List[dolfinx.cpp.fem.Constant_float64], subdomains: Dict[dolfinx::fem::IntegralType, dolfinx.cpp.mesh.MeshTags_int32], mesh: dolfinx.cpp.mesh.Mesh)

Invoked with: <cdata ‘uintptr_t’ 11305403016>, , [<dolfinx.cpp.fem.Function_float64 object at 0x2c60fae30>], , {<IntegralType.cell: 0>: None, <IntegralType.exterior_facet: 1>: array([192, 218, 241, 261, 278, 292, 303, 311, 316, 319], dtype=int32), <IntegralType.interior_facet: 2>: None, <IntegralType.vertex: 3>: None}, <dolfinx.mesh.Mesh object at 0x2ca7a0b30>

I had similar code that got the job done in legacy fenics, but I can’t figure this fenicsx stuff out! This was done with fenicsx/dolfinx 0.6.0, installed via conda-forge.

2

this should be

facets = locate_entities_boundary(mesh, dim=(mesh.topology.dim - 1),
                                       marker=neumann_boundary)
facet_tags = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, facets, np.full_like(facets, 1, dtype=np.int32))
n = FacetNormal(mesh)
ds2 = ufl.Measure("ds", domain=mesh, subdomain_data=facet_tags, subdomain_id=1)
flux = form(dot(grad(uh), n) * ds2)
3

Ahhh needed to use meshtags, very good. Thank you for the speedy reply!