Introduction

Mesh Generation

Discretization, also known as grid or mesh generation, is defined as the process of breaking up a physical domain into smaller sub-domains (usually called elements). Discretization is necessary in order to facilitate the numerical solution of partial differential equations. In other words, there are few analytic solutions for geometries of engineering interest thus the geometry must be divided into an aggregate of simpler pieces for analysis.

Automatic mesh generation has been an area of intense research for decades, and a tremendous amount of literature and numerous algorithms have been developed. There are three fundamental challenges in the field: robustness, mesh quality, and computational efficiency in generating the mesh. A recent survey indicated there were over 80 commercial and academic meshing products available, of which 39 automatically generated tetrahedral (“tet”) elements compared to 20 that performed unstructured hexahedral (“brick”) mesh generation. The current dominance of tetrahedral meshing can be attributed most notably to its ability to robustly mesh arbitrary, complex geometries. In addition, the use of tetrahedral elements often simplifies the process of adapting the mesh during simulation.

Once we have a geometric (solid) model representing a device of interest, we need to discretize the flow domain into little pieces (called elements) for simulation. This process is known as mesh generation.

We will review several important concepts prior to discussing mesh generation:

  1. For the applications discussed herein, a geometric model representing each energy domain of interest (fluid, mechanical, thermal, etc.) must exist. That is, if you are simulating the mechanical forces of a balloon expandable stent deployed in air, you need a geometric model representing at a minimum the stent geometry. However, if you want to simulate flow around a stent deployed in a vessel, you need a geometric model of the flow lumen (i.e. flow domain) not of the stent itself.

  2. There are multiple sources of geometric models. The most common source of geometry in the device community is likely from commercial CAD packages (SDRC I-deas, SolidEdge, SolidWorks, Pro-E, etc.). A second potential source of geometric models is digitized data (imaging data, etc.). Finally, a less common (but increasingly important) source of data is the results of a previous simulation. For example, if one simulated the deployment of a stent into a diseased artery, the results of this analysis (i.e. displacements) could be used for the initial geometry for a subsequent blood flow analysis around the deployed stent.

  3. As will be shown using a beam example, a simple geometry allows for great flexibility in choice of element type (e.g. tetrahedral vs. hexahedral elements). However, for many geometries of engineering interest, the complex geometry motivates the need for automatic mesh generation and often restricts the user to tetrahedral or possibly hex-dominate meshing. Pure hexahedral mesh generation for complex geometries is an area of active research.

In this manual, a finite-octree based tetrahedral mesh generator is used. The basic idea behind finite-octree methods is to decompose a complex geometry into simpler pieces and then mesh the individual pieces using a mesh generation technique (e.g. templates, Delaunay triangulation (see Figure 1), etc.). An example quadtree decomposition (the 2-D analog of finite-octree) is shown in Figure 2. A structured (i.e. tensor-product) quadtree grid that completely fills the space occupied by the geometric model (i.e. bounding box) is initially created based on a user-defined mesh density. Subdivision of the octants is then performed to reach a desired complexity of the geometry contained therein. Restricting the “level difference” or gradation between adjacent octants is essential to preserve mesh quality. After the geometry is decomposed, surface meshing is performed using projected 2-D Delaunay surface triangulation. Templates are used to create the interior volume mesh (i.e. octants contained completely inside of the geometric model), and 3-D Delaunay triangulation is used in the boundary octants (i.e. octants containing part of the geometric model boundary) to finish the meshing process.

Finally, we note that there are several techniques to evaluate the quality of a given discretization. For the exterior surface mesh, a visual inspection may provide useful information. However, it is impractical to visualize the individual tetrahedral elements for large meshes, so geometric-based mesh quality indicators are used. Specifically, three mesh quality indicators will be discussed (see Figure 3):

  • Minimum solid angle
  • Radius ratio
  • Aspect ratio

In the present work, a combination of these quality indicators is used by an iterative mesh optimization algorithm to improve the overall quality of the mesh. It should be noted that there is current research interest in generating error estimators that include both solution and geometric information that may lead to more accurate solutions while requiring ewer elements.

Figure 1. Delaunay criterion. The Delaunay criterion states that no other point in the triangulation can fall within the circumscribing sphere (circle in 2-D) of the points defining a simplex in the triangulation. Figure (a) shows a valid Delaunay triangulation of four points in $\mathbb{R}^2$ while (b) shows a non-Delaunay triangulation of the same four points. In 2-D, the Delaunay criterion minimizes the maximum interior angle producing an optimal triangulation for a given set of points.
Figure 2. Quadtree decomposition. The figure shows an example quadtree decomposition (directly analogous to octree decomposition in 3-D) that is used to divide the geometry into less complex individual pieces for automatic mesh generation.
Figure 3. Geometric mesh quality measures. Shown are 2-D geometric mesh quality measures with direct analogies in 3-D. The radius ratio (a) is the ratio of the radius of the maximum inscribed circle (sphere in 3-D) over the radius of the circumscribed circle (sphere in 3-D). The aspect ratio (b) is a ratio of the minimum height to the maximum base length. The maximum/minimum interior (dihedral in 3-D) angle is shown in ©. In this lab we will construct meshes automatically for two different geometries. The first example is that of an idealized vessel (i.e. cylinder). The second example will be an idealized stent, deployed in an idealized stenotic vessel with incomplete apposition.

SimVascular Meshers

This document describes how to use the Meshing software for discrete and continuous solids. SimVascular meshing includes both open source and commercial options. The open source meshing includes the open source libraries of TetGen and the Vascular Modeling Tool Kit(VMTK). These tools are combined to provide boundary layer, mesh refinement, and isotropic adaptive meshing. The commercial mesher in SimVascular is MeshSim and provides boundary layer, mesh refinement, and isotropic/anisotropic adaptive meshing. MeshSim is a very powerful tool and can provide high quality meshes for irregular and complicated domains.

TetGen Meshing

TetGen is an open source mesh generation software developed by Hang Si through WIAS in Berlin. TetGen is a tetrahedral mesh generator that uses 3D Delaunay Triangulation. Learn more about TetGen here.

Creating a TetGen Mesh from a Solid Model

This exercise will assume that you have created segmentations down the aorta and common iliac artery that you can now use to create a solid model.

To create a TetGen mesh (empty):

Right click the data node "Meshes" in the SV project in Data Manager
Click "Create Mesh" in the popup menu
Select Model: demo
Mesh Type: TetGen
Mesh Name: (use model name by default)
Click "OK"

Now a new data node “demo” for the mesh is created under the data node “Meshes” in Data Manager. Double click the data node “demo" and the tool “SV Meshing” automatically shows up. The new mesh is empty and has no data so far.

Isotropic Meshing

To generate an isotropic mesh, we either need to specify a desired number of elements or some property (e.g. maximum edge length for each element) of the mesh. In our case, we will specify the maximum allowable edge length for a given element for the entire mesh. We have the ability to specify our own mesh size or have SimVascular calculate a size for us.

Global Max Edge Size: click the button "Estimate" to let SimVascular to provide a value
Click the button "Run Mesher"
Click "Yes" to continue 

It may take a while fo meshing, which depends on the size of model and the estimated global max edge size. When your mesh has completed, the mesh shows in the 3D-view window of Display and a dialog pops up to show some mesh info (like the number of points, elements, faces, etc.). To see the mesh easier, hide the model temporarily.

HELPFUL HINT: The estimated edge size is about half of the radius of the smallest cap in the model.

If you zoom in, you will be able to see the individual mesh elements.

Now try to change the “Global Max Edge Size” value to generate a mesh with more elements. Note that as you specify lower global max edge size values, SimVasxcular will take longer to generate the mesh. The software might appear as not responding, but should still be working.

HELPFUL HINT: There is an option for faster meshing, to use it:

Toggle on "Fast Meshing"

It speeds up meshing with the same wall mesh, while cap mesh may not be desirable sometimes. Fast Meshing is automatically disabled if you turn on “Boundary Layer Meshing” or “Radius-Based Meshing”.

Boundary Layer Meshing

When simulating blood flow, interesting phenomenon can occur near the vessel walls. Under laminar flow, for example, boundary layers can form with large velocity and pressure gradients near the wall. It is advantageous to have increased mesh density in the areas of high gradients. If a preferential flow direction is known, you can often “elongate” elements in the direction of the flow without loss of accuracy to reduce computational costs. In the area of mesh generation applied to fluid flow, this is often referred to as boundary layer meshing. Let’s explain a few parameters used in boundary layer meshing.

* Portion of Edge Size: the portion of the edge size to use as the size for the initial boundary layer (typically this is a value between zero and one.
* Num Layers: the number of layers desired. Too many layers can cause self intersections on smaller vessels, so be careful.
* Layer Decreasing Ratio: the amount of the size decrease between successive boundary layers. This gradation factor is multiplied by the previous layer to get the thickness of the new layer. In order to decrease the layers by a factor of 2, apply a decreasing ratio of 0.5. 

We will now generate a boundary layer mesh from a solid model.

Make sure the model has correct types (wall,cap) assigned for all its faces
Go to Advance Options
Toggle on the checkbox "Boundary Layer Meshing"
Portion of Edge Size: in this case, we use 0.5.
Make sure Volume Meshing and Surface Meshing are on
Number of Layers: in this case, we use 4 layers.
Layer Decreasing Ratio: in this case, we would like each layer to be 0.6 of the previous layer. 
Make sure Volume Meshing and Surface Meshing are on
Click the button "Run Mesher"

WARNNING: Fast Meshing is automatically disabled if using boundary layer meshing.

The mesh generated will have the boundary layer mesh on the specified region. Zooming in, you can see that the number of layers in the volumetric mesh is four. The boundary layer extends all the way down the length of the surface named wall, and each layer is 0.8 of the previous layer.

It is important to note, once again, that you have generated a volumetric mesh. That is, the entire volume of the geometry has been filled with tetrahedral elements. Thus, the boundary layer meshing pattern that you see on the surface continues up the entire volume of the aorta. It is difficult to visualize thousands of elements at one time, so we only visualize this refinement on the exterior surface mesh.

Radius Based Meshing

When models contain surfaces of very different size scales, it is often desirable to scale the size of the mesh with the size of the vessel of the model. Radius-based mehsing will compute the centerlines and find the distance to these centerlines prior to meshing. These values will be normalized by the smallest value and then multiplied by the global max edge size given.

Toggle on the checkbox "Radius-Based Meshing".
Click the button "Run Mesher".

WARNNING: Fast Meshing is automatically disabled if using radius based meshing.

Note that a majority of our model has vessels of very similar size. In this case, a radius-based mesh does not do much; however, there are other models where this is a very desirable feature. The operation takes a bit long, so you can also try to change local meshing size introduced in the following section to reach the same goal. If you don’t use this operation, remember to turn it off before runing mesher.

Local Mesh Size Application

Often, a designer may have insight into locations of interest or regions of complex behavior and may want to locally increase the mesh density in those regions to improve the quality of the given finite element analysis results for a given number of nodes. In this example, we will assume the ends of a model are areas of interest and want to locally increase the mesh density near these faces. In SimVascular you can specify different mesh densities on a local geometric model face.

In this example, we are interested more in the flow in the right iliac of this model. we specify a mesh edge size of 0.1 on the right_iliac wall. We leave the global edge size unchanged.

Toggle off "Radius-Based Meshing"
Go to "Local Size"
All the faces are listed in the table.
Set loca size 0.1 for wall\_right_iliac
Or you can use table menu to set if there are many faces to set
Click the button "Run Mesher"

Spherical Refinement

When simulating blood flow, interesting phenomenon can occur around the site of vessel bifurcations. It is advantageous to have increased mesh density at and around bifurcations. In this section we will refine the mesh at the abdominal aorta and common iliac artery bifurcation using sphere refinement:

We will now refine a region with a sphere from a model of an Aorta.

Clear local size for wall\_right_iliac in "Local Size"  
Go to "Regional Refinement"
Toggle on "Sphere" and an sphere appears in the 3D-view window.
Move the sphere and change its size to make sure the bifurcation is include inside the sphere.
Click the button "Add" to add this region to the table
Provide a local size for it. Here we use 0.1
Click the button "Run Mesher"   

HELPFUL HINT: If you have multiple regions of interest for mesh refinement, you can move the sphere to another region, then add it to the table. If you want to adjust a sphere which is already added in the table, just select the sphere in the table. A red sphere will appear in its corresponding region in the 3D-view window. You can also set size or delete spheres from the table menu.

You can combine boundary layer, spherical refinment and local size to create a full mesh again.

When you are satisfied with the mesh, click the tool button “Save SV Project ” to save the mesh.

Advanced Flags

TetGen provide more advanced flags to operate meshing. These flags apply ONLY to the volumetric meshing operation.

 - O: Specity the number of times to optomize the mesh. This moves nodes to reach a better quality mesh.
 - q: Specify a quality measure for the mesh. This number corresponds to the ratio between the radius of the circumsphere of an element and the maximum edge size (See below for figure). This number can be anywhere from 2.0 to 1.1, and a lower number corresponds to a higher quality element. A mesh with a quality ratio 1.0 is not attainable and the mesher will run infinitely.  
 - T: Apply a tolerance to check whether a point lies on the surface or not (Default is $10^{-8}$).
 - Y: Preserve the exact surface mesh. If this parameter is used without a surface remeshing, make sure a mesh size that corresponds to the surface mesh is applied.
 - M: Do not merge facets that are coplanar or have very close vertices.
 - d: Detect for intersecting facets. This can be helpful if there are regions with close facets.  
 - C: Check the consistency of the final mesh.
 - Q: Output nothing to the terminal except for errors.
 - V: Print out more detailed information from TetGen. This information is viewable in the terminal.
 - Specify other mesh tags: The full TetGen documentation indicates other possible mesh flags. This is where those can be specified.

Adaptive Meshing

Adaptive meshing is used to create an adaptive mesh based on a simulation result.

Use simulation result from one single time step:

Go to Adapt
Options: Use One Step
Result File(vtu): all_results.vtu, from a previous simulation. It is created as you convert previous simulation results when "As Single File" is on. 
Step Number: Give a specifi time step number of the previous simulation
Error Reduction Factor: (0~1)Value multiplied by the average interpolation error in order to get a target uniform local error distribution. This should be a value between zero and one. A smaller factor will attempt to achieve a mesh with smaller error.
Global Min Edge Size: Specify a minimum target edge size. No edge size will be smaller than this size, even if the adaptor identifies that solution needs a edge length smaller than this.
Global Max Edge Size: Specify a maximum target edge size. No edge size will be larger than this size, even if the adaptor identifies that the solution does not require an edge length this small.
Click "Create Adapted Mesh"
Provide a name for the adapted mesh

Use simulation results from multiple time steps:

Go to Adapt
Options: Use Multiple Steps
Result File(vtu): all_results.vtu
Start Step Number: 
End Step Number:
Step Increment:
Click "Create Adapted Mesh"
Provide a name for the adapted mesh

In the end, a new adapted mesh is created under the “Meshes” in Data Manager. Also the new adapted solution is create in [proj-path]/Meshes with a name “adapted-restart.[step number].1”. You can use the adapted mesh and the adapted solution as IC file to create a new job to run simulation.

Unloading/Loading Mesh Data

Sometimes a created mesh (especially volume mesh) occupies very large amount of computer memory. In this case, after you save the mesh, you can unload the mesh data from a mesh to save memory while the mesh node is still in Data Manager.

With the same concern, when you open an existing SV project, mesh nodes are listed in Data Manager while the mesh data are not loaded for them.

To unload or load mesh data:

Right click a mesh node in Data Manager
Click "Load/Unload Surface Mesh" : load surface mesh if it has no mesh data; or unload both surface/volume mesh if it already has mesh data
Click "Load/Unload Volume Mesh" : load both volume/surface mesh if it has no mesh data; or unload volume mesh if it already has volume mesh

Exporting Mesh

To export a mesh to vtp/vtu files: a vtu file for the volume mesh, multiple vtp files for surface mesh, the whole wall surface, and individual faces:

Right click the mesh node "demo" in the SV project in Data Manager
Click "Export Mesh-Complete" in the popup menu
Select a directory and Click "Choose"

A new directory “demo-mesh-complete” is created and all the mesh files are saved to this folder.

TetGen Errors

It may be the case where TetGen cannot mesh the given PolyData and it will crash. This could be for a number of reasons including an open surface or intersecting facets. In most cases, TetGen is able to detect the error, and it outputs information about the event to the terminal. If you have built SimVascular from source, look for the error code in the output to the terminal. If you have not built from source, the best option is to try a smaller mesh edge size. It is likeley that the edge size specified is too large for certain parts of the model.

Code Message Description
1 Error: Out of Memory Not enough memory available to complete the mesh. Try increasing the mesh size
2 Please report this bug to Hang.Si@wias-berlin.de. Include the message above, your input data set, and the exact command line you used to run this program, thank you This error was caused by an unknown bug to TetGen
3 A self-intersection was detected. Program stopped. Hint: use -d option to detect all self-intersections. This failure was caused by an input error.
4 A very small input feature size was detected. Program stopped. Hint: use -T option to set a smaller tolerance. This failure was caused by a possible input error. For example, there are two segments nearly intersecting each. other. If you want to ignore this possible error, set a smaller tolerance by the -T switch, default is $10^−8$.

If the error code given is error code 2, please contact Hang Si at Hang.Si@wias-berlin.de. Otherwise, contact the SimVascular supporters or report a bug on the SimVascular website.

Inspecting Meshes

After creating a mesh in SimVascular, you can view the surface and make a visual inspection of the surface of the mesh. However, this is based purely on your idea of a “quality” mesh, and you can only see the surface of a volume mesh. In order to inspect a mesh more closely, we use open source software ParaView. ParaView contains filters to identify multiple quality measure such as volume and aspect ration.

To visualize and inspect a mesh in paraview, follow the steps below:

1) Download, Install, and Run ParaView. Go to File->Open->Choose the vtu file generated when saving the mesh.

2) Using the dropdown menus, select Filters->Mesh Quality.

3) Adjust the paramters of the Filter in the Properties tab on the left hand side of the ParaView screen.

4) Change Tet Quality Measure to the new desired mesh. Here we select Volume:

5) To visualize the metric, select the Representation to be “Volume”, and select the Array under Coloring to be “Quality”.

6) We can similarly do the same thing for many other metrics such as Minimum Dihedral Angle.

7) Apply any other qulity type to visualize the range and values of the corresponding metric.