Plotting subdomains/internal boundaries

mesh
1

Hello,
I am working with subdomains and I need to identify the internal boundary (interface between both subdomains) for integration. For the subdomains I am using a mesh function named “domains” and by running plot (domains) i can successfully see that the subdomains are marked correctly. Then, I am also using a mesh function to mark the interface between them, but I’m having trouble plotting this to make sure that I am doing it correctly. Could anyone help me on what is the easiest way to visualize this internal boundary? Thank you! This is how I am defining the mesh functions:

class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
omega = Omega() 
domains.set_all(0) 
omega.mark(domains, 1)  
dx = Measure('dx')(subdomain_data = domains)

and ploting domains gives me for example:

image

Then, for the boundary this is what I am doing:

class IntBound(SubDomain):
    def inside(self, x, on_boundary):
        return near(phi(x),.0) 

intBound = IntBound() 
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
int_boundary.set_all(0)
intBound.mark(int_boundary, 1)
ds_int = Measure("ds")(subdomain_data=int_boundary) #ds(1) integrates on internal boundary

But when I try plotting int_boundary I get

image

I am not sure if I am not marking correctly of if is a plot problem.

2

There are several issues with your code. The first one is that phi does not align with the grid, thus using the following code:

will not mark any facets.
This can be verified with the following code snippet (note usage of "dS" as these are interior facets):

intBound = IntBound() 
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
intBound.mark(int_boundary, 1)
dS_int = Measure("dS", domain=mesh, subdomain_data=int_boundary) 

print(assemble(Constant(1)*dS_int(0)), assemble(Constant(1)*dS_int(1)))

The second issue is that the dolfin interface to matplotlib does not support plotting facet markers and the following error should be printed in your terminal:

*** Warning: Matplotlib plotting backend does not support mesh function of dim 1. Continuing without plotting...

This can be bypassed by saving to pvd and visualizing in Paraview.

File("facets.pvd") << int_boundary

Next I present a code that resolves your problem:

from dolfin import *
import numpy as np

def phi(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - 0.2**2

class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0


mesh = UnitSquareMesh(30, 30)
domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

omega = Omega() 
omega.mark(domains, 1)  
dx_int = Measure("dx", domain=mesh, subdomain_data=domains)

print(assemble(Constant(1)*dx_int(0)), assemble(Constant(1)*dx_int(1)))
File("domains.pvd") << domains


# Interior facets between the two domains
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
tdim = mesh.topology().dim()
mesh.init(tdim-1, tdim)
facet_to_cell = mesh.topology()(tdim-1, tdim)
num_facets = facet_to_cell.size()
domain_values = domains.array()
facet_values = int_boundary.array()
for facet in range(num_facets):
    # Check if interior
    cells = facet_to_cell(facet) 
    if len(cells == 2):
        # Check if facet is on the boundary between two cells with different markers
        values = np.unique(domain_values[cells])
        if len(values) == 2:
            facet_values[facet] = 1        
    else:
        continue

dS_int = Measure("dS", domain=mesh, subdomain_data=int_boundary) 

print(assemble(Constant(1)*dS_int(0)), assemble(Constant(1)*dS_int(1)))
File("facets.pvd") << int_boundary
2 Likes
3

Thank you very much! That was very helpful!

4

Hello @dokken. Thank you again for your help for identifying the internal boundary. I am now using this internal boundary for integration, but I’m having trouble on what is the most efficient way to do that. I made a new post with this question and I wonder if you maybe could help me! This is the link Efficient integration on internal boundary

5

I am on vacation for the next 5 days, and do not have opportunity to answer posts in that time.

6

Thats ok, thank you for replying anyway! :slight_smile:

7

Hi dokken,

Apologies for bumping an older thread,
However, this code you have written here for identifying the subdomain boundary works well in series, but due to the way facet ownership works in parallel it seems to also mark the process boundaries when using MPI.

How could one extend this to work in parallel?

8

There seems to be a bug when using subdomains to mark cells in parallel with shared_vertex/shared_facet.
You can submit a bug at the bitbucket page. However, it is not likely to be fixed as DOLFIN is no longer developed.

A workaround would be to write the tags to XDMF in serial, and then read then in again when running in parallel.

In the development version of dolfin, dolfinx this works out of the box:

import numpy as np
import dolfinx
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 25, 25,
                                       ghost_mode=dolfinx.mesh.GhostMode.shared_facet)


def Omega(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - 0.2**2 > 0


tdim = mesh.topology.dim
indices = dolfinx.mesh.locate_entities(mesh, tdim, Omega)

num_cells = mesh.topology.index_map(tdim).size_local + \
    mesh.topology.index_map(tdim).num_ghosts
local_cell_marker = np.zeros(num_cells, dtype=np.int32)
local_cell_marker[indices] = 1
ct = dolfinx.mesh.MeshTags(mesh, tdim, np.arange(
    num_cells, dtype=np.int32), local_cell_marker)


mesh.topology.create_connectivity(tdim-1, tdim)
f_to_c = mesh.topology.connectivity(tdim-1, tdim)


num_facets = mesh.topology.index_map(tdim-1).size_local + \
    mesh.topology.index_map(tdim-1).num_ghosts

facets = []

for i in range(num_facets):
    cells = f_to_c.links(i)
    if len(cells) == 2:
        values = np.unique(local_cell_marker[cells])
        if len(values) == 2:
            facets.append(i)

all_facets = np.zeros(num_facets, dtype=np.int32)
all_facets[facets] = 1
ft = dolfinx.mesh.MeshTags(mesh, tdim-1, np.arange(num_facets, dtype=np.int32),
                           all_facets)

with dolfinx.io.XDMFFile(mesh.comm, "mf.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)
    ct.name = "Cell tag"
    xdmf.write_meshtags(ct)
    ft.name = "Facet tag"
    xdmf.write_meshtags(ft)


print(MPI.COMM_WORLD.rank, facets)
1 Like