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).
Can someone please help me in finding what am I doing wrong?
