# Presentation: 805462

## An All-Hex Meshing Strategy for Bifurcation Geometries in Vascular Flow Simulation

- Chaman Singh Verma*
- Mathematics and Computer Science Division
- Argonne National Laboratory
- Paul F. Fischer	Argonne National Laboratory
- Frank Loth	University of Illinois, Chicago
- Seung E. Lee	Massachusetts Institute of Technology
****Supported by NIH RO1*

## Outline:

- Introduction
- Basic Idea
- Methodology
- Performance Issues
- Results

## Bifurcation Geometry:

- Objective: Create all hex mesh with:
- Specified spacing (along the length and from the wall)
- Relatively few elements
- Respectable quality
- The algorithm should be reasonably fast.

##  Why all hex ?

- Targeting spectral element solution of vascular flows  – efficient at transitional Reynolds numbers****      ******	******	***
- Fast high-order operator evaluation ( O(N4) v/s O(N6) ) relies   on local tensor-product forms  *hexahedral elements*
- Spectral element meshes require orders of magnitude fewer elements than FEM meshes
*N**3** points per element, where N = [8-16]*
- T*ypical SEM meshes for bifurcation geometries require 1000 – 4000 elements, where FEM may require > 100000 elements.*

- *P.F. Fischer, G.W. Kruse, and F. Loth , *J. of Sci. Comput. ***17** 1, pp. 81-98 (2002)* *

## Limitations of Existing Approaches:

**Plastering **( Steven Owen etc. )
-      Based on *Advancing Front Technique (AFT). *Quality may suffer in the regions where front collide.
**Whisker Weaving ( Timorthy J. Tautges)**
**	**Also Based on AFT. Requires certain constraints on the surface mesh.
**Medial Axis Domain Decomposition:**
**	**Implementation of** **MAD in 3D is difficult and is very sensitive to the boundary shape.
**Octree **
**	**May have poor quality near the boundary.
*These methods are designed for fine granularity.*

## Our Goal:

- Not trying to solve the *open* *all–hex meshing problem.*
- Instead, we  are developing *domain-specific (vascular) mesh generation software*.
- Our approach is based on sweeping:

## Challenges for Sweeping Vascular Geometries

- Curvy and variable cross-sections.
- Bifurcations
- Surface descriptions come from medical images.
- 1000s of geometric models of identical topology.
*Need fast, accurate and reliable software.*

## Our Approach:  Conduction Based Sweeping (CBS)

- Develop a natural coordinate system that follows the
- undulating vessel and satisfies
*No two cross-sections intersect, and*
*Each cross section is orthogonal to the vessel wall.*
*Isosurfaces from the solution to a thermal conduction problem*
*(Laplace’s equation) satisfy these critera.*

*T=1*

*T=0*

* *

## Advantages: CBS Approach

- Very easy to implement with reusable software components( > 90% reuse )
- No intersection checks are required.
- Produce very high quality mesh.
- Mesh spacing and number of elements can be controlled.
- Numerically robust.

## Weakness of CBS Method:

- Designed specifically for tubes and bifurcation geometries.
- Suitable for those applications where end cross sections are allowed to modify or cut.

## 
Bifurcation 
Mesh Generation Procedure

**Bifurcation **
**Mesh Generation Procedure**

## Step 1: Create Unstructured Volume Mesh

- Use *Advancing Front Technique* to generate tetrahedral volume mesh for conduction problems.
- Quality of the volume mesh is not very important,   therefore, any reasonable mesh generator can be used.
- Assign boundary markers on the surface vertices
-    ( inflow, outflow, and surface ).

## Step 2: Solve Three Heat Conduction Problems

## Bifurcation: Domain Decomposition

**Bifurcation geometry is decomposed into six half-cylinder regions**
**Each half-cylinder is meshed by CBS**

## Step 3: Create Principal Chips:

-       Solution A            Solution B           Solution C                 Intersection

***Principal isosurfaces, emanating from insulated branches, define the domain decomposition***** **

*All isosurfaces from insulated branches converge exponentially to the principal isosurface.*

## Step 3: Create Principal Chips:

- Principal Chip divides the domain into six regions, each one can be topologically mapped to a half cylinder.

## Principal Chips:

- Uniqueness of common intersection can be proved in 2D case.
- In 3D, although not rigorously proven, experiments show that up-to machine precision, intersection take place along a common line.
- Exact intersection is not required for the mesh generation.

## Step 4: Create Branch Chips:

- Before Mesh Generation, consistently order all the chips
- With respect to the Primary Chips.

## Step 5: Split Chip Creation:

- Split Chips are created till the angle between left and right chip is less than 10 degrees.
- After split chip, simple chips are created till the end of the branch.

## Topology of Chip:

- Principal Chip’s Interface defines reference vertices  *South *and *North*
- All other chips are oriented with respect to PC.
- Three Curves (*left, right and Interface*) are parameterized [0,1]
*East, West* and *Center* are at u = 0.5 on their respective curves.

## Pseudo-Algorithm for Chip construction:

- 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 5: Apply Mesh Templates:

## Step 6: Projecting Quad-Mesh  On Chip

- New vertices are mapped onto the chip’s triangulated surfaces using nearest Euclidian neighbor.
- Triangulated surface and the boundary curves must be properly sampled.
***Taubin Smoothing***  is applied on both the triangulated and Quad mesh in order to remove noise.

## Step 7: Volume Mesh:

- Take two adjacent chips and create Hex element by joining two Quad elements.
- Ensure the positivity of each hex-element by calculating the *Jacobian* at each eight nodes.

## Step 8: Mesh Optimization and Quality Checks:

- The quality of final mesh further improved by mesh optimization software ***Mesquite ***
-             *Lori Freitag, LLNL, Pat Knupp, Sandia,...,Todd Munson, ANL*
- Quality Checks:
- Min-Max Aspect Ratio
- Min-Max Volume
- Worst Condition Number

## Results:

- Surface Triangulation, Tet Mesh and Hex Mesh for three different levels are available
- on the **BioMesh Project** site:     www.mcs.anl.gov/~csverma/BioMesh/biomesh.html

## 
			      Thanks

- 			      *Thanks*

## Surface Reconstruction:

|  |  |  |
| --- | --- | --- |

| Point Cloud | RBF Surface | Delaunay Surface |
| --- | --- | --- |

-  Non-Compact, thin plate RBF was used to construct implicit surface
-  Delaunay surface reconstruction algorithm (Dey, Ameta etc.)

## Performance Bottlenecks:

- From image acquisition to the final volume
- mesh generation, principal bottlenecks are:
- ( in decreasing order )
*Image Segmentation *
*Surface Reconstruction.*
*Solving Heat Conduction Problems.*
*Isosurface Extraction.*

## Image Segmentation:

- Images from Ultrasound may be fuzzy compared to CT-SCAN or MRI. In the case of free hand US, the arbitrary orientation may produce additional noise in the dataset.
- Full automation is difficult and expert intervention often provides the best solution.

## Further Information:

- Surface Triangulation, Tet Mesh and Hex Mesh for three different levels are available for free download from **BioMesh Project** site:
-                  http://www-unix.mcs.anl.gov/~csverma/BioMesh/biomesh.html
-      For further information, contact.
-                   	C.S. Verma or Paul Fischer
- 	            	Mathematics and Computer Science Division,
-                   	Argonne National Laboratory,
-                   	Argonne, IL 60439,

## Performance Results:

| Function | Time ( Seconds) |
| --- | --- |
| RBF Coefficients | 450  |
| RBF function on Box  | 90 |
| Marching Cube for Phi = 0 eval | 2 |
| Taubin Surface Fairing | 1 |
| Trimming end sections | 0.5 |
| Solution of Six heat Conduction Problems | 400  |
| Principal Chip Construction  | 1 |
| Split Chip Construction | 2 |
| Simple Chip Construction | 4 |
| Mesh Templates | 0.5 |
| Surface Smoothing | 0.001 |
| Volume Mesh | 0.0001 |

## Performance Improvement:(I) A Coarse Grid Solver:

- Marching Tetrahedra generate too many triangles that are needed to solve the heat conduction problems.
- We are not concerned with accuracy of heat conduction but in satisfying non-intersection and orthogonality criteria.
- A Coarse grid obtained from ***Independent Sets*** of Mesh Graph nodes is used to generate volume mesh.

## Performance Improvement: (II) Fast Iso-Surface Extraction:

- At each vertex, three scalar values corresponding to three solutions are stored.
- Sorted arrays are used to search for the *Active Mesh Entities* for isosurfaces.