1. Study Delaunay-Voronoi triangulation

  • Delaunay Triangulation

DelaunayTriangulation is a method used for constructing an irregulartriangulation out of a 2D point set, in a model fashion that thetriangles are well-formed.1 Regulartriangles are crucial for the accuracy of the surface representation.Before the commencement of any infrastructure design workdevelopment such as bridges, dams, roads, buildings, roads andenvironment-the affected area needs to be critically and empiricallysurveyed. It involves an intensive process of measuring the locationsof existing buildings, finding all delves and other stationaryobjects at the site. There are a number of vital detailsthat the Civil Engineer at site needs to look for in his/herobservations. The engineer should adequately survey the structure,shape and stability of the surface. Consistency or rather stabilityand composition of the existing surface are typically measured bydrilling and obtaining soil samples from the ground.

In measuring theshape of the surface, we use aerial photography, tacheometers, GPS(Global Positioning System) .2

The DelaunayTriangulation for a set P of points in a plane is atriangulation DT (P) so as point P is outside of thecircum-circle of any triangle in DT (P).3 Inanother way, all the minimum angles of all triangles are maximized inDelaunay Triangulation, which means that the process should avoid‘skinny edges’ at all cost.

The concept of Delaunay Triangulation brings about the idea of 3-Dand higher order dimensions if circum-circles appear in the notion.However, creating higher dimensions is almost non-existence orexclusive.

In Delaunay Triangulation, a collection of edges is a point setsatisfying an ‘empty circle’ property. 4 The cellnerves in a Voronoi diagram are counterpart to the triangulation in aDelaunay Triangulation.5

The data that the surveying methods yield is a random pointunprocessed data, and vast number of points on the surface has X, Yand Z coordinates that are correctly identified by the scheme. Overand over, an assumption is that the surface is 2.5 dimensional.However, for a particular 2-dimensional location (x, y) thereis only one point (x, y, z) that lies on the surface. A numberof discrepancies may arise depending on the surveying method applied.Consequently, this may also bring about the variation of number ofsample points in a given area. The accuracy and the type of sampleerror can vary as well.6For the reason, that the rawsurvey data merely consists of a set of points that lie on thesurface, another sophisticated yet clear method of constructing anestimation of the actual surface is considered necessary.

With no illustration of the surface, such responsibilities asobtaining Z of given (x, y) or scheming height line for aparticular Z, is rendered invalid. Precisely, we need a progressiverepresentation of the surface that is distinct from every point (x,y) in the region.7

In the infrastructure design, this surface illustration is usuallyreferred to as Digital Terrain Model (DTM). The initial digitalterrain models reference back to 1958 (based on the distinguishedscholar Professor Charles Miller of the Massachusetts Institute ofTechnology). Regardless, they gained popularity in mainstreaminfrastructure design industry in the 80’s. Closely relatedsoftware design to DTM is Digital Elevation Model (DEM). The digitalelevation model has a height value for each 2D point or a mean toobtain the height from the Digital Elevation Model exists.8

We have to know the operations involved in the representation so thatwe can decide on what to use. Predominantly, in a Civil Engineeringproject these operational procedures conform to the extract below:-

  • Height interpolation for random area points.

  • Solving the 2D-polyline, an intersection between surface representation and an XY-plane on an individual height line.

  • Computing the volume of the 2D-polyline.

  • Computing the volume of the 2D-polyline on a set XY area.

  • Computing the intersection between the surface representation and objects such as 3D-lines, 3D-planes or a further surface representation.

  • Actual visualization that includes polygon rendering and much more.9

Construction of a Delaunay Triangulation of any set points isachievable, provided that the set is not degenerate, i.e. ‘all thepoints may not lie in the same place or the same line’.10The Delaunay Triangulation is not exceptional when there are threeand above points on the circum-circles of any triangle in the giventriangulation.

TheDelaunay Triangulation has maximized the minimum angle of thetriangles. Hence, the minimum angle of the triangles of the DelaunayTriangulation is bigger than that of any other triangulation of thesimilar point set.

Thepoints’ set convex hull can be identified by the DelaunayTriangulation boundaries since it’s a maximal planar subdivisionthat is an important element of Delaunay Triangulation. Despite that,the triangulation will occasionally be bounded so that there aresections with concavities on areas with no point of interest or whereinput data is unavailable.

Thetypical definition of triangulation utilizes a point set in 2D. It isalso probable to produce the Delaunay Triangulation to extradimensions For instance, the three-dimensional space is called the‘DelaunayTetra-hedralyzation’where no other pointresides in the sphere defined by a tetrahedron.

Obtainingthe two-dimensional Delaunay Triangulation of a set of point issimilar to projecting the set on a 3D-paraboloid (x+y + z) 22.Afterthat, we need to find the convex hull of that point set. Extradimensions can be generated using this idea.

  • Voronoi diagram

TheDelaunay Triangulation is the twofold of the Voronoidiagram. Witha plane and a set ofpoints, the Voronoi depiction breaks the plane into sections with allthe points that are nearer to the point that distinct the region morethan any other set point. The diagram edges are equidistant from thetwo points, and the nodes are also equidistant from three oradditional points.11

Figure 1.Voronoitessellation

  • Algorithms

Thereare six essentially classes of algorithms for constructing a DelaunayTriangulation from a point set.

  1. Sweep-line algorithms that sweep the plane with a line and add edges to the triangulation as the line moves.

  2. Divide and conquer, split the point set to smaller subsets until they are trivial to triangulate and merge the subsets.

  3. Greedy algorithms that algorithms that start with one edge and incrementally construct the triangulation by adding one Delaunay triangle at a time.

  4. Refining algorithms that first form a non-Delaunay Triangulation and then refine it with edge-flips until it satisfies the Delaunay criterion.

  5. Incremental algorithms that start with a trivial triangulation and incrementally add points to it while retaining the Delaunay property.

  6. Convex hull based algorithms that take note of the fact the Delaunay Triangulation of a point set in R2 is equivalent to the convex hull of the same points set projected onto a paraboloid in R3.

Mostof the algorithms for computing Delaunay Triangulations depend onquick operations for detecting if a point is within a triangle`scircum-circle and a proper data structure for triangles and edges. Intwo dimensions, one method to identify if point&nbspD&nbspliesin the circum-circle of&nbspA,&nbspB,&nbspC&nbspisto evaluate the&nbspavailabledeterminants.12

  • Understanding Delaunay Triangulation

Themost fundamental property of a Delaunay Triangulation is the Delaunaycriterion, often referred to as the empty circum-circle criterion in2D triangulations. In a simple illustration below, the T1circum-circle is empty no interior point shown. Similarly, the T2circum-circle has no interior points as well as well as being empty.

Figure 2.DelaunayTriangulation model.

A non-DelaunayTriangulation is as illustrated below. In the diagram, we can seeinterior points V1 and point V3 where T1 and T2 circum-circlesections are occupied by the interior points.13

Figure3.Non-Delaunay Triangulation model.

DelaunayTriangulation has a ‘well shaped’ characteristic since the emptycircum-circle property there is a selection of internal large angledtriangles over small angled ones. In Figure 3 above, we can see thatthe vertices V2 and V4 have sharp edges. V1 and V2 can be joined byan edge and at the same time replacing V2 and V4 edges hence creatinga Delaunay Triangulation from a non-Delaunay Triangulation. In thisview, maximization of the minimum angles has occurred. In the nextpicture illustration, we can see the 3-D Delaunay Triangulation,composed of 2 tetrahedral objects. The unique tetrahedron exhibitsthe empty circum-sphere criterion.

Figure4.3-D Delaunay Triangulation.

MATLABis a programming language that simplifies the triangulation in termsof solving PDE (Partial Differential Equation) in a complex domain.Various functions to produce Delaunay triangulations do exist.

  • Delaunayn and Delaunay functions.

  • The Delaunay triangulation class.

2-Dand 3-D triangulations are constructed by Delaunay Triangulationclass functions. Triangulation class can support the engineering ofDelaunay algorithms. They class functions restrict with the DelaunayTriangulation. It also holds a framework to the convex hullconstruction as well as the Voronoi diagrams.

  • Rundown

TheDelaunay function exploits the fundamental triangulation of dataanalysis. The data is adequately inclusive for the sole purposeefficiency.

TheDelaunay Triangulation class is the most resourceful in triangulationbased advancements. More so, it’s exceedingly applicable in otherside based operations that include.

  • Obtaining a triangulation for triangles or a query point enclosed by tetrahedral objects.

  • Using triangulation to search a closest neighbor point.

  • Removing or adding points by triangulation modification.

  • Computation of the convex hull and Voronoi diagrams.

  • Polygonal triangulation and the choice removal of exclusive domain triangles.

  • Edges constraint.

2-Dand 3-D geometric domains are illustrated by triangulation incomputer graphics and applied in Civil Engineering. The followingfigures can help to understand the concept more accurately. 14

Figure5.A simple geographical mapping.

Thefigure 5 above shows a simple map that deploys DelaunayTriangulation. In Figure 5 below, it’s clearly noted that thetriangulation depicts a complex polygonal feature.

Figure6.Complex polygon triangulation

Inother words, the complex polygon has been condensed down to a clusterof individual yet simple triangular polygons. Geometric basedalgorithms and graphical applications utilize these polygons.

Also,the 3-D geometric domain can also be shown in a triangulation format.The figure below illustrates the concept visually.

Figure7.A 3-D geometric triangulation.

Everyindividual section of the convex hull is a triangle. The 3-D spaceillustrates the set of points.14


  1. Study mesh quality criteria

Thequality of a mesh serves a significant role in stability and accuracyof various numerical computational analysis.1Theproperties concerning mesh quality are as follows.

  • Skewness.

  • Smoothness.

  • Node point.

  • Distribution.

Studyingthe quality of the grid is so vital in any particular domain of amesh quality type than ignoring it. Such different cell types in amesh include.

  • Tetrahedral.

  • Hexahedral.

  • Polyhedral.

Withthat in mind, we can evaluate different quality criteria with regardsto mesh cell types. Such criteria are as highlighted below.

  • Face squish on polyhedral meshes.

  • Cell equal-volume skew on tri and or tetra elements.

  • Evaluating aspect ratio on all meshes.

  • Cell squish on all meshes.

  1. Aspect Ratio

Ina broader sense, the aspect ratio is described as the measure of acell stretching. It is the ratio of ‘maximum distance between thecentroid of a cell and face centroids over the minimum distancebetween the cell nodes.’2 Aflaw in the grid quality causes an error to occur in the FLUENTconsole of the sample mesh. However, it’s still possible to carrythe task in a successful way.

Figure8.Aspect ratio calculation with measurements

Wecan use a text command to check the grid quality. The box belowillustrates the command text.

grid quality


Grid Quality

Criteria quality fortriangular/mixed cells

Maximum cell squish = 4.61001e – 01

Maximum cell skewness = 4.48776e – 01

Maximum aspect ratio = 5.23830e + 00

  1. Smoothness

A truncationerror is the variationbetween the partial derivatives in the overridingequations and their uniqueestimations. Quickchanges in the cellsize between borderingcells convertinto larger truncation errors.A&nbspFLUENTallowance&nbspprovidesthe potential toadvance the smoothness by meshrefining based on avariation of cell volume or cell volume gradient.Refiningthe grid based oncell size variationis an application of Volume Adaptation. 3Changesin cell size should be steady so that accuracy and convergencemaintain in a particular mesh.Inaccurate and non-convergent meshes are corrected by gridimprovement after compromising volume adaptation in options refining.Achieved by changes in volume of the cell or/and volume changesbetween cells and their neighbouring cells.4

Refining the grid based on size of the volumeis commonly used in large cells deletion or refining the mesh on anepic scale. In this criterion, unusually bigger sized cells are improved whenever they fail to attain the set threshold mark.

The grid’s smoothness is enhanced by markingthe grid based on cell volume alterations. By looping the faces anddoing a ratio comparison of the cell neighbours to the face, we canperform a computational volume change. The Figure 9 belowdemonstrates the idea.5

Figure9.Change in cell volume-cell volume ratios.

Inthe figure above, the ratio V1/V2 and the proportion V2/V1 iscompared to a primary entry value. Whenever V2/V1 is greater than theprimary entry value, then C2 is evident to refinement.

WhereV1 is volume 1 and V2 is volume 2.

  • Volume Adaptation Approach

Inorder for me to describe the volume adaptation concept in depth, Ihave decided to use the following figures as quite simple andelaborate in the scope. The turbulent jet is computed by the mesh.6

Figure10.A jet mesh before volume adaptation.

T-Gridemployed a local marking and refining so that a fine mesh formed.However, it was coarse in a different place. Consequently, a rapidyet significant change in cell volume was produced at the jet’sedge.

Inorder to improve on mesh quality, volume adaptation should utilizethe criterion basis that the maximum volume/size of the cell shouldnot be more than 50 % that is half the total cell volume.

Figure11.A jet mesh after volume adaptation-cell volume change.

Thesharpness of the interface between the refined jet region and theneighbouring mesh is somehow reduced or absent. The final output isthe figure 11 above, done after adequately smoothing the mesh.Additionally, the minimum volume of the cell for adaptation is alsorestricted.

  • Gradient Adaptation Approach

Thegradient adaptation approach allows one to define cells or modify thegrid according to the gradient, twist, or it may also be a value ofthe preferred field variables.

In a digital solution, an adaptive gridrefinement is done so that a numerical error-free digital solutionforms. Unluckily, an undeviating error evaluation for point-insertionadaptation approach is challenging since it’s almost impossible toproduce an accurate approximation and error modification in theadapted grids. We can assume a hypothesis that a maximum erroroccurs in high gradient approaches so that the real time availablephysical elements of the changing field flow can exploit the gridadaptation method.

Additionally, this scope has covered twoapproaches for grid adaptation in FLUENT.

  1. Gradient approach: FLUENT is used in the calculation of the Euclidian norm having a gradient function of a preferred solution variable that is multiplied by the characteristic length. For instance, the gradient function in 2D assumes the subsequent layout below:

= r/2

-indicates the errorpresent.

-Cell Area.

r -Gradient volume weight.

-Undividable Laplacian of f.

f- FieldVariable.

A complete volume weighting is in conjunctionwith the typical value of the gradient volume weight. A zero valueeliminates the volume weighting. In addition to that, any values inthe range of 0-1 will make use of proportional volume weighting. Thisadvancement has an application in smooth solution problems.7

  1. Equal-value Approach: This approach does not rely on derivatives. As a substitute, the equal values of the required field variable f can regulate the adaptation. As a result, the function assumes the formula below.

= f

– Error indicator.

The cellvolume is squared or cubed to obtain a length scale value. Thepotential of additional accurate solutions is cultivated by theintroduction of length scale values so that the resolution of bothstrong and weak stress maintains in the process. However, we canchange the gradient volume weight so that the volume weighting can beregulated or chopped out. The gradient adaptation function canmanipulate any of the field variables in contouring. Physical andgeometric features of the numerical solutions are scalar functioninclusive. Hence, we can choose a cell volume field adaptation so asto reduce the spontaneous cell variations. Conventional adaptation tophysical features is velocity. However, the cell volume field ispreferred over the speed feature. 8

A supersonicflow against a cylindrical observation can be a case solution studyfor ‘steady gradient adaptation.’9

Figure12.Before Adaptation

In Figure 12above, we can see a coarse mesh regardless of the number of enoughcells to define the cylinder shape sufficiently.10Thecoarseness of the mesh ahead of the cylinder can’t resolve theshock wave&nbspthatappears in front of the cylinder. So, in this view, pressure is anappropriate variable for gradient adaptation since there will be apressure jump across the shock wave.

Figure13.Gradient Adaptation

&nbspImmediatelyafter a number of slope adjustments, the grid mesh depicts Figure 13above. In this manner, the shock wave is satisfactorily resolved.Additionally, the gradient adjustments can be computed in softwarelike dialog box panel following a set of algorithms in thebackground. The Figure 14 below illustrates the concept.

Figure14.Gradient Adaptation Panel

  1. Cell Shape

The form ofthe cell plays a vital role in the precision of the numericalelucidation. The skewness and the Aspect Ratio of a cell are some ofthe key considerations in the cell shape description. 11

  1. Flow Fluid Dependency

The simulation of a given flow field assumes aconcurrence with the impact of the cell shape, smoothness, accuracyand stability, and resolution. For instance, cells with more skewnessare accepted in weak flow regions. On the other hand, such cell typescan be so detrimental in sections with highly active flowgradients.12

The finding of ‘active flow gradients’ isusually impossible. As a result, it is also an urge to establish afirst-grade mesh quality in the whole domain region.

  • Additional concepts of the mesh quality criteria

In theresearch included here, I found out that FLUENT&nbspfunctioncanutilize triangular grids or quadrilateral cells Grid or maybe both ofthem in 2-Dimensional approaches. FLUENT is also capable of usingtetrahedral grids, hexahedral, polyhedral cells grids in 3 Dimension.The sole purpose will determine the mesh type option. When pickingthe mesh type, the issues below are carefully sought after.13

  1. Setup time:

In engineering, a lot of flow problems aresolved every day. Such problems are complex geometries affiliated.Conceptualization of complex or block-structured grids for suchcomplicated problems is prolonged or unattainable. For that reason,setup time for complex geometries is the threshold in creatingunstructured grids encompassing simple networks triangular ortetrahedral cells. Somewhat unsophisticated geometries are nottime-wasteful. Hence, setup time can disregard with eitherapproach.14

Besides, we can save time in regenerating agiven mesh. Prior creation of structured codes is essential tokeeping up with the time. Thus, it is an incentive for using complexgeometries in the FLUENT&nbspsimulationfunction.

  1. Computational expense:

In complex geometries, a triangular ortetrahedral mesh can be created with very smaller quantity cells thanthe corresponding mesh comprising of quadrilateral or hexahedralobjects. The difference is because a simple unstructured mesh permitsthe clustering of cells in pre-defined regions of the flow domain.Structured hexahedral meshes will automatically push cells to settlein areas where they have no or less importance. Unstructuredquadrilateral or hexahedral meshes place most of the benefits ofsimple meshes for relatively complex geometries. Quadrilateral orhexahedral objects allow larger aspect ratio than simply structuredcells. The skewness of the cell is ever affected when and if thesimple structured elements have a large aspect ratio causinginaccuracies and non-convergence that is unwanted. Consequently, afairly simpler geometry in which the flow is in agreement with thegeometry shape can attract a mesh of higher aspect ratio unstructuredcells. The mesh will automatically have far fewer cells over the useof unsophisticated cells.15

A lower cell count (lower than the defaultmesh) will result if a full domain of a tetrahedral mesh converts toa polyhedral mesh. The outcome is a coarser mesh, but the convergencecan maintain at a faster rate and thus minimizing on thecomputational expense.

  1. Numerical Diffusion:

A chief source of inaccuracies inextra-dimensional circumstances is numerical diffusion commonlyreferred to as ‘false diffusion’. The flow is not a realoccurrence, despite the fact that the effect on a flow calculation iscomparable to that of increasing the actual flow coefficient.


  1. Study CVT method of improving mesh quality.

The CentroidalVoronoi Tessellation provides the finest division for generatingpoints comparative to a specified density function and consequentlymaking a high-quality mesh. In my regard, I have studied thealgorithms for the building of a constrained Centroidal VoronoiTessellation-triangulation based from an original tetrahedral mesh ofa three-dimensional field. Lloyd’s iteration principle provides aplatform that can construct a controlled CVT mesh from the typicaluniversal initial mesh optimization. In addition, an affair betweenthe density function and the regulated sizing can finesse thisoperation. Mesh quality optimization also employs edge flippingmethod.1

Agood number of engineering problems require the generation ofhigh-quality mesh to solve the numerical challenges in the wholeprocess. The range of issues includes the fluid flows and structureinvestigation. Complex challenges make use of auto-unstructuredtetrahedral mesh generation for 3-D domains and thus leading to asuccessful problem-solving, in the long run. 2

Simple meshgeneration employs a number of adequate operations to solve theseproblems. They include:

  1. Voronoi-Delaunay methods.

  2. Octree technique.

  3. Advancing Front Method.

However, meshes with structuredtetrahedral and that are Delaunay based have become a commonphenomenon in finding a solution to complex challenges. This conceptis rendered possible due to their ability to perform hardarithmetical hypothesis and systemic execution.3

In the constructionof any initial mesh type, advanced optimization is given priority.Based on that idea, there is a need to employ such intensive effortstremendously for the realization of high-quality mesh generation. Forus to optimize on structured tetrahedral meshes, we need to utilizethe three algorithm based, existing categories.4

  1. Geometric optimization– This type of optimization frequently refers to as mesh smoothing. Node connectivity repositions the vertices. Mesh topology remains. One of the smoothing techniques is fine Laplacian smoothing.

  2. Topological optimization– A restricted reconnection assures that edges exchange so that there is a node connectivity change for a particular set of vertices.

  3. Vertex insertion and deletion– Mesh refining or mesh coarsening is carried out as a measure for optimization. The vertices are inserted or removed.

The techniques aboveare localized based. Their iterative combination can help to producea high-quality mesh. Nevertheless, the techniques are useless inglobal mesh enhancement. In the application, the global parameters ofthe whole mesh quality are supposedly inclusive of the meshconsistency to a pre-defined sizing field. Automatically, this refersto the function that at a given vertex, it provides the size of theedges linking the vertex. In this sense, the enhanced mesh will lackoptimization since a local operation is not enough to optimize thesample mesh. In experimental local optimization the combination isstill unable to delete all insignificant elements. In addition,global optimization techniques typically enhance the quality of themesh as well as conforming to a particular sizing. In a nutshell, wecan comfortably suggest that such challenges in high-grade meshgeneration are a stepping stone to global refinement.5

Figure 14.AComplex Model’s tetrahedral mesh and its cutting view.

  • Conclusion remarks

Briefly, theproposed tetrahedral mesh generation and optimization methods arebased on the optimization elements of the Centroidal Voronoitessellations and numerically dependable. They offer very rapidquality enhancement hence meeting the specified objective: globallyoptimizing the point distribution and the topological structure ofthe mesh through the construction of the CVT. Mesh quality furtherimproves with local optimization techniques.


  1. Prepare a MATLAB code for triangulation.

  • (please see the Appendix in this page, for this code)


Baker,and J. Timothy, ‘Mesh Generation:Art or Science?’: Progress in Aerospace Science,Elseiver Ltd, Vol. 4, Issue 1, 2009, pp. 29-63, Available from:Google Scholar Books, (retrieved on 10 February 2015).

Barzins,and Martins, ‘Mesh Quality:a Function ofGeometry’, ErrorEstimates or both?Engineering withComputers, 1999,pp. 236- 247, Available from: Google Scholar Books, (retrieved on 10February 2015).

C.B, Barber., Dobkin, D. P and Huhdanpaa, Quickhullalgorithm for convex hulls, ‘TheGeometry Centre’, ACM Trans, Mathematical software22,&lt , 1996, (retrieved 10 February2015).

Du,and M. Gungzburger, ‘Grid Generation and Optimization based onCentroidal Voronoi Tessellations’, AppliedMathematics and Computation Journal, vol.133.2, 2002, pp. 591-607.

Du, Qiang, and D. Wang, ‘Tetrahedral Mesh Generation and Optimizationbased on Centroidal Voronoi Tessellations’, InternationalJournal Numerical Methods in Engineering, vol.56.9, 2003, pp. 1355-1373.

Guibas,‘Primitives for the Manipulation of General Sub-division and theComputation of Voronoi’, ACMTransactions on Graphics,Vol. 4, Issue 2,1985, p. 107, (retrieved 10 February 2015).

Hjelle,O., and M. Daehlen, Mathematicsand Visualization Series: Triangulation andApplications,Berlin, SpringerHeidelberg, 2006, pp. 73-93.

Herva,V., ‘Using Delaunay Triangulation in Infrastructure DesignSoftware’, Masters’ Thesis, Helsinki University of Technology,2009, p. 11.

Kreveld,M.V., and M. Overmars, ‘Computational Geometry’, in Mark (ed.),Algorithms andApplications, Utrecht,Chapter 7, (retrieved 10 February 2015).

‘MeshQuality Criteria’, ANSYFLUENT, 2006, p.6.2.2, &lt,(retrieved 10 February 2015).

Nurmen,and Kimmo, ‘Digital Aerial Photographs in Municipal GeographicalInformation System’, Master’s Thesis, Helsinki University, 2002

Pasquale,G., and G. Pinelli, AGround Penetrating Radar for Soil Pattern Recognition:Geo-science andRemote Sensing Symposium Proceedings,IEEE International, 1998, Vol. 2, p. 853.

PD,Zavattieri., Dari, and Buscaglia, ‘Optimization Strategies inUnstructured Mesh Generation’, International Journal forNumerical Methods in Engineering, vol. 39, 1996, pp. 2055-2071.

Samir,J. Echaabi, and M. Hattabi, ‘Numerical Algorithm andAdaptive Meshing for Simulating the Effect of Variation inResin Transfer Molding Process’, Composites Part

B:Engineering, no. 5,2011, pp. 1015-1028.

Shewchuk,and J. Richard, ‘LectureNotes on Delaunay Mesh Generation’,The Department of Electrical Engineering and Computer Science,University of California at Berkeley, 1999,(retrieved10 February 2015).


function[Npts,Nelm,p,ne,n,nbe] = trgl6_octa (Ndiv)




%Copyright by C. Pozrikidis, 1999

%All rights reserved


%This program is to be used only under the

%stipulations of the licensing agreement



%Triangulation of the unit sphere

%into six-node quadratic triangles

%by subdividing a regular octahedron





% Ndiv …. level of discretization of icosahedron

% Nvid = 0 gives 20 elements


% Npts …. number of nodes

% Nelm …. number of surface elements


% x(i,j), y(i,j), z(i,j) ….


% Cartesian coordinates of local node j (1…6)

% on element i (1…Nelm)


% p(i,j) …. Cartesian coordinates of global node i

% where j=1,2,3, with

% x = p(i,1) y = p(i,2) z = p(i,3)


% n(i,j) …. global node lable of local node number j on element i,

% where j=1,…,6


% ne(i,j) … ne(i,1) is the number of elements touching global nodei.

% ne(i, 2:ne(i,1)+1) are the corresponding element labels


% nbe(i,j) .. label of element sharing side j of element i

% where j = 1, 2, 3



%Begin with the zeroth-level

%discretization (8 elements)


%Nodes are set manually on the unit sphere


Nelm= 8


% vertex nodes in the upper half of xz plane


x(1,1)=0.0D0% first element










x(5,1)=1.0D0% fifth element










x(6,1)=0.0D0 % sixth element










x(2,1)=-1.0D0% second element










% vertex nodes in the lower half xz plane


x(4,1)=0.0D0 % fourth element










x(8,1)=1.0D0 % eighth element










x(7,1)=0.0D0% seventh element










x(3,1)=-1.0D0% third element










%compute the mid-points of the three sides

%of the 20 first-generation elements


%midpoints are numbered 4, 5, 6



x(i,4)= 0.5D0*(x(i,1)+x(i,2))

y(i,4)= 0.5D0*(y(i,1)+y(i,2))

z(i,4)= 0.5D0*(z(i,1)+z(i,2))

x(i,5)= 0.5D0*(x(i,2)+x(i,3))

y(i,5)= 0.5D0*(y(i,2)+y(i,3))

z(i,5)= 0.5D0*(z(i,2)+z(i,3))

x(i,6)= 0.5D0*(x(i,3)+x(i,1))

y(i,6)= 0.5D0*(y(i,3)+y(i,1))

z(i,6)= 0.5D0*(z(i,3)+z(i,1))



%project the nodes onto the unit sphere




rad= sqrt(x(k,l)^2+y(k,l)^2+z(k,l)^2)

x(k,l)= x(k,l)/rad

y(k,l)= y(k,l)/rad

z(k,l)= z(k,l)/rad







%compute the local element node coordinates

%for discretization levels 1 through Ndiv



nm= 0% count the new elements arising by sub-division

%four element will be generated during each pass

forj=1:Nelm % over old elements


%assign corner points to sub-elements

%these will become the &quotnew&quot elements


nm= nm+1

xn(nm,1)= x(j,1) % first sub-element

yn(nm,1)= y(j,1)

zn(nm,1)= z(j,1)

xn(nm,2)= x(j,4)

yn(nm,2)= y(j,4)

zn(nm,2)= z(j,4)

xn(nm,3)= x(j,6)

yn(nm,3)= y(j,6)

zn(nm,3)= z(j,6)

xn(nm,4)= 0.5D0*(xn(nm,1)+xn(nm,2))

yn(nm,4)= 0.5D0*(yn(nm,1)+yn(nm,2))

zn(nm,4)= 0.5D0*(zn(nm,1)+zn(nm,2))

xn(nm,5)= 0.5D0*(xn(nm,2)+xn(nm,3))

yn(nm,5)= 0.5D0*(yn(nm,2)+yn(nm,3))

zn(nm,5)= 0.5D0*(zn(nm,2)+zn(nm,3))

xn(nm,6)= 0.5D0*(xn(nm,3)+xn(nm,1))

yn(nm,6)= 0.5D0*(yn(nm,3)+yn(nm,1))

zn(nm,6)= 0.5D0*(zn(nm,3)+zn(nm,1))

nm= nm+1

xn(nm,1)= x(j,4)% second sub-element

yn(nm,1)= y(j,4)

zn(nm,1)= z(j,4)

xn(nm,2)= x(j,2)

yn(nm,2)= y(j,2)

zn(nm,2)= z(j,2)

xn(nm,3)= x(j,5)

yn(nm,3)= y(j,5)

zn(nm,3)= z(j,5)

xn(nm,4)= 0.5D0*(xn(nm,1)+xn(nm,2))

yn(nm,4)= 0.5D0*(yn(nm,1)+yn(nm,2))

zn(nm,4)= 0.5D0*(zn(nm,1)+zn(nm,2))

xn(nm,5)= 0.5D0*(xn(nm,2)+xn(nm,3))

yn(nm,5)= 0.5D0*(yn(nm,2)+yn(nm,3))

zn(nm,5)= 0.5D0*(zn(nm,2)+zn(nm,3))

xn(nm,6)= 0.5D0*(xn(nm,3)+xn(nm,1))

yn(nm,6)= 0.5D0*(yn(nm,3)+yn(nm,1))

zn(nm,6)= 0.5D0*(zn(nm,3)+zn(nm,1))

nm= nm+1

xn(nm,1)= x(j,6) % third sub-element

yn(nm,1)= y(j,6)

zn(nm,1)= z(j,6)

xn(nm,2)= x(j,5)

yn(nm,2)= y(j,5)

zn(nm,2)= z(j,5)

xn(nm,3)= x(j,3)

yn(nm,3)= y(j,3)

zn(nm,3)= z(j,3)

xn(nm,4)= 0.5D0*(xn(nm,1)+xn(nm,2))

yn(nm,4)= 0.5D0*(yn(nm,1)+yn(nm,2))

zn(nm,4)= 0.5D0*(zn(nm,1)+zn(nm,2))

xn(nm,5)= 0.5D0*(xn(nm,2)+xn(nm,3))

yn(nm,5)= 0.5D0*(yn(nm,2)+yn(nm,3))

zn(nm,5)= 0.5D0*(zn(nm,2)+zn(nm,3))

xn(nm,6)= 0.5D0*(xn(nm,3)+xn(nm,1))

yn(nm,6)= 0.5D0*(yn(nm,3)+yn(nm,1))

zn(nm,6)= 0.5D0*(zn(nm,3)+zn(nm,1))

nm= nm+1

xn(nm,1)= x(j,4)% fourth sub-element

yn(nm,1)= y(j,4)

zn(nm,1)= z(j,4)

xn(nm,2)= x(j,5)

yn(nm,2)= y(j,5)

zn(nm,2)= z(j,5)

xn(nm,3)= x(j,6)

yn(nm,3)= y(j,6)

zn(nm,3)= z(j,6)

xn(nm,4)= 0.5D0*(xn(nm,1)+xn(nm,2)) % mid points

yn(nm,4)= 0.5D0*(yn(nm,1)+yn(nm,2))

zn(nm,4)= 0.5D0*(zn(nm,1)+zn(nm,2))

xn(nm,5)= 0.5D0*(xn(nm,2)+xn(nm,3))

yn(nm,5)= 0.5D0*(yn(nm,2)+yn(nm,3))

zn(nm,5)= 0.5D0*(zn(nm,2)+zn(nm,3))

xn(nm,6)= 0.5D0*(xn(nm,3)+xn(nm,1))

yn(nm,6)= 0.5D0*(yn(nm,3)+yn(nm,1))

zn(nm,6)= 0.5D0*(zn(nm,3)+zn(nm,1))

end% end of old-element loop


%number of elements has increased

%by a factor of four


Nelm= 4*Nelm


%relabel the new points (xn -&gt x)

%and place them in the master list




x(k,l)= xn(k,l)

y(k,l)= yn(k,l)

z(k,l)= zn(k,l)

%—project onto the unit sphere

rad= sqrt(x(k,l)^2+y(k,l)^2+z(k,l)^2)

x(k,l)= x(k,l)/rad

y(k,l)= y(k,l)/rad

z(k,l)= z(k,l)/rad

xn(k,l)= 0.0D0% zero just in case

yn(k,l)= 0.0D0

zn(k,l)= 0.0D0




end % of discretization-level loop





%Generate a list of global nodes by looping

%over all elementsand adding nodes not found

%in the current list.


%Fill in the connectivity table n(i,j)

%containing node numbers of element points 1-6



%six nodes of the first element are

%entered mannualy


p(1,1)= x(1,1)

p(1,2)= y(1,1)

p(1,3)= z(1,1)

p(2,1)= x(1,2)

p(2,2)= y(1,2)

p(2,3)= z(1,2)

p(3,1)= x(1,3)

p(3,2)= y(1,3)

p(3,3)= z(1,3)

p(4,1)= x(1,4)

p(4,2)= y(1,4)

p(4,3)= z(1,4)

p(5,1)= x(1,5)

p(5,2)= y(1,5)

p(5,3)= z(1,5)

p(6,1)= x(1,6)

p(6,2)= y(1,6)

p(6,3)= z(1,6)

n(1,1)= 1 % first node of first element is global node 1

n(1,2)= 2 % second node of first element is global node 2

n(1,3)= 3 % third node of first element is global node 3

n(1,4)= 4

n(1,5)= 5

n(1,6)= 6 % sixth node of first element is global node 6

Npts= 6


%loop over further elements


%Iflag = 0 will signal a new global node


fori=2:Nelm % loop over elements

forj=1:6 % loop over element nodes

Iflag= 0





Iflag = 1 % the node has been recorded previously

n(i,j)= k % the jth local node of element i

%is the kth global node





if(Iflag==0)% record the node

Npts= Npts+1% one more global node

p(Npts,1)= x(i,j)

p(Npts,2)= y(i,j)

p(Npts,3)= z(i,j)

n(i,j)= Npts % the jth local node of element i

%is the new global node



end% of loop over elements


%Generate connectivity table: ne(i,j)

%for elements touching global node i


% ne(i,j) … ne(i,1) is the number of elements touching

%the ith global node


% ne(i,2:ne(i,1)+1) are the corresponding element labels







ne(i,j)= 0




%loop over global nodes



ne(i,1)= 0

Icount= 1

forj=1:Nelm % loop over elements

fork=1:6 % loop over element nodes

if(abs(p(i,1)-x(j,k))&lt= eps)

if(abs(p(i,2)-y(j,k))&lt= eps)

if(abs(p(i,3)-z(j,k))&lt= eps)

Icount= Icount+1

ne(i,1)= ne(i,1)+1

ne(i,Icount)= j






end % of loop over global nodes


%Generate connectivity table nbe(i,j)


% nbe(i,j) .. label of element sharing

% the jth side of the ith element

% for j = 1, 2, 3







nbe(i,j)= 0




%loop over elements


fori=1:Nelm % loop over elements

jcount= 1

forj=4:6% loop over mid-points

fork=1:Nelm % test element

if(k~=i)% not a self-element

forl=4:6 % loop over mid-points




nbe(i,jcount)= k






end% of test element


jcount= jcount+1



end % of loop over elements


%project points p(i,j) onto the unit sphere



rad= sqrt(p(i,1)^2+p(i,2)^2+p(i,3)^2)

p(i,1)= p(i,1)/rad

p(i,2)= p(i,2)/rad

p(i,3)= p(i,3)/rad









%generate a display a network of 6-node

%triangles descending from the icosahedron


ndiv= 2

[npts,nelm,p,ne,n,nbe]= trgl6_octa (ndiv)


%plot elements subdivided into 4 triangles















px1= p(j1,1) px2 = p(j2,1) px3 = p(j3,1)

px4= p(j4,1) px5 = p(j5,1) px6 = p(j6,1)

py1= p(j1,2) py2 = p(j2,2) py3 = p(j3,2)

py4= p(j4,2) py5 = p(j5,2) py6 = p(j6,2)

pz1= p(j1,3) pz2 = p(j2,3) pz3 = p(j3,3)

pz4= p(j4,3) pz5 = p(j5,3) pz6 = p(j6,3)









11G. De Pasquale and G. Pinelli, A Ground Penetrating Radar for Soil Pattern Recognition: Geo-science and Remote Sensing Symposium Proceedings, IEEE International, 1998, Vol. 2, p. 853.

2 ibid, .p .855

3Samir, J. Echaabi, and M. Hattabi, ‘Numerical Algorithm and Adaptive Meshing for Simulating the Effect of Variation in Resin Transfer Molding Process’, Composites Part

B: Engineering, no. 5, 2011, pp. 1015-1028.

4M .V Kreveld and M. Overmars, ‘Computational Geometry’, in Mark (ed.), Algorithms and Applications, Utrecht, Chapter 7, (retrieved 10 February 2015).

5Barber, C. B, Dobkin, D. P and Huhdanpaa, Quickhull algorithm for convex hulls, ‘The Geometry Centre’, ACM Trans, Mathematical software 22,&lt , 1996, (retrieved 10 February 2015).

6Nurmen and Kimmo, ‘Digital Aerial Photographs in Municipal Geographical Information System’, Master’s Thesis, Helsinki University, 2002


8O. Hjelle and M. Daehlen, Mathematics and Visualization Series: Triangulation andApplications, Berlin, Springer Heidelberg, 2006, pp. 73-93.


10Kreveld and Overmars, loc.cit

11V. Herva, ‘Using Delaunay Triangulation in Infrastructure Design Software’, Masters’ Thesis, Helsinki University of Technology, 2009, p. 11.

12Guibas, ‘Primitives for the Manipulation of General Sub-division and the Computation of Voronoi’, ACM Transactions on Graphics, Vol. 4, Issue 2, 1985, p. 107, (retrieved 10 February 2015).

13 Shewchuk and J. Richard, ‘Lecture Notes on Delaunay Mesh Generation’, The Department of Electrical Engineering and Computer Science, University of California at Berkeley, 1999, (retrieved 10 February 2015).


21 ‘Mesh Quality Criteria’, ANSY FLUENT, 2006, p. 6.2.2, &lt, (retrieved 10 February 2015).





6ibid. , p. 26.8.2

7 ‘Mesh Quality Criteria’, loc. cit.




11‘Mesh Quality Criteria’, op.cit, 6.2.2

12Barzins and Martins, ‘Mesh Quality: a Function of Geometry’, Error Estimates or both? Engineering with Computers, 1999, pp. 236- 247, Available from: Google Scholar Books, (retrieved on 10 February 2015).

13Baker and J. Timothy, ‘Mesh Generation: Art or Science?’ : Progress in Aerospace Science, Elseiver Ltd, Vol. 4, Issue 1, 2009, pp. 29-63, Available from: Google Scholar Books, (retrieved on 10 February 2015).


31Du , Qiang, and D. Wang, ‘Tetrahedral Mesh Generation and Optimization based on Centroidal Voronoi Tessellations’, International Journal Numerical Methods in Engineering, vol. 56.9, 2003, pp. 1355-1373.


3Du and M. Gungzburger, ‘Grid Generation and Optimization based on Centroidal Voronoi Tessellations’, Applied Mathematics and Computation Journal, vol. 133.2, 2002, pp. 591-607.


5Zavattieri, PD. Dari, and Buscaglia, ‘Optimization Strategies in Unstructured Mesh Generation’, International Journal for Numerical Methods in Engineering, vol. 39, 1996, pp. 2055-2071.

Related Posts

© All Right Reserved