"Remarking" boundary facets (changing markers)

mesh
1

Dear all,
With the amazing functionality offered by the LEoPart library I am able to track “strain-history” of particles inside a complex domain.
image

In the picture above, one can distinguish the particles inside the red box that will leave through Outlet1 and the particles inside the blue box leaving through Outlet2. These particles will carry information of their own “strain-history”.

For my application I am interested in the cumulative strain-history of particles leaving a specific boundary. For example, I would like to know the average of strain-histories of particles leaving through outlet1 and compare the average strain-histories of particles leaving through outlet2. With the effect of the stenosis, I expect higher strain-histories carried by particles inside the red box.

LEoPart functionality enables specification of “open” boundaries, which means particles will be deleted if they “escaped” through a facet on an open boundary. To this end, one should specify an “open boundary” by marking it with integer “2” in a MeshFunction object.

Currently, I have a model with 1 inlet, 2 outlets and a wall. From a .h5-file I obtain a MeshFunction object carrying the markers for these boundaries: the inlet is marked as “1”, outlet 1 as “2”, outlet 3 as “3” and wall is marked as “4”. I use these markers to assign boundary conditions and calculate flow at the specific outlets, for example.

To enable LEoPart “open” boundary functionality, I would require a second MeshFunction object with different markers.
Is it possible to create a second meshfunction object and to reassign markers conveniently?
For my application I would need to reassign facets carrying markers “2” and “3”, so that they both carry marker “2” to become “open” in LEoPart.

With kind regards,

2

In addition to the above, I found a MWE example for a similar problem in another post:

from fenics import *

mesh = UnitSquareMesh(3,3)

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0, DOLFIN_EPS) and on_boundary
    
class Rest(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
    
domains     = MeshFunction('size_t', mesh, mesh.topology().dim(), 0) 
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
new_boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1, 0)
Rest().mark(domains,2) ; Rest().mark(boundaries,2)
Left().mark(domains,1) ; Left().mark(boundaries,1)
Rest().mark(domains,2) ; Rest().mark(new_boundaries,2)
Left().mark(domains,1) ; Left().mark(new_boundaries,1)

for i in range(len(new_boundaries)):
    if new_boundaries==1:
        new_boundaries[i]=3

File('boundaries.pvd') << boundaries # save PVD for visualization

The only difference with my case is that I do not have class definitions for the boundaries of interest, as they are “too” complicated for capture. Instead, I use a hdf5 file from which I read boundary markers:

mesh = Mesh()
hdf = HDF5File(mesh.mpi_comm(),mesh_file, "r")  

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)

hdf.read(boundaries, "/boundaries")

new_boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
hdf.read(new_boundaries, "/boundaries")

For my application, I can do:

for i in range(len(new_boundaries)):
    if new_boundaries[i]==3:
        new_boundaries[i]=2

I tried to implement this for a complex mesh with LEoPart.
This works and the MeshFunction correctly behaves with LEoPart “open” or “closed” boundary conditions.

3

I can now confirm that the above works.