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

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

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.

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.

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

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.

Thank you very much @dokken

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

Im not sure what you mean by:

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