# Issues with marking subdomains for 3D geometry

Hi all,
I am new to FEA in general and have been working on a steady state heating model for a 3D cylinder. For this case, I have created a simple cylinder mesh, with radius R and thickness H, defined the boundaries and then tried to use MeshFunction to mark the boundaries, but I cannot seem to mark the barrel of the cylinder at all.

My code is shown below.

from fenics import *
import matplotlib.pyplot as plt
from dolfin import *
import numpy as np

H = 0.06
R = 450e-3/2
cylinder = Cylinder(Point(0,0,0),Point(0,0,H),R,R)
geometry = cylinder

msh = generate_mesh(geometry,42)
V = FunctionSpace(msh,'P',2)
T = TestFunction(V)      # v
dT = TrialFunction(V)    # u

tol = 1e-14

class All(SubDomain):
def inside(self, x, on_boundary):
return on_boundary

class Barrel(SubDomain):
tol = 1e-14
def inside(self,x,on_boundary):
return np.isclose(x[0] ** 2 + x[1] ** 2, R**2) and on_boundary

# Along HR surface (z=0)
class HR(SubDomain):
tol = 1e-14
def inside(self,x,on_boundary):
return near(x[2],0,tol) and on_boundary

# Along AR surface (z=H)
class AR(SubDomain):
tol = 1e-14
def inside(self,x,on_boundary):
return near(x[2],H,tol) and on_boundary

# Marking boundaries
boundary_markers = MeshFunction("size_t", msh, msh.topology().dim()-1,0)
bAll = All()
bBar = Barrel()
bHR = HR()
bAR = AR()
bAll.mark(boundary_markers, 0)
bBar.mark(boundary_markers, 1)
bHR.mark(boundary_markers, 2)
bAR.mark(boundary_markers, 3)
File('Domain.pvd') << boundary_markers

From my understand, the first face (z=0) should be marked 2, the second 3 and the barrel 1. But my barrel stays 0 (or stays whatever the marking of all facets in the mesh are). I am not too sure what I am doing wrong here. Any help would be so great!

I have tried both of the solutions to those problems. Increasing the tolerance for near does not help, and neither does further refining my mesh.

Using this, I can mark the barrel easily.

boundary_markers = MeshFunction("size_t", msh, msh.topology().dim()-1)
CompiledSubDomain("near(x[2], 0.06, tol) && on_boundary", tol=1e-14).mark(boundary_markers,1)

I’m a bit confused. Seems like you are marking the end face as 1, and everything else as 0. Not the barrel only as 0.

Can you mark the area which is the barrel? I can’t find it.

The barrel is marked by the following in my code.

class Barrel(SubDomain):
tol = 1e-14
def inside(self,x,on_boundary):
return np.isclose(x[0] ** 2 + x[1] ** 2, R**2) and on_boundary

The barrel is the remaining boundary of the cylinder minus the two end faces.