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