5
$\begingroup$

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 = 
  AssociationThread[
   Range[Length[triangleList]] -> Values[triangleList]];


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

 tempBuffer = 
  AssociationThread[Range[Length[tempBuffer]] -> Values[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}, 
    Heads -> False][[All, 1, 1]];
 tempBuffer = KeyDrop[tempBuffer, rows];

 tempBuffer = 
  AssociationThread[Range[Length[tempBuffer]] -> Values[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 = 
  AssociationThread[
   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 = 
  AssociationThread[Range[Length[triangles]] -> Values[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"]
$\endgroup$
10
$\begingroup$

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.

If you want a book then "Delaunay Mesh Generation", by Cheng, Dey and Shewchuck (CRC Press) is a starting point.

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.

I am looking forward to see your mesher ;-)

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.