Syntax around ufl.conditional

General
1

I have solved poisson problem in 3D and want to do further analysis that includes integration over a choice boundary (surface) to obtain area where values in \nabla u \cdot \hat{n} are less than or equal to \textbf{tol}

pval = \int_{\partial\Omega_1}(\nabla \cdot \hat{n} < tol)\text{ds(1)}

I thought I could use conditionals for this, but I cannot seem to figure out what the input types are for the functions from the ufl documentation ufl package — Unified Form Language (UFL) 2021.1.0 documentation.
Below is the minimum non-working example.

#!/usr/bin/env python3

import logging
import numpy as np
import ufl

from dolfinx import cpp, fem, io, mesh
from mpi4py import MPI
from petsc4py import PETSc


if __name__ == '__main__':
    logger = logging.getLogger()
    logger.setLevel(logging.DEBUG)
    formatter = logging.Formatter('%(levelname)s:%(asctime)s:%(message)s')
    fh = logging.FileHandler('transport.log')
    fh.setFormatter(formatter)
    logger.addHandler(fh)
    comm = MPI.COMM_WORLD

    logger.debug("Loading tetrahedra (dim = 3) mesh..")
    domain = mesh.create_box(comm, [[0, 0, 0], [10, 10, 10]], [20, 20, 20], mesh.CellType.tetrahedron)
    domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
    # meshtags = mesh.meshtags(domain, 2, ft.indices, ft.values)
    # Dirichlet BCs
    V = fem.FunctionSpace(domain, ("Lagrange", 2))
    u0 = fem.Function(V)
    with u0.vector.localForm() as u0_loc:
        u0_loc.set(1)

    u1 = fem.Function(V)
    with u1.vector.localForm() as u1_loc:
        u1_loc.set(0)

    left_boundary = mesh.locate_entities_boundary(domain, dim=(domain.topology.dim - 1),
                                                  marker=lambda x: np.isclose(x[2], 0.0))
    right_boundary = mesh.locate_entities_boundary(domain, dim=(domain.topology.dim - 1),
                                                  marker=lambda x: np.isclose(x[2], 10.0))
    insulated_boundary = mesh.locate_entities_boundary(domain, dim=(domain.topology.dim - 1),
                                                   marker=lambda x: np.logical_not(np.isclose(x[2], 0), np.isclose(x[2], 10)))
    lmarker = 1
    rmarker = 2
    imarker = 3

    left_bc = fem.dirichletbc(u0, fem.locate_dofs_topological(V, 2, left_boundary))
    right_bc = fem.dirichletbc(u1, fem.locate_dofs_topological(V, 2, right_boundary))
    n = ufl.FacetNormal(domain)
    x = ufl.SpatialCoordinate(domain)
    # ds = ufl.Measure("ds", domain=domain, subdomain_data=meshtags)

    # Define variational problem
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)

    # bulk conductivity [S.m-1]
    kappa = fem.Constant(domain, PETSc.ScalarType(0.1))
    f = fem.Constant(domain, PETSc.ScalarType(0.0))
    g = fem.Constant(domain, PETSc.ScalarType(0.0))

    a = ufl.inner(kappa * ufl.grad(u), ufl.grad(v)) * ufl.dx
    L = ufl.inner(f, v) * ufl.dx + ufl.inner(g, v) * ufl.ds

    options = {
               "ksp_type": "gmres",
               "pc_type": "hypre",
               "ksp_rtol": 1.0e-12
               }

    model = fem.petsc.LinearProblem(a, L, bcs=[left_bc, right_bc], petsc_options=options)
    logger.debug('Solving problem..')
    uh = model.solve()
    
    # Save solution in XDMF format
    with io.XDMFFile(comm, "u.xdmf", "w") as outfile:
        outfile.write_mesh(domain)
        outfile.write_function(uh)

    # # Update ghost entries and plot
    uh.vector.ghostUpdate(addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD)
    logger.debug("Post-process calculations")
    # Post-processing: Compute derivatives
    W = fem.VectorFunctionSpace(domain, ("Lagrange", 1))
    grad_expr = fem.Expression(-kappa * ufl.grad(uh), W.element.interpolation_points())
    grad_h = fem.Function(W)
    cdf_fun = fem.Function(W)
    cdf_fun2 = fem.Function(W)
    grad_h.interpolate(grad_expr)

    with io.XDMFFile(comm, "gradu.xdmf", "w") as file:
        file.write_mesh(domain)
        file.write_function(grad_h)

    logger.debug("Post-process Results Summary")
    # distribution at surfaces
    EPS = 1e-30
    def check_condition(v1, check_value=1):
        v2 = lambda x: check_value * (x[0] + EPS) / (x[0] + EPS)
        cdf_fun.interpolate(v2)
        return ufl.conditional(ufl.le(v1, cdf_fun), v1, cdf_fun)
    tol = 0.5
    new_v = fem.Expression(check_condition(ufl.inner(grad_h, n), tol), W.element.interpolation_points())
    cdf_fun2.interpolate(new_v)
    # lpvalue = fem.assemble_scalar(fem.form(ufl.inner(cdf_fun2, n) * ds(lmarker)))
    # rpvalue = fem.assemble_scalar(fem.form(ufl.inner(cdf_fun2, n) * ds(rmarker)))
    # logger.info((lpvalue, rpvalue))

I appreciate any tips for this, including pointing me to any of the tutorials that use ufl conditionals on vector function space. I found this example that uses non-VectorFunctionSpace: Ufl.conditional doesn't return precise results when using ('CG', 1)

2

The error you get is:

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/root/shared/mwe_cond.py", line 99, in <module>
    new_v = fem.Expression(check_condition(ufl.inner(grad_h, n), tol), W.element.interpolation_points())
  File "/root/shared/mwe_cond.py", line 96, in check_condition
    cdf_fun.interpolate(v2)
  File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 377, in interpolate
    self._cpp_object.interpolate(np.asarray(u(x), dtype=self.dtype), cells)
RuntimeError: Interpolation data has the wrong shape/size.

which is quite clear as:
v2 is a scalar function

while cfd_fun is in a vector space.

Changing the space of v2 and cfd_fun to V resolves that issue.

However, Expression does not support FacetNormal, as it is not a quantity defined in the interior of a cell (which is required for computing a function at the interpolation points).

Your integral (a scalar value) can be computed as follows:

    print(fem.assemble_scalar(fem.form(check_condition(ufl.inner(grad_h, n), tol)*ufl.ds)))

for the code above I use the whole boundary, but you can use a ds-measure with markers.

3

Thank you. That works!