# Where to get the detailed code of MMA for grid mesh

I know DelaunayMesh can triangulate the specified area, but I can't see the code and partition idea of this built-in function. I want to know how to get the source code of grid generation or where to get the detailed code of MMA grid generation.

I've done a related code before, but I want to get more detailed code on Meshing.

isInCircleRight[circleCenter_?ListQ, R_, point_?ListQ] :=
Module[{maxPointInCircleX},
maxPointInCircleX = circleCenter[[1]] + R;
If[point[[1]] >= maxPointInCircleX, True, False]]
(*SetAttributes[drop,HoldFirst];
drop[as_?AssociationQ,list_?ListQ]:=Delete[as,Transpose@{list}]*)
x = {{211.18}, {211.28}, {211.1}, {211.38}, {211.45}, {211.23}, \
{211.31}, {210.98}, {211.13}, {211.39}, {211.06}, {211.51}, \
{211.19}, {211.26}, {211.32}, {211.1}, {211.01}, {211.39}, \
{211.48}, {211.17}, {211.21}, {211.25}, {211.29}, {210.91}, \
{211.57}, {211.15}, {211.19}, {211.23}, {211.27}, {211.31}, \
{211.34}, {211.06}, {211.47}, {211.13}, {211.17}, {211.21}, \
{211.25}, {211.29}, {211.32}, {211.53}, {210.96}, {211.11}, \
{211.15}, {211.19}, {211.23}, {211.27}, {211.3}, {211.34}, \
{211.38}, {211.04}};
y = {{-47.97}, {-47.97}, {-47.98}, {-47.98}, {-48.}, {-48.01}, \
{-48.01}, {-48.02}, {-48.02}, {-48.02}, {-48.03}, {-48.03}, \
{-48.04}, {-48.04}, {-48.05}, {-48.06}, {-48.07}, {-48.07}, \
{-48.07}, {-48.08}, {-48.08}, {-48.08}, {-48.08}, {-48.09}, \
{-48.09}, {-48.1}, {-48.1}, {-48.1}, {-48.1}, {-48.1}, \
{-48.1}, {-48.11}, {-48.11}, {-48.12}, {-48.12}, {-48.12}, \
{-48.12}, {-48.12}, {-48.12}, {-48.12}, {-48.14}, {-48.14}, \
{-48.14}, {-48.14}, {-48.14}, {-48.14}, {-48.14}, {-48.14}, \
{-48.14}, {-48.15}};

data = MapThread[List, {x // Flatten, y // Flatten}];

Xmax = Max[x]; Xmin = Min[x]; Ymax = Max[y]; Ymin = Min[y];

upPoint = {Xmin - 0.01 + (Xmax + 0.02 - (Xmin - 0.01))/2,
Ymax + 0.02 + 0.4};
h = upPoint[[2]] - (Ymax + 0.02);
w = upPoint[[1]] - (Xmin - 0.01);
arc = ArcTan[h/w];
h1 = upPoint[[2]] - (Ymin - 0.01);
w1 = h1/arc;
leftPoint = {upPoint[[1]] - w1, Ymin - 0.01 - 0.1};
rightPoint = {upPoint[[1]] + w1, Ymin - 0.01 - 0.1};

Graphics[{Line[{{leftPoint[[1]], leftPoint[[2]]}, {rightPoint[[1]],
rightPoint[[2]]}}],
Line[{{rightPoint[[1]], rightPoint[[2]]}, {upPoint[[1]],
upPoint[[2]]}}],
Line[{{upPoint[[1]], upPoint[[2]]}, {leftPoint[[1]],
leftPoint[[2]]}}]}, Axes -> True, AxesOrigin -> {210.5, -48.3}]

pointSet = SortBy[data, {#[[1]], -#[[2]]} &];
pointNum = Length[pointSet];

triangleList = <||>;
AssociateTo[triangleList, 1 -> <|1 -> leftPoint|>];
AssociateTo[triangleList[1], 2 -> rightPoint];
AssociateTo[triangleList[1], 3 -> upPoint];

tempBuffer = <||>;
triangles = <||>;
DtriangleNum = 0;
For[i = 1, i <= pointNum, i++,
triangleNum = Length[triangleList];
tempBuffer = <||>;
(*ps={};*)
dellists = {};

For[j = 1, j <= triangleNum, j++,
circ = Circumsphere[Values[triangleList[j]]];
{circleCenter, R} = {First@circ, Last@circ};
inCircle = Element[pointSet[[i]], Disk[circleCenter, R]];
If[inCircle === True,
temp = <||>;
temp[1] = {triangleList[j, 1], triangleList[j, 2]};
temp[2] = {triangleList[j, 1], triangleList[j, 3]};
temp[3] = {triangleList[j, 2], triangleList[j, 3]};
len = Length[tempBuffer];
tempBuffer = Join[tempBuffer, KeyMap[# + len &, temp]];

dellists = {dellists, j} // Flatten;,
inRight = isInCircleRight[circleCenter, R, pointSet[[i]]];
If[inRight === True,
DtriangleNum = DtriangleNum + 1;

triangles[DtriangleNum] = triangleList[j];]]
];

triangleList = KeyDrop[triangleList, dellists];
triangleList =
Range[Length[triangleList]] -> Values[triangleList]];

tempBuffer = Map[Append[#, 0] &, tempBuffer];

tempBuffer =
If[Length[tempBuffer] != 0,
For[k = 1, k <= (Length[tempBuffer] - 1), k++,
For[ s = k + 1, s <= Length[tempBuffer], s++,
If[(tempBuffer[k] === tempBuffer[s]),
tempBuffer[k] = ReplacePart[tempBuffer[k], 3 -> 1];
tempBuffer[s] = ReplacePart[tempBuffer[s], 3 -> 1];]]]

];

rows = Position[tempBuffer, _?(Last[##] != 0 &), {1},
tempBuffer = KeyDrop[tempBuffer, rows];

tempBuffer =

triangleNum = Length[triangleList];
(*Print[triangleNum];*)
For[w = 1, w <= Length[tempBuffer], w = w + 1,

triangleList[triangleNum + w] = <||>;
triangleList[triangleNum + w, 1] = pointSet[[i]];
triangleList[triangleNum + w, 2] = tempBuffer[[w, 1]];
triangleList[triangleNum + w, 3] = tempBuffer[[w, 2]];];
triangleList =
Range[Length[triangleList]] -> Values[triangleList]];

]

triangles =
Join[triangles, KeyMap[# + Length[triangles] &, triangleList]];

n = Length[triangles];
del2list = {};
For[i = 1, i <= n, i++, flag = 0;
row = (Position[triangles[i][[;; , 1]], leftPoint[[1]]][[All, 1,
1]]);
If[row != {},
If[(First@triangles[i][[row, 2]]) === leftPoint[[2]], flag = 1,
flag = 0]];
row = (Position[triangles[i][[;; , 1]], rightPoint[[1]]][[All, 1,
1]]); If[row != {},
If[(First@triangles[i][[row, 2]]) === rightPoint[[2]], flag = 1,
flag = 0]];

row = (Position[triangles[i][[;; , 1]], upPoint[[1]]][[All, 1, 1]]);
If[row != {},
If[(First@triangles[i][[row, 2]]) === upPoint[[2]], flag = 1,
flag = 0]]; If[flag === 1, del2list = {del2list, i} // Flatten]];
triangles = KeyDrop[triangles, del2list];

triangles =

n = Length[triangles];

Graphics@Table[
Line[Append[Values@triangles[i], First[(Values@triangles[i])]]], {i,
1, n}]
DelaunayMesh[data,
MeshCellStyle -> {{1, All} -> {Thick, Red}, {0, All} -> {Black,
PointSize[0.02]}}, PlotTheme -> "Lines"]


Mathematica uses Triangle as a 2D mesh generator. The interface to triangle is in TriangleLink. Triangle comes with source code and you can find that here. Triangle is also used as a basis for DelaunayMesh, if I am not mistaken.
Should you want to write your own mesh generator you can easily make use of that in the finite element package. There are also other packages that make their own mesh for FEM, for example Pinti's MeshTools or DistMesh, which is a simple but high quality mesh generator that is implemented in the FEMAddOns. Both of these return an ElementMesh. That can be used for NDSolve, NIntegrate and NDEigensystem etc.