Subdomain based on cell index

1

Hi,
I want to do some post-processing (calculating an integration) on some elements that are intersected by a specified curved line. Is it possible to loop over all cells and based on some conditions mark some cells and make those marked cells as a subdomain for later use? I am wondering if you could provide some hints regarding the matter.
Thanks

2

I would suggest having a look at: Marking subdomains of a mesh — DOLFIN documentation

1 Like
3

Hi @dokken, If intend to print out or basically just know all the cells in a mesh associated with (or fall within) a marked subdomain how do I go about it? A post I saw here was about 7 years ago and the suggested solution doesn’t work anymore.

Thanks.

4

You first need to specify what version of FEniCS you are running (legacy FEniCS, FEniCSx).
it would also help to have an actual example of how you create your subdomain.

5

I’m using the Legacy version of fenics. Here is how I created my subdomain.

Assuming I already defined the mesh, the finite element, and my function spaces

Here’s a short version

nx = ny = 256
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)
materials = MeshFunction("size_t", mesh, mesh.topology().dim()) # over cells f

# iterate on the cells of the mesh
for cell_no in range(len(materials.array())):

	# define the cell
	cell = MeshEntity(mesh, 2, cell_no)

	# define the barycenter of the cell
	cell_bc = cell.midpoint()

	# create subdomain class of radius r for each cell on the mesh
	class Interaction_range(SubDomain):
		def inside(self, x, on_boundary):
			return on_boundary and (x[0]-cell_bc[0])*(x[0]-cell_bc[0]) + (x[1]-cell_bc[1])*(x[1]-cell_bc[1]) <= r*r + tol

	# mark the interaction range
	subdomain_1 = Interaction_range()
	subdomain_1.mark(materials, 1)
     
    # integrate on the subdomain
    nonlocal_interaction = chi_f*(0.5*dot(grad(f), g)*dx(1)
    
    
	
	# remark as zero after each iteration
	subdomain_1.mark(materials, 0)

A broader version is below


from __future__ import print_function
from fenics import *

nx = ny = 256
mesh = RectangleMesh(Point(-1, -1), Point(1, 1), nx, ny)

# Define function space for system of concentrations
P2 = FiniteElement('P', triangle, 2)
P0 = FiniteElement('DG', triangle, 0)

element = MixedElement([P1, P1, P1, P1])
element2 = MixedElement([P2, P2, P2, P2])
V = FunctionSpace(mesh, element2)
V1 = FunctionSpace(mesh, P1)
V0 = FunctionSpace(mesh, P0)
V2 = FunctionSpace(mesh, P2)

# Define test functions
v_g, v_f = TestFunctions(V) 

# Define functions
u = Function(V)

# initial condition
u_0 = Expression(('0.1',
                  '0.4*((0.5+0.5*tanh(20*x[0]-3))+(0.5+0.5*tanh(-20*x[0]-4)))'), degree=2)

u_n = interpolate(u_0, V)

# Split system functions to access components
g, f = split(u)

func_g = dot(grad(g), grad(v_g))*dx +g*v_g*dx 
func_f = dot(grad(f), grad(v_f))*dx + f*v_f*ds

materials = MeshFunction("size_t", mesh, mesh.topology().dim()) # over cells f
bd = MeshFunction("size_t", mesh, mesh.topology().dim()-1) # over facets for boundary
ds = Measure('ds', domain=mesh, subdomain_data=bd)
dx = Measure('dx', domain=mesh, subdomain_data = materials)

# define DG order 0 element, one value per cell since we want to assign a constant value to each cell
DG_Function = Function(V0)

# range of interaction kernel
r = 0.0079 

tol = 1E-14

# iterate on the cells of the mesh
for cell_no in range(len(materials.array())):

	# define the cell
	cell = MeshEntity(mesh, 2, cell_no)

	# define the barycenter of the cell
	cell_bc = cell.midpoint()

	# create subdomain class of radius r for each cell on the mesh
	class Interaction_range(SubDomain):
		def inside(self, x, on_boundary):
			return on_boundary and (x[0]-cell_bc[0])*(x[0]-cell_bc[0]) + (x[1]-cell_bc[1])*(x[1]-cell_bc[1]) <= r*r + tol

	# mark the interaction range
	subdomain_1 = Interaction_range()
	subdomain_1.mark(materials, 1)
    
     # integrate on the subdomain
    nonlocal_interaction = chi_f*(0.5*dot(grad(f), g)*dx(1)

    # append each integration result (a float) to the corresponding F
    DG_Function.vector()[cell_no]=assemble(nonlocal_interaction)
	
	# remark as zero after each iteration
	subdomain_1.mark(materials, 0)

# project it to a V1 space so we can add it to the rest and the overall shape function will be obtained
Nonlocal_interaction=interpolate(DG_Function, V1)

# Add to the rest of the weak form
func_f+=Nonlocal_interaction*v_f*dx

# Adding the functions together
f_u = func_g + func_f

Thank you in anticipation

6

You already have all the information you need in the code you have here.
materials.array() is an array of size number of cells on your process.
You can use numpy.flatnonzero(materials.array() == some_tag) to get the cell indices of those cells that are marked with a given tag.

7

Thank you very much @dokken

8

This works perfectly. Is there perhaps also a way to get the indices of the elements in the subdomain marked?

9

Im not sure what you mean by:

10

My apologies. You already answered me. I didn’t know numpy.flatnonzero(materials.array() == some_tag) returns the indices of the finite elements in the marked subdomain, I assumed it returned the indices of the rectangles in the fine mesh marked by the subdomain.

Thanks a lot @dokken

1 Like