Assign boundary markers with an indicator function

[enricofacca]

Hello everyone.
Is there a combination of simple Firedrake functions to mark portions of the boundary of a mesh as the intersection of the boundary and the support indicator function? (a similar question was #2860)

In the following mwe, I tried to use the function RelabeledMesh to mark the portion of the boundary the upper-right corner of a square (x,y) = (1 x [0.5,1]) U ([0.5,1] x 1).

from firedrake import *

mesh = UnitSquareMesh(10, 10)
x, y = SpatialCoordinate(mesh)

# Define indicator function suitable for RelabeledMesh
W = FunctionSpace(mesh, "HDiv Trace", 0)
indicator = Function(W,name="boundary")
indicator.interpolate(conditional(And(x > 0.5, y > 0.5), 1, 0))
relabeled_mesh = RelabeledMesh(mesh, [indicator], [99])

# Verify with a Poisson Solve
V = FunctionSpace(relabeled_mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

a = dot(grad(u), grad(v)) * dx
L = Constant(1.0) * v * dx
bc = DirichletBC(V, 0.0, 99)

u_sol = Function(V,name="u_sol")

# Apply BC only on the newly marked region     
try:
    solve(a == L, u_sol, bcs=[bc])
    print("Solve successful using new boundary ID.")
except Exception as e:
    print(f"Solve failed: {e}")

outfile = VTKFile("solution.pvd")
outfile.write(u_sol)

However, zero Dirichlet boundary conditions are imposed on all nodes for x>0.5, y>0.5 including those inside the domain.

[connorjward]

What if you replace

And(x > 0.5, y > 0.5)

with

Or(And(x > 0.99, y > 0.5), And(x > 0.5, y > 0.99))

[connorjward]

How is the relevant part of the boundary identified? I think the usual approach is to mark it in the mesh generator (e.g. netgen, gmsh) and this will get passed through to the DMPlex.

[enricofacca]

I am using [this tools of pygalmesh .
I have to two numpy indicators. One for the domain and another describing a structure that intersects the boundary of the domain in the portion I am interested to.
I can generate the mesh I need from the first array, store it in a .msh file, and load it in Firedrake.

Your suggestion is to store the information I need using pygalmesh tools.

[connorjward]

Your suggestion is to store the information I need using pygalmesh tools.

Yes. Those numpy arrays will not make sense inside of Firedrake.

[enricofacca]

I was able to "interpolate" the second array using the instruction descibed here .
Thus I was trying to mark the boundary faces (and nodes).
So the question I asked.