Compute minimum and maximum values of expression over boundary

Legacy DOLFIN
1

Hello everyone,

I was trying to define an expression that contains the normal n = FacetNormal(mesh). I need to compute the minimum and maximum values of the expression over a boundary region. I have tried with the following code:

def compute_SSF(n):
    # define SSF
    SSF_expr = n[0] + n[1]

    # create SSF vector from expression
    SSF_function_space = FunctionSpace(mesh, "DG", 0)
    SSF_fun = Function(SSF_function_space)
    SSF_fun.interpolate(SSF_expr)

    SSF_local = SSF_fun.vector().get_local()

    SSF_min = SSF_local.min()
    SSF_max = SSF_local.max()

    # collect from MPI processes
    SSF_min = comm.allreduce(SSF_min, op=MPI.MIN)
    SSF_max = comm.allreduce(SSF_max, op=MPI.MAX)

    SSF_rescaled_expr = (SSF_expr - SSF_min) / (SSF_max - SSF_min)

    return SSF_rescaled_expr

n = FacetNormal(mesh)
SSF = compute_SSF(n)

but it raises the following error

AttributeError: 'Sum' object has no attribute '_cpp_object'

Is there a way to evaluate the expression over a generic boundary ds(tag) to get the values? Thank you in advance!

2

I thought I could get this easily done with scifem (v0.19.0):

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 13)

Vh = dolfinx.fem.functionspace(mesh, ("Lagrange", 2))
uh = dolfinx.fem.Function(Vh)
uh.interpolate(lambda x: x[0] + x[1] * (1 - x[1]))

tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

nh = ufl.FacetNormal(mesh)
expr = uh * nh


# Define submesh of the facets you would like to compute the maximum over
facet_submesh, entity_map = dolfinx.mesh.create_submesh(mesh, fdim, exterior_facets)[0:2]
num_submesh_cells = facet_submesh.topology.index_map(fdim).size_local
parent_facets = entity_map.sub_topology_to_topology(
    np.arange(num_submesh_cells, dtype=np.int64), inverse=False
)


# Transfer data to an appropriate function space on the facet submesh
Vh_facet = dolfinx.fem.functionspace(facet_submesh, ("DG", 2, (mesh.geometry.dim,)))
exterior_facet_as_integration_entities = dolfinx.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.exterior_facet, mesh.topology, parent_facets
).reshape(-1, 2)
facet_expr = dolfinx.fem.Expression(expr, Vh_facet.element.interpolation_points)
facet_values = facet_expr.eval(mesh, exterior_facet_as_integration_entities)
expr_func = dolfinx.fem.Function(Vh_facet)
if hasattr(expr_func._cpp_object, "interpolate_f"):
    expr_func._cpp_object.interpolate_f(
        facet_values.reshape(-1, 2).T.copy(), np.arange(num_submesh_cells, dtype=np.int32)
    )
else:
    expr_func._cpp_object.interpolate(
        facet_values.reshape(-1, 2).T.copy(), np.arange(num_submesh_cells, dtype=np.int32)
    )

# Use scifem to compute maximum
import scifem

value, point_of_extrema = scifem.compute_extrema(ufl.dot(expr_func, expr_func), kind=max)
print(value, point_of_extrema)
with dolfinx.io.VTXWriter(facet_submesh.comm, "facet_expr.bp", [expr_func]) as writer:
    writer.write(0.0)

returning

1.5625 [1.  0.5]

However, due to the fact these degrees of freedom being DG, there seems to be a permutation missing, so I would need to think about this a bit longer.
I think one can get it right by getting the permutations for the same CG space, but I would have to test this. You can spot the issue on the top surface here, where it should be linearly increasing.

image
image807×671 29.8 KB

3

Here is a bit of a complicated workaround, where one transfer the function itself, and the facet normal individually to the facet submesh before performing the optimization. Here tested with a vector space to ensure that components are mapped correctly

from mpi4py import MPI
import dolfinx
import ufl
import numpy as np
import scifem


def cell_locator(x):
    return x[0] > 0.5 - 1e-10


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 13)

degree = 3
Vh = dolfinx.fem.functionspace(mesh, ("Lagrange", degree, (mesh.geometry.dim,)))
uh = dolfinx.fem.Function(Vh)
uh.interpolate(
    lambda x: (x[1] * (1 - x[1]), x[0] + (1 - x[1]) * x[1]),
    cells0=dolfinx.mesh.locate_entities(mesh, mesh.topology.dim, cell_locator),
)

tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_connectivity(fdim, tdim)
exterior_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)

nh = ufl.FacetNormal(mesh)


# Define submesh of the facets you would like to compute the maximum over
facet_submesh, entity_map = dolfinx.mesh.create_submesh(mesh, fdim, exterior_facets)[0:2]
num_submesh_cells = facet_submesh.topology.index_map(fdim).size_local
parent_facets = entity_map.sub_topology_to_topology(
    np.arange(num_submesh_cells, dtype=np.int64), inverse=False
)
exterior_facet_as_integration_entities = dolfinx.fem.compute_integration_domains(
    dolfinx.fem.IntegralType.exterior_facet, mesh.topology, parent_facets
).reshape(-1, 2)


# Transfer uh to the facet submesh
V_facet = dolfinx.fem.functionspace(facet_submesh, ("DG", degree, (mesh.geometry.dim,)))
u_facet = dolfinx.fem.Function(V_facet)

# WORKAROUND: Use equivalent continuous space to get permutation on facets
facet_expr = dolfinx.fem.Expression(uh, V_facet.element.interpolation_points)
facet_values = facet_expr.eval(mesh, exterior_facet_as_integration_entities)

tmp_space = dolfinx.fem.functionspace(mesh, ("Lagrange", degree, (mesh.geometry.dim,)))
cell_info = mesh.topology.get_cell_permutation_info()
ft = facet_submesh.basix_cell()
for i in range(exterior_facet_as_integration_entities.shape[0]):
    perm = np.arange(facet_values.shape[1], dtype=np.int32)
    tmp_space.element.basix_element.permute_subentity_closure_inv(
        perm,
        cell_info[exterior_facet_as_integration_entities[i, 0]],
        ft,
        int(exterior_facet_as_integration_entities[i, 1]),
    )
    facet_values[i] = facet_values[i][perm]


facet_uh = dolfinx.fem.Function(V_facet)
if hasattr(facet_uh._cpp_object, "interpolate_f"):
    facet_uh._cpp_object.interpolate_f(
        facet_values.reshape(-1, mesh.geometry.dim).T.copy(),
        np.arange(num_submesh_cells, dtype=np.int32),
    )
else:
    facet_uh._cpp_object.interpolate(
        facet_values.reshape(-1, mesh.geometry.dim).T.copy(),
        np.arange(num_submesh_cells, dtype=np.int32),
    )

# Transfer normal vector onto facet submesh
N = dolfinx.fem.functionspace(facet_submesh, ("DG", 0, (facet_submesh.geometry.dim,)))
nh_facet = dolfinx.fem.Function(N, name="FacetNormal")
n_expr = dolfinx.fem.Expression(nh, N.element.interpolation_points)
n_values = n_expr.eval(mesh, exterior_facet_as_integration_entities)
if hasattr(nh_facet._cpp_object, "interpolate_f"):
    nh_facet._cpp_object.interpolate_f(
        n_values.reshape(-1, facet_submesh.geometry.dim).T.copy(),
        np.arange(num_submesh_cells, dtype=np.int32),
    )
else:
    nh_facet._cpp_object.interpolate(
        n_values.reshape(-1, facet_submesh.geometry.dim).T.copy(),
        np.arange(num_submesh_cells, dtype=np.int32),
    )

with dolfinx.io.VTXWriter(facet_submesh.comm, "facet_expr.bp", [facet_uh]) as writer:
    writer.write(0.0)
with dolfinx.io.VTXWriter(facet_submesh.comm, "facet_normal.bp", [nh_facet]) as writer:
    writer.write(0.0)


# Use scifem to compute maximum


expr = facet_uh[0] * nh_facet
max_expression = ufl.dot(expr, expr)

value, point_of_extrema = scifem.compute_extrema(max_expression, kind=max)
print(value, point_of_extrema)

Here you see the facet normal (in blue) and uh mapped onto the facets of the mesh:

image
image807×671 36.4 KB

We get the max value of the x component of uh as:

0.06250000000000003 [1.  0.5]

which is 0.25**2

4

Thank you. It was not as straightforward as I was expecting it to be. I am trying to replicate your procedure in legacy fenics now. Some commands differ but most of the structure should be the same. Thank very much for your help again!

5

This is due to the fact that the FacetNormal is not a numerical set of values, but a symbolic quantity computed at runtime within assembly/interpolation.
Furthermore, as you want the max on the boundary of the cell, and if your space is not piecewise linear, it is not sufficient to look at the maximal values at the degrees of freedom.

I’d say that doing the optimization process in legacy FEniCS is a bit more complicated, due to the fact it has limited support of evaluation of UFL expressions at runtime (without being part of an integral).

6

OK thank you. I have also the possibility to use a Function variable instead of FacetNormal, maybe that could simplify a bit the task