Replicating mesh generated using Python script with the one in SV gui

[Jimmy Azarnoosh]

Hello,

We generated a mesh in SV gui as shown below that has a pretty nice element size around the stenosis.

fig1.JPG (47.11 KiB) Viewed 790 times

Parameters assigned for this mesh are:

fig2.JPG (27.64 KiB) Viewed 790 times

We tend to do the same using Python code. Our current code, attached here, generated the mesh that is not desirable.

fig3.JPG (64.44 KiB) Viewed 790 times

What changes are needed to make to the code to replicate the mesh from gui?

Thank you,
Jimmy

Code: Select all

def create_mesh(directory, model, centerlines, face_names, edge_size, boundary_layer=True):
  print("Generating mesh ...")
  mesher = sv.meshing.TetGen()

  if isinstance(model, type(sv.modeling.OpenCascade())):
    polydata = model.get_polydata()
    model = sv.modeling.PolyData()
    model.set_surface(polydata)
    face_ids = model.compute_boundary_faces(60.0)
  
  #caps = model.identify_caps()
  #print(caps)
  #wall_id = [i + 1 for i in range(len(caps)) if not caps[i]]
  inlet_id = [face_id for face_id, face_name in face_names.items() if 'in' in face_name]
  outlet_ids = [face_id for face_id, face_name in face_names.items() if 'in' not in face_name]
  outlet_names = {face_id:face_name for face_id, face_name in face_names.items() if face_name != 'in'}
  wall_id = [i for i in range(1, len(face_names) + 1) if i not in list(face_names.keys())]
  print("Face ids:", model.get_face_ids())
  print("Wall id {}".format(wall_id))
  print("Inlet id {}".format(inlet_id))
  print("Outlet ids & names {}".format(outlet_names))

  # Set the model for the mesher.
  mesher.set_model(model)
  mesher.set_walls(wall_id)

  # If the model has faces already computed (i.e. has 'ModelFaceID' array) then
  # don't call this, the face IDs will no longer match the original face IDs.
  #mesher.compute_model_boundary_faces(angle=60.0)
  face_ids = mesher.get_model_face_ids()
  #print(face_ids)
  # Set meshing options.
  #mesher.compute_model_boundary_faces(45)  
  #faces = mesher.get_face_polydata(0) 

  print("\nSet meshing options ... ")
  options = sv.meshing.TetGenOptions(global_edge_size=edge_size, surface_mesh_flag=True, volume_mesh_flag=True)
  options.radius_meshing_on = True
  options.boundary_layer_inside = True
  options.radius_meshing_scale = 0.4 

  if boundary_layer:
    mesher.set_boundary_layer_options(number_of_layers=3, edge_size_fraction=0.5, layer_decreasing_ratio=0.8, constant_thickness=False)
    #options.no_bisect = False

  # Generate the mesh. 
  mesher.generate_mesh(options)

  # Get the mesh as a vtkUnstructuredGrid. 
  mesh = mesher.get_mesh()
  print("Mesh:")
  print("  Number of nodes: {0:d}".format(mesh.GetNumberOfPoints()))
  print("  Number of elements: {0:d}".format(mesh.GetNumberOfCells()))

  # Write the mesh.
  #mesh_series = sv.meshing.Series()
  #mesh_series.write(os.path.join(directory, 'Meshes', 'mesh.msh'))
  mesh_directory = os.path.join(directory, 'Meshes')
  mesh_path = os.path.join(mesh_directory, 'mesh.vtu')
  print(mesh_path)
  mesher.write_mesh(mesh_path)
  model_name = 'mesh'
  write_msh(mesh_directory, model_name, options)
  write_XML(mesher.get_surface(), os.path.join(mesh_directory, 'mesh.vtp'))

  areas = []
  for face_id in outlet_ids:
    face_poly = mesher.get_face_polydata(face_id)
    areas.append(compute_face_area(face_poly))

  write_files(mesher, directory, wall_id, inlet_id, outlet_ids, outlet_names)

  print("Successully generated the mesh *****\n")
  return mesher, areas, wall_id, inlet_id, outlet_ids
  

[David Parker]

HI Jimmy,

You need to have centerlines for radius-based meshing. The centerlines geometry can be read in from a file and passed to sv.meshing.TetGenOptions

Code: Select all

options.radius_meshing_centerlines = centerlines

See https://github.com/SimVascular/SimVascu ... meshing.py

The is also an options.radius_meshing_compute_centerlines option you could try but I just tested this and it failed to correctly determine the caps for the model I was using, not sure why.

Cheers,
Dave

[Jimmy Azarnoosh]

Hi Dave,

Thanks for your help. We actually got it. options.radius_meshing_centerlines = centerlines didn't work, though. But we added two lines that worked and got the mesh similar to sv gui.

Code: Select all

options.mmg = False
options.radius_meshing_compute_centerlines = True

The only issue is that changing the value of edge size doesn't change the mesh quality. We used large and small edge size and saw minor differences in the number of elements. Is there anything we are missing?


Thanks,
Jimmy

[David Parker]

Hi Jimmy,

Maybe if radius-based meshing is enabled then edge size is not used or it is only used if mmg is enabled to create a surface mesh.

Cheers,
Dave

[Jimmy Azarnoosh]

Hi Dave,

The issue was with radius_meshing_scale. The default value for this is 0.1. To fix this issue, we used our assigned edge size for radius_meshing_scale and got the same mesh that sv GUI generates. Here's the updated options in the code:

Code: Select all

  print("\nSet meshing options ... ")
  options = sv.meshing.TetGenOptions(global_edge_size=edge_size, surface_mesh_flag=True, volume_mesh_flag=True)
  options.radius_meshing_on = True
  options.boundary_layer_inside = True
  options.radius_meshing_scale = edge_size # default is 0.1
  options.use_mmg = False
  options.radius_meshing_compute_centerlines = True

Thank you for your help

Cheers,
Jimmy

[David Parker]

Hi Jimmy,

Very useful information, thanks!

Cheers,
Dave

[Jimmy Azarnoosh]

Hi Dave,

Just a follow-up question. I used the same mesh setup for the models shown below using a global max edge size of 0.015. RCA (left one) looks okay with a reasonable surface element size. But for LCA (the right one), large elements are visible. I think that's due the the stenosis as I marked in the figure. If there is no stenosis I would expect to get a similar mesh to RCA. The solution that we found is to use a smaller global max edge size. The question is how we can determine what edge size we better use in case we have stenosis?


We tried local size on the surface wall but failed for LCA. Is it possible to assign the maximum element size to prevent elements larger than a certain size? What do we need to add to the code for that?

Thank you so much,
Jimmy

[David Parker]

Hi Jimmy,

Are you using radius-based meshing for this? I would think that that would correctly mesh the stenosis.

You could also try sphere refinement.

Cheers,
Dave

[Jimmy Azarnoosh]

Hi Dave,

Sorry for the late reply. Yes, we used radios-based meshing that correctly mesh the stenosis. We actually want to get rid of large elements near aortic. I think the only parameter to control element size is global max edge size. When we have stenosis, this value should be smaller to get a reasonable element size for large vessels. Anyway, I think we are good. Thanks for your help


Best,

Jimmy