TRIANGULATION OF A DOMAIN BY USE OF MATLAB FUNCTIONS
1
TRIANGULATIONOF A DOMAIN BY USE OF MATLAB FUNCTIONS

Study DelaunayVoronoi 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 wellformed.^{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 andenvironmentthe 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 thecircumcircle 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 3Dand higher order dimensions if circumcircles appear in the notion.However, creating higher dimensions is almost nonexistence 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 2dimensional 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.^{6}For 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 2Dpolyline, an intersection between surface representation and an XYplane on an individual height line.

Computing the volume of the 2Dpolyline.

Computing the volume of the 2Dpolyline on a set XY area.

Computing the intersection between the surface representation and objects such as 3Dlines, 3Dplanes 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’.^{10}The Delaunay Triangulation is not exceptional when there are threeand above points on the circumcircles 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 threedimensional space is called the‘DelaunayTetrahedralyzation’where no other pointresides in the sphere defined by a tetrahedron.
Obtainingthe twodimensional Delaunay Triangulation of a set of point issimilar to projecting the set on a 3Dparaboloid (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.

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

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

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

Refining algorithms that first form a nonDelaunay Triangulation and then refine it with edgeflips until it satisfies the Delaunay criterion.

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

Convex hull based algorithms that take note of the fact the Delaunay Triangulation of a point set in R^{2}^{ }is equivalent to the convex hull of the same points set projected onto a paraboloid in R^{3}.
Mostof the algorithms for computing Delaunay Triangulations depend onquick operations for detecting if a point is within a triangle`scircumcircle and a proper data structure for triangles and edges. Intwo dimensions, one method to identify if point D liesin the circumcircle of A, B, C isto evaluate the availabledeterminants.^{12}

Understanding Delaunay Triangulation
Themost fundamental property of a Delaunay Triangulation is the Delaunaycriterion, often referred to as the empty circumcircle criterion in2D triangulations. In a simple illustration below, the T1circumcircle is empty no interior point shown. Similarly, the T2circumcircle has no interior points as well as well as being empty.
Figure 2.DelaunayTriangulation model.
A nonDelaunayTriangulation is as illustrated below. In the diagram, we can seeinterior points V1 and point V3 where T1 and T2 circumcirclesections are occupied by the interior points.^{13}
Figure3.NonDelaunay Triangulation model.
DelaunayTriangulation has a ‘well shaped’ characteristic since the emptycircumcircle 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 nonDelaunay Triangulation. In thisview, maximization of the minimum angles has occurred. In the nextpicture illustration, we can see the 3D Delaunay Triangulation,composed of 2 tetrahedral objects. The unique tetrahedron exhibitsthe empty circumsphere criterion.
Figure4.3D 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.
2Dand 3D 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.
2Dand 3D 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 3D geometric domain can also be shown in a triangulation format.The figure below illustrates the concept visually.
Figure7.A 3D geometric triangulation.
Everyindividual section of the convex hull is a triangle. The 3D 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.^{1}Theproperties 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 equalvolume skew on tri and or tetra elements.

Evaluating aspect ratio on all meshes.

Cell squish on all meshes.

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
GridQuality
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

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 FLUENTallowance providesthe 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. ^{3}Changesin cell size should be steady so that accuracy and convergencemaintain in a particular mesh.Inaccurate and nonconvergent 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 volumecell 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.
TGridemployed 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 adaptationcell 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 errorfree digital solutionforms. Unluckily, an undeviating error evaluation for pointinsertionadaptation 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.

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 01 will make use of proportional volume weighting. Thisadvancement has an application in smooth solution problems.^{7}

Equalvalue 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.^{10}Thecoarseness of the mesh ahead of the cylinder can’t resolve theshock wave thatappears 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
 Immediatelyafter 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

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}

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 afirstgrade mesh quality in the whole domain region.

Additional concepts of the mesh quality criteria
In theresearch included here, I found out that FLUENT functioncanutilize triangular grids or quadrilateral cells Grid or maybe both ofthem in 2Dimensional 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}

Setup time:
In engineering, a lot of flow problems aresolved every day. Such problems are complex geometries affiliated.Conceptualization of complex or blockstructured 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 nottimewasteful. 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 simulationfunction.

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 predefined 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 nonconvergence 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.

Numerical Diffusion:
A chief source of inaccuracies inextradimensional 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.
^{2}

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 highquality mesh. In my regard, I have studied thealgorithms for the building of a constrained Centroidal VoronoiTessellationtriangulation based from an original tetrahedral mesh ofa threedimensional 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 ofhighquality mesh to solve the numerical challenges in the wholeprocess. The range of issues includes the fluid flows and structureinvestigation. Complex challenges make use of autounstructuredtetrahedral mesh generation for 3D domains and thus leading to asuccessful problemsolving, in the long run. ^{2}
Simple meshgeneration employs a number of adequate operations to solve theseproblems. They include:

VoronoiDelaunay methods.

Octree technique.

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 highquality mesh generation. Forus to optimize on structured tetrahedral meshes, we need to utilizethe three algorithm based, existing categories.^{4}

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.

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

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 highquality 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 predefined 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 highgrade 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.
^{3}

Prepare a MATLAB code for triangulation.

(please see the Appendix in this page, for this code)
Bibliography
Baker,and J. Timothy, ‘Mesh Generation:Art or Science?’: Progress in Aerospace Science,Elseiver Ltd, Vol. 4, Issue 1, 2009, pp. 2963, 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,<http://www.qhull.org/.> , 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. 591607.
Du, Qiang, and D. Wang, ‘Tetrahedral Mesh Generation and Optimizationbased on Centroidal Voronoi Tessellations’, InternationalJournal Numerical Methods in Engineering, vol.56.9, 2003, pp. 13551373.
Guibas,‘Primitives for the Manipulation of General Subdivision 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. 7393.
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, <http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1093.htm>,(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:Geoscience 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. 20552071.
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. 10151028.
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).
Appendix
function[Npts,Nelm,p,ne,n,nbe] = trgl6_octa (Ndiv)
%
%FDLIB BEMLIB
%
%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 sixnode quadratic triangles
%by subdividing a regular octahedron
%
%SYMBOLS:
%—
%
% 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 zerothlevel
%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
y(1,1)=0.0D0
z(1,1)=1.0D0
x(1,2)=1.0D0
y(1,2)=0.0D0
z(1,2)=0.0D0
x(1,3)=0.0D0
y(1,3)=1.0D0
z(1,3)=0.0D0
%—
x(5,1)=1.0D0% fifth element
y(5,1)=0.0D0
z(5,1)=0.0D0
x(5,2)=0.0D0
y(5,2)=0.0D0
z(5,2)=1.0D0
x(5,3)=0.0D0
y(5,3)=1.0D0
z(5,3)=0.0D0
%—
x(6,1)=0.0D0 % sixth element
y(6,1)=0.0D0
z(6,1)=1.0D0
x(6,2)=1.0D0
y(6,2)=0.0D0
z(6,2)=0.0D0
x(6,3)=0.0D0
y(6,3)=1.0D0
z(6,3)=0.0D0
%—
x(2,1)=1.0D0% second element
y(2,1)=0.0D0
z(2,1)=0.0D0
x(2,2)=0.0D0
y(2,2)=0.0D0
z(2,2)=1.0D0
x(2,3)=0.0D0
y(2,3)=1.0D0
z(2,3)=0.0D0
%—
% vertex nodes in the lower half xz plane
%—
x(4,1)=0.0D0 % fourth element
y(4,1)=0.0D0
z(4,1)=1.0D0
x(4,2)=0.0D0
y(4,2)=1.0D0
z(4,2)=0.0D0
x(4,3)=1.0D0
y(4,3)=0.0D0
z(4,3)=0.0D0
%—
x(8,1)=1.0D0 % eighth element
y(8,1)=0.0D0
z(8,1)=0.0D0
x(8,2)=0.0D0
y(8,2)=1.0D0
z(8,2)=0.0D0
x(8,3)=0.0D0
y(8,3)=0.0D0
z(8,3)=1.0D0
%—
x(7,1)=0.0D0% seventh element
y(7,1)=0.0D0
z(7,1)=1.0D0
x(7,2)=0.0D0
y(7,2)=1.0D0
z(7,2)=0.0D0
x(7,3)=1.0D0
y(7,3)=0.0D0
z(7,3)=0.0D0
%—
x(3,1)=1.0D0% third element
y(3,1)=0.0D0
z(3,1)=0.0D0
x(3,2)=0.0D0
y(3,2)=1.0D0
z(3,2)=0.0D0
x(3,3)=0.0D0
y(3,3)=0.0D0
z(3,3)=1.0D0
%–
%compute the midpoints of the three sides
%of the 20 firstgeneration elements
%
%midpoints are numbered 4, 5, 6
%–
fori=1:Nelm
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))
end
%—
%project the nodes onto the unit sphere
%—
fork=1:Nelm
forl=1:6
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
end
end
%—
if(Ndiv>0)
%—
%—
%compute the local element node coordinates
%for discretization levels 1 through Ndiv
%—
fori=1:Ndiv
nm= 0% count the new elements arising by subdivision
%four element will be generated during each pass
forj=1:Nelm % over old elements
%—
%assign corner points to subelements
%these will become the "new" elements
%—
nm= nm+1
xn(nm,1)= x(j,1) % first subelement
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 subelement
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 subelement
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 subelement
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 oldelement loop
%–
%number of elements has increased
%by a factor of four
%–
Nelm= 4*Nelm
%—
%relabel the new points (xn > x)
%and place them in the master list
%—
fork=1:Nelm
forl=1:6
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
end
%–
end % of discretizationlevel loop
%
end
%
%—
%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 16
%—
%
%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
fork=1:Npts
if(abs(x(i,j)p(k,1))<=eps)
if(abs(y(i,j)p(k,2))<=eps)
if(abs(z(i,j)p(k,3))<=eps)
Iflag = 1 % the node has been recorded previously
n(i,j)= k % the jth local node of element i
%is the kth global node
end
end
end
end
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
end
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
%–
%—
%initialize
%—
fori=1:Npts
forj=1:7
ne(i,j)= 0
end
end
%—
%loop over global nodes
%—
fori=1:Npts
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))<= eps)
if(abs(p(i,2)y(j,k))<= eps)
if(abs(p(i,3)z(j,k))<= eps)
Icount= Icount+1
ne(i,1)= ne(i,1)+1
ne(i,Icount)= j
end
end
end
end
end
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
%–
%—
%initialize
%—
fori=1:Nelm
forj=1:3
nbe(i,j)= 0
end
end
%—
%loop over elements
%—
fori=1:Nelm % loop over elements
jcount= 1
forj=4:6% loop over midpoints
fork=1:Nelm % test element
if(k~=i)% not a selfelement
forl=4:6 % loop over midpoints
if(abs(x(i,j)x(k,l))<=eps)
if(abs(y(i,j)y(k,l))<=eps)
if(abs(z(i,j)z(k,l))<=eps)
nbe(i,jcount)= k
end
end
end
end
end
end% of test element
if(nbe(i,jcount)~=0)
jcount= jcount+1
end
end
end % of loop over elements
%—
%project points p(i,j) onto the unit sphere
%—
fori=1:Npts
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
end
%
%done
%
return
closeall
clearall
%==============
%generate a display a network of 6node
%triangles descending from the icosahedron
%==============
ndiv= 2
[npts,nelm,p,ne,n,nbe]= trgl6_octa (ndiv)
%
%plot elements subdivided into 4 triangles
%
figure(1)
holdon
axisequal
boxon
axisoff
%plot3(p(:,1),p(:,2),p(:,3),`ko`)
fori=1:nelm
j1=n(i,1)
j2=n(i,2)
j3=n(i,3)
j4=n(i,4)
j5=n(i,5)
j6=n(i,6)
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)
%patch([px1,px4,px6],[py1,py4,py6],[pz1,pz4,pz6],`w`)
%patch([px2,px5,px4],[py2,py5,py4],[pz2,pz5,pz4],`w`)
%patch([px3,px6,px5],[py3,py6,py5],[pz3,pz6,pz5],`w`)
%patch([px4,px5,px6],[py4,py5,py6],[pz4,pz5,pz6],`w`)
patch([px1,px4,px2,px5,px3,px6,px1],…
[py1,py4,py2,py5,py3,py6,py1],…
[pz1,pz4,pz2,pz5,pz3,pz6,pz1],`w`)
end
1^{1}G. De Pasquale and G. Pinelli, A Ground Penetrating Radar for Soil Pattern Recognition: Geoscience and Remote Sensing Symposium Proceedings, IEEE International, 1998, Vol. 2, p. 853.
^{2 }ibid, .p .855
^{3}Samir, 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. 10151028.
^{4}M^{ }.V Kreveld and M. Overmars, ‘Computational Geometry’, in Mark (ed.), Algorithms and Applications, Utrecht, Chapter 7, (retrieved 10 February 2015).
^{5}Barber, C. B, Dobkin, D. P and Huhdanpaa, Quickhull algorithm for convex hulls, ‘The Geometry Centre’, ACM Trans, Mathematical software 22,<http://www.qhull.org/.> , 1996, (retrieved 10 February 2015).
^{6}Nurmen and Kimmo, ‘Digital Aerial Photographs in Municipal Geographical Information System’, Master’s Thesis, Helsinki University, 2002
^{7}ibid.
^{8}O. Hjelle and M. Daehlen, Mathematics and Visualization Series: Triangulation andApplications, Berlin, Springer Heidelberg, 2006, pp. 7393.
^{9}ibid.
^{10}Kreveld and Overmars, loc.cit
^{11}V. Herva, ‘Using Delaunay Triangulation in Infrastructure Design Software’, Masters’ Thesis, Helsinki University of Technology, 2009, p. 11.
^{12}Guibas, ‘Primitives for the Manipulation of General Subdivision 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).
^{14}ibid.
2^{1 }‘Mesh Quality Criteria’, ANSY FLUENT, 2006, p. 6.2.2, < http://aerojet.engr.ucdavis.edu/fluenthelp/html/ug/node1093.htm>, (retrieved 10 February 2015).
^{2}ibid.
^{3}ibid.
^{4}ibid.
^{5}ibid.
^{6}ibid. , p. 26.8.2
^{7 }‘Mesh Quality Criteria’, loc. cit.
^{8}ibid.
^{9}ibid.
^{10}ibid.
^{11}‘Mesh Quality Criteria’, op.cit, 6.2.2
^{12}Barzins 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).
^{13}Baker and J. Timothy, ‘Mesh Generation: Art or Science?’ : Progress in Aerospace Science, Elseiver Ltd, Vol. 4, Issue 1, 2009, pp. 2963, Available from: Google Scholar Books, (retrieved on 10 February 2015).
^{ 14}ibid.
3^{1}Du , 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. 13551373.
^{2}ibid.
^{3}Du and M. Gungzburger, ‘Grid Generation and Optimization based on Centroidal Voronoi Tessellations’, Applied Mathematics and Computation Journal, vol. 133.2, 2002, pp. 591607.
^{4}ibid.
^{5}Zavattieri, PD. Dari, and Buscaglia, ‘Optimization Strategies in Unstructured Mesh Generation’, International Journal for Numerical Methods in Engineering, vol. 39, 1996, pp. 20552071.