From 3D surfaces to CFD meshes

This tutorial demonstrates how to process a surface model (obtained like in this tutorial or with other techniques) to generate a computational mesh for use in CFD.

We start off by assuming you have a 3D surface model of a vascular segment with blobby closed ends, like the ones generated using level sets or Lagrangian deformable models. If you already have a model with open inlets and outlets, just skip the Opening the surface section.

Smoothing the surface

Image segmentation can result in bumpy surfaces, especially if the image quality is not high and one didn't use any curvature term in level sets evolution. Since artifactual bumps in the surface can result in spurious flow features and affect wall shear stress distributions, one may want to increase surface smoothness prior to building the mesh.

vmtksurfacesmoothing is a wrapper around vtkWindowedSincPolyDataFilter, which implements one of the non-shrinking filters by Taubin.

By typing vmtksurfacesmoothing --help we see that there are two main parameters controlling the amount of smoothing: passband, which is the cut-off spatial frequency of the low pass filter, and iterations, which is the number of smoothing passes. For more details, visit the vtkWindowedSincPolyDataFilter doxygen page.

For typical vessels, the following should be ok

 vmtksurfacesmoothing -ifile foo.vtp -passband 0.1 -iterations 30 -ofile foo_sm.vtp

If you want more smoothing, try using a passband of 0.01. Be careful not to kill surface features by smoothing too much. Also, watch out the apex of bifurcations, since its curvature may decrease resulting in a shallower apex and affecting the simulated hemodynamics.

If you want to compare the smoothed and original versions, to make sure that no shrinking occured and the main surface features were preserved, try this

 vmtksurfacereader -ifile foo.vtp --pipe vmtksurfacesmoothing -iterations 30 -passband 0.1 
 --pipe vmtkrenderer --pipe vmtksurfaceviewer -display 0 --pipe 
 vmtksurfaceviewer -i @vmtksurfacereader.o -color 1 0 0 -display 1

You'll see the original surface in red and the smoothed surface in gray.

A potential side effect of the filter is that it slightly shifts the position of the surface in space. If this occurs, try increasing iterations

Opening the surface

Under the assumption that you generated the surface using a deformable model, it's likely that your surface is closed at inlets and outlets, with a blobby appearance. We'll now proceed by opening the surface by clipping the blobby endcaps.

You have at least three options.

The first option is to use vmtksurfaceclipper. If you call

 vmtksurfaceclipper -ifile foo.vtp -ofile foo_cl.vtp

a rendering window will show up.

Press i to start the interaction. A cube will appear (like in vmtkimagevoiselector). Position the cube in such a way that the portion of the surface you want to clip lies inside the cube (it would be easier with a cyberglove, uh?). Try to position one face of the cube perpendicular to the vessel. You don't need to be extremely precise, just do the best you can. It'll get easier with time.

Press the space bar to proceed with clipping.

Press i again if you want to clip another piece, or q if you want to quit.

The second option is to use Paraview. The procedure is not dramatically different from what you do with vmtksurfaceclipper, but the interface is surely friendlier.

The third option is to clip endcaps automatically. No 3D interaction involved. Endcap clipping can be performed using vmtkendpointextractor, which needs centerlines to be computed beforehand with vmtkcenterlines. Let's set up the corresponding pype

 vmtksurfacereader -ifile foo.vtp --pipe vmtkcenterlines --pipe 
 vmtkendpointextractor --pipe vmtkbranchclipper --pipe 
 vmtksurfaceconnectivity -cleanoutput 1 --pipe vmtksurfacewriter -ofile foo_ct.vtp

When prompted, select one source point on one inlet or outlet and press q. Then select as many target points as the rest of inlets/outlets and press q.

If you look at the output file, you'll see that the endcaps have been removed. If you want to clip a larger extent of the endcaps, just use the vmtkendpointextractor option numberofendpointspheres (the default is 2, increase it to 3 or 4 and see what happens).

Increasing the number of surface triangles (optional)

Sometimes, with low-resolution images, small vessels or stenoses, the number of triangles defining the surface can be low. This can cause potential problems to some mesh generation algorithms. The solution to this is to subdivide the mesh using a smooth subdivision scheme like Butterfly (the original points are preserved and new ones are added in such a way that at the limit it would produce a C1-continuous surface) or Loop (the original points are displaced but the limit surface is C2 continuous).

If you need to subdivide your surface, do the following

 vmtksurfacesubdivision -ifile foo.vtp -ofile foo_sb.vtp -method butterfly

Adding flow extensions

Flow extensions are cylindrical extensions added to the inlets and outlets of a model. They are important for ensuring that the flow entering and leaving the computational domain is fully developed, so that fully developed boundary conditions aren't forcing the solution in the actual vessel.

Adding flow extensions is a typical problem in CFD modeling. Stitching cylindrical flow extension to an inlet or outlet of a realistic vessel is not a trivial task, and may result in worsening the reproducibility and adding operator-dependence to the modeling procedure. A fully automatic procedure can solve this often overlooked problem preserving reproducibility and speeding up the modeling phase considerably.

Once again, adding flow extensions relies on centerlines. This time centerlines can use open inlet and outlet profiles for the definition of seed and targets.

Let's see how this is done

 vmtksurfacereader -ifile foo.vtp --pipe vmtkcenterlines -seedselector openprofiles
 --pipe vmtkflowextensions -adaptivelength 1 -extensionratio 20 -normalestimationratio 1 -interactive 0
 --pipe vmtksurfacewriter -ofile foo_ex.vtp

The adaptivelength argument of the vmtkflowextensions script is a boolean flag which enables computing the length of each flowextension proportional to the mean profile radius. The proportionality factor is set through extensionratio. The normalestimationratio argument controls how far into the centerline the algorithm looks for computing the orientation of the flow extension.

In the previous line, the flag -interactive 0 was specified. This means that vmtkflowextensions will not prompt the user about what inlet or outlet to extend, but it will perform the task on all the available open boundaries. The default behavior is -interactive 1, which prompts the user about which boundaries to extend through a graphical window. Once the rendering has started and you have determined what extensions to generate, press q and you'll be asked to list the ids of the desired boundaries. Remember that you can pipe more than one vmtkflowextensions script one after the other if you need to perform the task repeatedly.

Generating the mesh

Surface remeshing (vmtk)
Volume meshing (Tetgen)

Both are bundled in vmtkmeshgenerator, which takes care of capping the surface, assigning entity ids to wall and caps, remeshing the surface and meshing the volume.

Uniform element mesh

 vmtkmeshgenerator -ifile foo.vtp -ofile foo.vtu -edgelength 0.5

where -edgelength is the absolute nominal length of a surface triangle edge.

Radius-adaptive element mesh

 vmtksurfacereader -ifile foo.vtp --pipe vmtkcenterlines -endpoints 1 -seedselector openprofiles 
 --pipe vmtkdistancetocenterlines -useradius 1 --pipe vmtkmeshgenerator -elementsizemode 
 edgelengtharray -edgelengtharray DistanceToCenterlines -edgelengthfactor 0.3 -ofile foo.vtu

Converting mesh elements to quadratic

 vmtklineartoquadratic -ifile foo.vtu -ofile foo_q.vtu -rfile foo.vtp -entityidsarray CellEntityIds

Scaling the mesh

Medical images are often in mm. Typically you want your computations and results to be either in cgs or SI systems, which means scaling the model in cm (1e-1) or in m (1e-3). Let's scale our mesh in cm and store it in the usual vtu format.

 vmtkmeshscaling -ifile foo.vtu -scale 0.1 -ofile foo_scaled.vtu

Inspecting the mesh entities

 vmtkmeshboundaryinspector -ifile foo.vtu -entityidsarray CellEntityids

Exporting the mesh for your solver

Fluent:

 vmtkmeshwriter -ifile foo.vtu -entityidsarray CellEntityIds -ofile foo.msh

Dolfin:

 vmtkmeshwriter -ifile foo.vtu -entityidsarray CellEntityIds -ofile foo.xml

libMesh:

 vmtkmeshwriter -ifile foo.vtu -entityidsarray CellEntityIds -ofile foo.xda

LifeV:

 vmtkmeshwriter -ifile foo.vtu -entityidsarray CellEntityIds -ofile foo.lifev

The Netgen way

Capping the surface

Before meshing we may need to generate flat circular endcaps at the open inlet/outlet profiles. This depends on the meshing tool we use. Netgen requires it.

 vmtksurfacecapper -ifile foo.vtp -ofile foo_cp.vtp -interactive 0

Again, -interactive 0 means that all open profiles will be capped. If you need just a few of them capped, set -interactive 1 (or omit it, since 1 is the default), and you'll be queried what profiles you want capped.

Exporting the surface

We can now export our surface in stl format, which Netgen will import.

 vmtksurfacereader -ifile foo.vtp --pipe vmtksurfacewriter -ofile foo.stl

or equivalently

 vmtksurfacewriter -ifile foo.vtp -ofile foo.stl

(there's usually more than one way to perform tasks).

Remember that using PypeS you can pipe the previous steps and perform them in one shot without writing and reading the intermediate files each time.

Meshing with Netgen

Netgen is a great open source mesh generation tool. With Netgen you can generate high-quality meshes from stl files in a matter of minutes. In the following I'll give you some advice on how to get started with Netgen.

First of all download Netgen and follow the instructions for installing it on your system. When you're done, fire it up with ng.

A nice tkInter-style GUI will pop up. You're encouraged to read the manual, in order to be able to explore the capabilities of this mesh generation tool. For now, let's go straight to the point.

Click on File->Load Geometry. Select your foo.stl file. Netgen will read it and try to figure out feature edges (i.e. particular edges whose connected triangles form an angle greater than a specified threshold - feature edges are preserved during meshing, in such a way that they'll become edges of the final mesh). Typically, in a realistic vessel, the only feature edges are located at the boundary between the vessel wall and the inlet/outlet endcaps (although exceptions could exist). If surface triangles are reasonably dense and the surface is smooth, Netgen will automatically identify only boundary inlet/outlet edges as feature edges. However, it's always a good practice to check if this is true and eventually adjust the threshold feature angle manually.

Do inspect feature edges, click on Geometry->STL Doctor and select the Edit Edges tab. You'll see that edges appear in yellow on the 3D surface. You can interact with it to check for the presence of spurious feature edges. In order to reduce the number of feature edges by increasing the feature angle threshold, move the first slider build edges with yellow angle. As you move it, you'll see weaker feature edges disappear. Since flow extensions end with approximately 90-degrees angles, feature edges at the boundary between vessel wall and inlet/outlet endcaps will be particularly stable. Once you're satisfied, click on Build Edges and close the STL Doctor window.

We're ready to mesh now. Click on Mesh->Meshing Options, and select a Mesh granularity or fine or very fine (you can play with all the possibilities and see what happens, anyway). If you want Second order elements (as I do) activate that option. If you click on the Mesh Size tab, you'll see a few parameters that will affect the final mesh. Again, take some time to play with them. Notice that changing the mesh granularity will basically select a predefined set of parameters in the Mesh Size tab. Just a hint, STL - chart distance has a great influence on the resulting mesh, as also Elements per curvature radius does.

At this point, click on the Generate Mesh button on the toolbar of the main application window, and meshing will start. You'll see the mesh being generated real time in the rendering window. The main phases the mesher will go through are Analyze Geometry, Mesh Edges, Mesh Surface, Optimize Surface, Mesh Volume, Optimize Volume. You can change the behavior of all these phases by playing with the meshing options.

As the mesh gets generated, you'll notice a very nice feature of Netgen: element density is related to the curvature radius, which is very useful since smaller vessels will be meshed with a denser mesh.

Let's take care of defining the boundary entities now. Click on Mesh -> Edit Boundary Conditions and click on the prev and next buttons. You'll see that Netgen has assigned a different bc property to each group of elements bounded by a closed feature edge. This way we have our inlet, outlets and wall automatically defined by the feature edges and identified by a bc property (a bug in the last version of Netgen will prevent the rendering window to update surface colors according to the selected boundary face on some systems - at least on mine! No worries, it's just a display problem, Netgen will get the faces right). If you want, you can change the bc property value assigned to each face.

Well, that's what you need to know to get started. Netgen offers a lot of more advanced features, but you can explore them once you get accustomed with the basics. The last step is to export the mesh. Netgen offers several possibilities of exporting the mesh for particular solver packages. Click on File -> Export Filetype to see them all. Alternatively, you can export the mesh in the Netgen neutral file format, which is particularly easy to read with a home-made program. vmtk is able to import Netgen neutral files, as shown in the next section.

Once you're done selecting the solver package, click on File -> Export Mesh, input the file name and you're all set!

Scaling the mesh

 vmtkmeshscaling -ifile foo.neu -scale 0.1 -ofile foo.vtu

Generating an input file for a CFD solver

This step greatly depends on the solver you use. I've got a vmtk script that generates an input file to the solver I use (newtetr), which is not available under an open source licence. Contact if you're interested in this solver.

As reported above, Netgen will export the mesh for a number of solver packages out there. Otherwise, contact me for information on how to write an input file generator for a particular solver to include in vmtk. Writing mesh I/O stuff is one of the most boring activities I know, so any help would be greatly appreciated...