Hello all,

I am currently solving a problem in which I need to apply a boundary condition on only particular parts of a boundary face. I have explained the details below, with the code that I’m working with.

The domain is a unit cube ranging from x = 0 to 1, y = 0 to 1 and z = -1 to 0.

I have a function that is defined over the domain. I have assigned it to random values as shown in the code below.

Now, based on the values of this function on the Top Boundary Face (z=0), I need to appropriately mark this boundary face.

That is, if this function is less than zero on the Top boundary, I want to mark those parts of the Top Boundary as 2.

If this function is greater than or equal to zero on the Top boundary, then I want to mark those parts of the Top Boundary as 1.

(For my specific case, this function needs to be defined over the entire domain)

Please see the code below.

```
from dolfin import *
import numpy as np
N = 30
mesh = UnitCubeMesh.create(N, N, N//2, CellType.Type.hexahedron)
# mesh is mapped to a [-1;0] domain along z
mesh.coordinates()[:, 2] = -mesh.coordinates()[:, 2]
# Top boundary: z = 0
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[2], 0.) and on_boundary
# exterior facets MeshFunction
facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
Top().mark(facets, 1) # Marking the entire top boundary as 1 initially
ds = Measure('ds', subdomain_data=facets)
# Function based on which the Top boundary (z=0) needs to be marked
V2 = FunctionSpace(mesh, "CG", 1)
indicator = Function(V2)
dim = V2.dim()
N = mesh.geometry().dim()
coor = V2.tabulate_dof_coordinates().reshape(dim,N)
fx_dofs = V2.dofmap().dofs()
x = coor[:, 0] # x-coordiantes
y = coor[:, 1] # y-coordinates
fx_x, fx_y = x[fx_dofs], y[fx_dofs] # x, y of components
# there are a few operations I need to do based on the x any y co-ordinates,
# but I have simplied it to a random function (as shown below) for now
fx = np.random.randn(*fx_x.shape)
# Final indicator function
indicator.vector()[fx_dofs] = fx
print(fx)
```

I want to do something like this,

if(indicator function on top boundary >= 0) --> mark those co-ordinates on the Top boundary as 1.

if(indicator function on top boundary < 0) --> mark those co-ordinates on the Top boundary as 2.

I am then hoping to use the marked regions of the boundary - ds(1) and ds(2) in the Variational Formulation.

Thank you all in advance for taking time to help me with the issue !