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

**URL:** <https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575>\
**Category:** variational formulation\
**Created:** [November 19, 2020, 5:35pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575 "2020-11-19T17:35:53Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sabarish\_Gopi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sabarish_gopi/32/980_2.png) [@Sabarish\_Gopi](https://fenicsproject.discourse.group/u/Sabarish_Gopi)\
**Post date:** [November 19, 2020, 5:35pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/1 "2020-11-19T17:35:53Z")

</div>

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.

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

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 19, 2020, 6:29pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/2 "2020-11-19T18:29:17Z")

</div>

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](https://github.com/FEniCS/dolfinx/) for such problems.

---

<div class="post-metadata">

**Author:** ![Sabarish\_Gopi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sabarish_gopi/32/980_2.png) [@Sabarish\_Gopi](https://fenicsproject.discourse.group/u/Sabarish_Gopi)\
**Post date:** [November 19, 2020, 7:36pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/3 "2020-11-19T19:36:52Z")

</div>

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?

Thank you in advance for your help.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [November 19, 2020, 8:38pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/4 "2020-11-19T20:38:22Z")

</div>

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:

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

class Topboundary(SubDomain):
    def inside(self, x, on_boundary):
        return x[2] > 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[0]", 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[0]
        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)

```

---

<div class="post-metadata">

**Author:** ![Sabarish\_Gopi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sabarish_gopi/32/980_2.png) [@Sabarish\_Gopi](https://fenicsproject.discourse.group/u/Sabarish_Gopi)\
**Post date:** [November 19, 2020, 9:18pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/5 "2020-11-19T21:18:42Z")

</div>

Dear @dokken,

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

---

<div class="post-metadata">

**Author:** ![shaoyaoqian](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/shaoyaoqian/32/4777_2.png) [@shaoyaoqian](https://fenicsproject.discourse.group/u/shaoyaoqian)\
**Post date:** [September 25, 2023, 1:11pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/6 "2023-09-25T13:11:37Z")

</div>

> [@dokken](#):
>
> `entity_closure_dofs`

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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [September 25, 2023, 2:25pm UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/7 "2023-09-25T14:25:37Z")

</div>

Entity\_dofs lists the degrees of freedom that are associated with the entity.  
See for instance: [DefElement](https://defelement.com/elements/examples/interval-lagrange-equispaced-3.html)  
or  
[DefElement](https://defelement.com/elements/examples/triangle-lagrange-equispaced-3.html)  
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).

---

<div class="post-metadata">

**Author:** ![shaoyaoqian](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/shaoyaoqian/32/4777_2.png) [@shaoyaoqian](https://fenicsproject.discourse.group/u/shaoyaoqian)\
**Post date:** [September 26, 2023, 2:24am UTC](https://fenicsproject.discourse.group/t/marking-the-parts-of-a-boundary-face-based-on-function-values-defined-over-the-mesh/4575/8 "2023-09-26T02:24:09Z")

</div>

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!