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

variational formulation
1

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 !

2

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.

3

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.

4

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[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)
2 Likes
5

Dear @dokken,

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

6

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.

7

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
8

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!