# Marking the Parts of a Boundary Face based on Function Values defined over the Mesh

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)

``````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, 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 !

You need to clarify a few things. Say that we have a one of the top facets, which would be a quadrilateral.
This consists of four dofs, one at each vertex.
What do you want to happen if 1 or 2 of the dofs has a value >=0 and the others have value < 0?

Please also note that the support for hexahedral elements are very limited in dolfin, and I would advise you to use dolfin-x for such problems.

Dear @dokken,

Thank you so much for the reply and clarification. I am now able to understand the picture.

My understanding based on your explanation is that we cannot mark individual co-ordintaes on the top boundary for using in the variational formulation; we can only mark the facets - and each facet has four dofs. Is my understanding correct?

I have initially marked the entire top boundary as 1.

I would like to now do the following: if the function is >= 0 at more than 2 dofs of a top boundary facet, then I want mark those top boundary facets as 2. Can you please help me with a code for this?

Here is an implementation of such a function.
In my code, I only add the facet if there are 2 or more markers satisfying my constraint:

``````from dolfin import *
import numpy as np

class Topboundary(SubDomain):
def inside(self, x, on_boundary):
return x > 1-DOLFIN_EPS

mesh = UnitCubeMesh(3,3,3)
tdim = mesh.topology().dim()
fdim = tdim - 1
mesh.init(tdim, fdim)
mesh.init(fdim, tdim)

marker = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
Topboundary().mark(marker, 1)

c_2_f = mesh.topology()(tdim, fdim)
f_2_c = mesh.topology()(fdim, tdim)

V = FunctionSpace(mesh, "CG", 1)
u = interpolate(Expression("x", degree=1), V)

facet_indices = marker.where_equal(1)

petsc_vec = u.vector().vec()
marked_facets = []
with petsc_vec.localForm() as loc:
petsc_array = loc.array
for facet in facet_indices:
cells = f_2_c(facet)
assert(len(cells)==1)
cell = cells
facets = c_2_f(cell)
local_index = np.argwhere(facets==facet)[0,0]
local_dof_indices = V.dofmap().entity_closure_dofs(mesh, fdim, [local_index])
facet_dofs = V.dofmap().cell_dofs(cell)
num_marked =  0
values = petsc_array[facet_dofs]
for value in values:
if value > 0.5:
num_marked += 1

if num_marked > 1:
marked_facets.append(facet)
marker.array()[marked_facets] = 2
File("marker.pvd") << marker
print(marked_facets)
``````
2 Likes

Dear @dokken,

Thank you so much for sharing the code and helping me with the problem !

Hello dokken!

I am trying to understand your code. Could you tell me the difference between `entity_closure_dofs` and `entity_dofs`. `entity_closure_dofs` worked here but `entity_dofs` failed. Besides, I found that `num_entity_dofs(0)` returns 2 and `num_entity_dofs(dim)` with dim=1,2 returns 0 for 2D mesh. `num_entity_closure_dofs` works. I have read the source code, still confusing about these functions.

Entity_dofs lists the degrees of freedom that are associated with the entity.
See for instance: DefElement
or
DefElement
for examples where the associated entity for each basis function is listed.

The entity closure dofs loops through all entities of a lower dimension (for a facet it loops through vertices, edges and facets) and finds all dofs associated with any of these entities.

For a cell it loops through (vertices, edges, facets and the cell itself).
For an edge in 3D it would only loop through (vertices and the edge itself).

1 Like

Thank you very much!

I undstand it. In 3D, `num_entity_dofs` returns [1,0,0,0] for P1 scalar element and it returns [1,1,0,0] for P2 scalar element.`num_entity_closure_dofs` returns [1,2,3,4] for P1 scalar element and it returns [1,3,6,10] for P2 scalar element.

Thanks again!