BioMesh  Project

An All-Hex Meshing Strategy for Bifurcation Geometries

Software Download
User's Guide

Single Bifurcation Model

Multiple Bifurcation Model

Data Source: Prof. Francis Loth, UIC

 Data Source : Arizona State University


Bifurcation is very common in all natural flow carriers (example arteries, veins etc.).  Most of the numerical flow simulation using  Navier-Stokes equations require a high quality discretization of computational domain into 2D (Triangles, Quadrilaterals or Polygons) or 3D simplexes ( tetrahedra, hexahedral, prism, or hybrid).  The selection of element shape is decided by the solver, accuracy and resource requirements. Unstructured tetrahedral mesh, although, much simpler to generate,  may require large number of elements to resolve very small scale phenomenon compared to structured or unstructured hexahedral  mesh.
In this project, we implement a robust algorithm to generate All-Hex elements in bifurcation and tubes arbitrary shapes  arising in subject specific computational hemodynamics modelling. The key components of our approach are the use of a natural co-ordinate system, derived from the solutions of Laplace's equations that follows the tubular vessels. We demonstrate with many different geometric models, that a very high quality hex mesh can be generated in a fast and accurate way.

Data Sets:

For the following geometric models, you may download
1.   Surface Triangulation
2.   Tetrahedral Volume Mesh
3.   Hexahedral Volume Mesh ( Level I, II and III )

Qmesh1 Qmesh Qmesh Qmesh Qmesh Qmesh Qmes

Acknowledge:   All the visualization have been done using ParaView Software.


Table of Contents

1.    Introduction
2.    Procedure
2.1     Surface Reconstruction
2.2     Classify Surface Vertices
2.3     Tetrahedral Volume Mesh Generation
2.4     Solving Heat Conduction Problems
2.5     Constructing Primary Chips
2.6     Constructing Split and Simple Chips
2.7     Applying Mesh Templates on each Chip
2.8     Smoothing QuadMesh on each Chip.
2.9     Create Volume Mesh
2.10   Apply Volume Mesh Optimization and Obtain Quality Measures
3.    Performance Results

A general pipeline starting from image acquisition from various modalities till the mesh generation for the numerical simulation is as follows.  One of the main goal of this project is to automate the entire process with minimal human intervention.


Given a general bifurcation geometry as shown below, our algorithm will generate a high quality all hexahedral mesh in the tubular domain satisfying two requirements:

  1. No two cross sections intersect.
  2. Each cross section is orthogonal to the vessel wall.
  3. Number of hexahedral elements must be small and predictable.


       Fig:  (I)   Domain decomposition of a vascular geometry into six regions 
(II)  Each region is mapped to a canonical  half cylinder

Some existing general purpose all-hex mesh algorithms

  1. Whisker Weaving ( Timorthy J. Tautges )
  2. Plastering ( Steven J. Owen )
  3. Medial Axis Decomposition
  4. Octree

All the above algorithms have been used in varieties of applications. Compared to a general purpose approach, our domain specific approach has the following features.

1.    Specifically designed for tubes and bifurcations,
2.    Easy to implement  with reusable software components ( > 90 % reuse)
3.    Produce very high quality mesh.
4.    Number of mesh elements and spacing can be controlled.
5.    Numerically robust. 

We have been using these geometries to study turbulence in human carotid arteries which develop atherosclerosis causing stokes. For more details, read an article from Dr. Paul Fischer entitled "Stroke Busters in Turbulent Blood

Domain Decomposition:

The fundamental problem with the bifurcation problem is related to "Domain Decomposition". How can we divide the three vascular tubes into three domain, each of the domain should be topologically equivalent to a one-dimensional tube. The following pictures describe the decomposition performed by our method.

                                                        Fig:  Domain Decomposition of vascular geometry into three tubes.

Basic Idea:

Consider an undulating tube with varying diameter. Let us solve the Steady thermal conduction problem with the following boundary conditions
  1.      T =  Ta   At the opening cap;
  2.      T  = Tb   At the end cap
  3.      On the side walls  DT/DN = 0 where N is outward normal to the boundary surface.
Assume that Tb > Ta  and now construct the  isosurfaces of  monotonically increasing sequence of temperatures Ta = a1 < a2 < a3 .... =Tb. These isosurfaces will be non-intersecting and orthogonal to the wall.  In other words, these isosurfaces create a natural coordinate system.

                                           Fig:  Isosurface in an undulating tube creates a natural coordinate system.

We want to extend this idea create all-hex volume mesh in bifurcation geometries.  In the following sections, we will explain the entire process. For greater details, refer to our publication.

Short Description of the Procedure:

Step I:   Surface Reconstruction

Let us assume that image segmentation and registration is already performed and that produces a noisy point cloud.  The first step in the generation of the volume mesh is to produce a high quality surface representation. Among the various schemes, Implicit Surface Reconstruction using Radial Basis Functions ( RBF) and Surface Reconstruction using Delaunay Triangulation (Dey, Ameta etc.) are the most popular choices.  For detail descriptions,  refer to SRPC.

The following pictures shows surface reconstruction using RBF and Delaunay Triangulation (DT). Perhaps the biggest  advantage of RBF is in filling gaps in an undersampled dataset.  As shown in the following picture II, RBF can also fill a gap in some unwanted regions. This necessitates visualization, manual clipping and therefore, hinders full automation. 

 Input Point Cloud
     RBF Surface
   Delaunay Surface
                                                                           Fig: Surface Reconstruction using Point Cloud

The following table shows that system of linear equations arising from thin plate RBF are ill conditioned which become worse as the number of sample points increases. 

     Matrix Size Condition Number
2.5645 x 10^5
1.1683 x 10^6
1.5976 x 10^7
Surface reconstruction using DT is extremely fast and provably correct.  Perhaps the only disadvantage of this approach comes while handling undersampled dataset. However, it is possible to detect undersampling in the dataset while creating DT, therefore, surface repair tools can be applied on the undersampled regions.

At a final step in surface reconstruction, mesh optimization schemes  1)  Edge Flipping 2)  Non-shrinking surface smoother (Taubin Smoothing, Mean Curvature flow, or Implicit Fairing )  are applied to remove additional noise in the dataset. For details, refer to  SurfSmooth                                     

Step II:   Classify Surface Vertices:

The output of step one is a clean, noise free, smooth 2D manifold surface triangulation. In order to solve the Laplace equations ( described in the next section ),  we have to classify surface vertices into three categories i.e. surface (1) , inflow  (2) , and outflow (3 & 4)  ( Assuming that in the following picture, the inflow is  at the bottom and two  branches are outflow).  At least one vertex must be classified as inflow and two vertices as outflow.

In order to automate identification of inflow and outflow vertices from the surface model, a Reeb graph (or Medial Axis) is constructed which has the shape of "Y". ( Reeb graphs is very simple to construct compared to Medial Axis). The end point of the hand is "inflow" and the fork ends are "outflow" vertices.

Surface Model

     Reeb Graph
                                                    Fig:  (I) Surface Triangulation (II)  Identifying end points  with approximate MD.

Step 3: Tetrahedral  Volume Mesh Generation

Using the surface triangulation as constraint, generate a tetrahedral volume mesh. (Many Delaunay triangulation based software may insert new vertices on the surface). Since the quality of the volume mesh is not critical for the final hex mesh generation, therefore, any suitable tetrahedral mesh generation ( of course without slivers ) can be used.  In our case, we used SolidMesh software based on Advancing  Front Technique (AFT) to obtain the tetrahedral volume mesh.

If the surface is not properly sampled or has holes, we use RBF Implicit surface reconstruction algorithm to create a smooth 2D manifold surface. In this procedure, too many triangles are generated that are needed to solve the heat conduction problems and subsequently mesh  generation. Therefore, a coarser grid is obtained by  finding "Independent Set" from the mesh graph.  "Topological Noise" is the biggest problem in the surface reconstruction through implicit polygonization. The holes can be identified by calculating Euler Characteristic of the surface. At present, we don't have robust implementation to remove the "Topological Noise", therefore, we experimented with different spacing in MC and were able to get zero genus surfaces.

Step 4:  Solving Heat Conduction Problems:

As described in the paper the central theme of our approach is to create a natural co-ordinate systems by solving three different Laplace Equations with the following Dirichlet and Neumann Boundary Conditions.  For this purpose, standard Finite Element Formulation ( FEM) is used.  The symmetric, positive definite system of linear equations arising from FEM is solved using Preconditioned Conjugate Gradient method. ( We used Matrix Template Library MTL for Matrix objects and Iterative Template Library ITL for PCG) . The color plots of the temperature distribution are shown in the following pictures.

Branch A:     Insulated
Branch B:     -Ta
Branch C:     1-Ta
Branch A:       1-Tb
Branch B:        Insulated
Branch C:       -Tb
Branch A:        -Tc
Branch B:         1-Tc
Branch C:        Insulated
                                                                                  Fig:     Temperature distribution

Observe that there is rotational symmetry in the boundary conditions, which can be exploited to reduce the code size. How should we select the Dirichlet Boundary condition values ?  We can choose any  two  values, but any arbitary values may not give a same value on the insulated cap for all the three cases ( It is must and the reason will be explained later). Therefore, for each case, we solve  Laplace solver twice in order to ensure that the end cap of each insulated boundary is maintained at the same value. For each case, in the first run,  we specify Dirichlet values 0 and 1 at two ends. This will create some temperature ( let us say alpha) at the cap of insulated boundary. In the second run,  Dirichlet Values are shifted by alpha, so that the temperature at the center of the insulated end cap becomes zero.

The Isosurface emanating from the center of insulated caps are called "Principal Isosurfaces".                         

Step 5:  Constructing Primary Chips:

Obtain the three iso-surfaces at T = 0 from each solution. Ideally, with the exact arithmetic, we should get a smooth surface from the centre of insulated cap and right upto the end near the bifurcation, but as shown in the pictures, iso-surface away from the centre are jig jagged and sometimes there may be pockets of isolation.  This behaviour is because of finite numerical precision, but our region of importance (ROI) is the "Bifurcation Centre", where all three iso-surface are smooth.


                                                        Fig:   Isosurfaces from the solution of heat conduction problems.

IsoSurface A
IsoSurface B
IsoSurface C
Intersection of IsoSurfaces
                              Fig: Creating principal chips with principal isosurfaces which emanates from the insulated BC branch.
Use implicit surface cutter to remove the part of the surface that is not required.        


Zoomed Picture


Primary Chips

                       Topology of A Chip
                           Fig:  (I) A zoomed picture of Principal Chips ( A and B).  (II)  Topology of each chip. (III)  A half cylinder

"South Vertex"  of the primary chip acts as a reference point for the proper orientation of the subsequent chips. Three curves (left, right, and Interface boundary) are parameterized with south (u = 0.0) and north (u= 1.0). East, West and Center vertices are then determined at u = 0.5 along these curves.

Step 6: Constructing Split and Simple Chips:

As we move away from the bifurcation region, the angle between the left and right chip decreases ( as shown in the zoomed picture below) . When the angle between right and left chip become less than 10 degrees, ordinary chips are  constructed from the isosurface in the branch.


Chips in Artery

Zoomed Region Near the Bifurcation
Note:  Simple chips near the end caps may be distorted or cut depending upon the geometry and boundary conditions applied on selected vertices. There are two options  (1)  Reject the chips that are not proper  (2) Start with an artificially extended branch and then discard all the chips in the extended branch.

Pseudo-codes for constructing Split and Simple chips are given below. Note that Isosurface extraction and implicit surface cutter are the only two software modules that are required to create chips. Also note that, isosurface extraction from tetrahedra mesh will be called large number of times and it is the most expensive compared to other calls, therefore, any optimization can have significant improvement in the efficiency of the algorithm.

Algorithm: ConstructSplitChip ( tetmesh, branchID,  isoVal, center0)

             // Input Parameters:
             // branchID   :  The branch in which which Split chips are created.
             // IsoVal       :   Temperature value for which the right surface of the chip is constructed.
             //                       The left side depends on the location of the right ( because of alignment).
             // center0      :   Approximate center of the previously generated chip.

             bA = ( branchID+0)%3
             bB  = ( branchID+1)%3
             bC  = ( branchID+2)%3
              //  Search for right half of the chip separated by distance "width"
              while(1) {
                     rightmesh  = getSurface( tetmesh, isoval, bA)
                     cutmesh    = getCutSurface( rightmesh, 0.0, bC)  // Implicit Cut the surface
                     interface    = getInterface(cutmesh)                      // Extract the Cut boundary contour
                     center1      = getCenter( cutmesh)                        // Get the center of interface contour
                     if( dist( center0, center1) > width) break;
               discard( rightmesh, -1)                         // Remove the left portion of the righmesh.
               addInterface( chip,  interface)               // Add the interface
               addSurface( chip, rightmesh, 1)            // Remaining mesh constitutes right half of the chip.

               // Create left half of the chip, such that interface align together.
               isovalB  = getIsoValue( center1, bB)         
               leftmesh = getSurface( tetmesh, isoValB, bB)
               cutmesh = getCutSurface( leftmesh, 0.0, bC)
               discard( leftmesh, 1)                            // Remove the right side of the leftmesh.
               addSurface( chip, leftmesh, -1)          // Remaining mesh constitutes left half of the chip.

Algorithm : ConstructSimpleChip( tetmesh, branchID, isoVal, center0)
                          trimesh = getSurface( tetmesh, isoVal, branchID)
                          addSurface( chip, trimesh)

Step 7:  Apply Mesh Templates  on each Chip:

Each Chip can be divided into four quadrants and user defined mesh templates as shown below can be applied on each chip. Vertices are always projected on the chip surface constructed by the isosurface. With increasing level of subdivisions, better approximation of the outer surface is obtained.


             Level-0 Mesh Template
               Level-I Mesh Template

             Level-II Mesh Template
                                                                                 Fig:   Examples of Mesh Templates.

Step 8:  Smoothing QuadMesh on each Chip:

Isosurface generate  triangulated surface which acts as a background mesh on which the quadrilateral mesh is superimposed. The quality of quadrilateral mesh depends on the node distribution on the background mesh. As the number of levels increases, two or more points may project on the same node of the background mesh which will result in zero edge length of some quadrilaterals.  Therefore, it is essential to apply smoothing, after mesh template is applied. Experiments shows that, 4-5 steps of Taubin's smoothingis sufficient to get very good quadrilateral mesh.

Step 9:  Create Volume Mesh:

After consistently orienting each chip with respect to the primary chip, volume mesh is constructed using quad-mesh from  two consecutive chips. As the number of chips and quad-mesh is deterministic, total number of hex-elements is also known. In order to ensure that each hex is positive, Jacobian at each eight nodes is checked.

Step 10:  Apply Volume Mesh Optimization and Obtain Quality Statistics.

The quality of the volume mesh may be further improved by using standalone Mesh Optimization Tools ( such as Mesquite).  Volume mesh, however, can not be improved beyond a limit unless the surface mesh is also improved and/or topological flipping operations are carried out. (At present, we have not implemented these optimization)

Out of many quality measures, we measure only few to ensure that the correctness of the final volume mesh.
1.    MinMax Volume.         
2.    MinMax Aspect Ratio
3.    Worst Condition Number.

Performance Results:

Time ( in seconds)
Vertex Classification

Heat Conduction Problems

 Chips Construction
 1. Primary Chips
 2. Split Chips
 3. Simple Chips

Mesh Templates
1.  Level  I
2.  Level  II
3.  Level  III

Mesh Smoothing

Volume Meshing

Volume Mesh Optimization

Quality Checks


Presently, the source code is not ready for public distribution.  Please contact  Dr. Paul Fischer (E-mail: for  further information.


File Format:

All the files are in ASCII and the data format is as follows. For order of connection, refer to MeshDatabase

#Nodes Number of Nodes
x1  y1  z1
x2  y2  z1
xn  yn  zn

#Triangles Number of Triangle
v11   v12  v13 
v21   v22  v23 
vn1  vn2   vn3
#Quads Number of Quads
v11 v12 v13 v14
v21 v22 v23 v24
vn1 vn2 vn3 vn4
#Hex Number of Hex
v11 ...           v18
v21 ...           v28
vn1 ...            vn8

Note:  All vertex index starts from 0 ( not 1 as in FORTRAN) and are contiguous.


 14th IMR Conference, San Diego Presentation 13th Sep, 2005  Download Presentation.
 (The paper is also available in the Proceedings of the 14th International Meshing Roundtable )


        Chaman Singh Verma, Paul F. Fischer, Seung E. Lee, F. Loth
        14th International Meshing Roundtable Conference-11-14 Sep, 2005, San-Deigo, California
        Download Paper

Contact :

Paul Fischer
Mathematics and Computer Science Division
Argonne National Lab
Francis Loth
UIC College of Engineering,
Chaman Singh Verma
Mathematics and Computer Science Division
Argonne National Lab.

Last Modified: Wed Sep 31 17:03:06 CDT 2006