# Marking Subdomains in a cylindrical 3D mesh

**URL:** <https://fenicsproject.discourse.group/t/marking-subdomains-in-a-cylindrical-3d-mesh/3460>\
**Category:** mesh\
**Created:** [June 7, 2020, 2:48pm UTC](https://fenicsproject.discourse.group/t/marking-subdomains-in-a-cylindrical-3d-mesh/3460 "2020-06-07T14:48:20Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Sambeet\_Panigrahi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sambeet_panigrahi/32/1259_2.png) [@Sambeet\_Panigrahi](https://fenicsproject.discourse.group/u/Sambeet_Panigrahi)\
**Post date:** [June 7, 2020, 2:48pm UTC](https://fenicsproject.discourse.group/t/marking-subdomains-in-a-cylindrical-3d-mesh/3460/1 "2020-06-07T14:48:20Z")

</div>

I have a cylindrical 3D tetrahedral mesh (of radius = 5,height=10) which I imported from gmsh using meshio.Now I want to mark the subdomain another cylinder of same height but radius =1. I then use the following code to accomplish this.

```
from fenics import *
from mshr import *
from math import *
from dolfin import *
# Create Geometry
import meshio
msh = meshio.read("Hellogmsh.msh")
for cell in msh.cells:
    if cell.type == "triangle":
        triangle_cells = cell.data
    elif cell.type == "tetra":
        tetra_cells = cell.data

for key in msh.cell_data_dict["gmsh:physical"].keys():
    if key == "triangle":
        triangle_data = msh.cell_data_dict["gmsh:physical"][key]
    elif key == "tetra":
        tetra_data = msh.cell_data_dict["gmsh:physical"][key]
tetra_mesh = meshio.Mesh(points=msh.points, cells={"tetra": tetra_cells})
triangle_mesh =meshio.Mesh(points=msh.points,
                           cells=[("triangle", triangle_cells)],
                           cell_data={"name_to_read":[triangle_data]})
meshio.write("mesh.xdmf", tetra_mesh)

meshio.write("mf.xdmf", triangle_mesh)

# Subdomains

def boundary(x,on_boundary):
    return on_boundary

import numpy as np
tol = 1E-14

class Copper(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(x[0]*x[0] + x[1]*x[1] ) <= 1 + tol

    
class Air(SubDomain):
    def inside(self, x, on_boundary):
        return np.sqrt(x[0]*x[0] + x[1]*x[1] ) >= 1 - tol

# Read the mesh
mesh = Mesh ( ) # Create mesh object
with XDMFFile("mesh.xdmf") as infile:
    infile.read(mesh)
mvc = MeshValueCollection( "size_t",mesh,3)
info(mesh) # Print some useful mesh info

materials = MeshFunction('size_t',mesh,mesh.topology().dim()-1)
Subdom_cu = Copper()
Subdom_air = Air()
materials.set_all(0)
Subdom_cu.mark(materials, 1)

file = File("/home/fenics/shared/Cylinder/Materialsgmsh.pvd")
file << materials

```

The radius of the cylinder = sqrt(x^2 + y^2) \< 1.That is the logic here.But when I see this in paraview, instead of marking an inner smaller cylinder, its marking just a line(very faint red line in the pcture).

 ![snap_paraview](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/b400ae7e6024c45dcf1afb608426c7499d371dfe.png)

Can someone please help me in finding what am I doing wrong?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [June 7, 2020, 2:51pm UTC](https://fenicsproject.discourse.group/t/marking-subdomains-in-a-cylindrical-3d-mesh/3460/2 "2020-06-07T14:51:53Z")

</div>

You should refine your mesh, as you only have about 10 elements in diameter, way too few to resolve a cylinder of radius one.

---

<div class="post-metadata">

**Author:** ![Sambeet\_Panigrahi](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/sambeet_panigrahi/32/1259_2.png) [@Sambeet\_Panigrahi](https://fenicsproject.discourse.group/u/Sambeet_Panigrahi)\
**Post date:** [June 7, 2020, 3:01pm UTC](https://fenicsproject.discourse.group/t/marking-subdomains-in-a-cylindrical-3d-mesh/3460/3 "2020-06-07T15:01:39Z")

</div>

Thanks that solved it! I used a radius of 4 instead (flexible constraint in my problem)and I get the cylinder.

 ![snap_2223](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/bce882664cafb414f4622dbfaf39062b9ecdac5d.png)