I'm working in a finite element mesh generator. I built this function which generates an 8 node mesh (polynomials of order 2) without any interior node:
(*Generate Grid Mesh of dimensions axb with nx divisions in x and ny \
divisions in y*)
GenerateGridMesh[aa_, bb_, nx_, ny_, order_] :=
Block[{x = 0., y = 0., dx, dy, meshnodes = {}, i, j,
meshtopology = {}, allcoords, k, topolsz, l, data, c, a, b},
k = 0;
meshnodes = {};
dx = aa/(2 nx);
dy = bb/(2 ny);
For[i = 1, i <= 2 ny + 1, i++,
If[OddQ[i] == True,
For[j = 1, j <= 2 nx + 1, j++,
AppendTo[meshnodes, {x, y}];
x += dx ;
];
,
For[k = 1, k <= nx + 1, k++,
AppendTo[meshnodes, {x, y}];
x += 2 dx ;
];
];
x = 0;
y += dy;
];
meshtopology = {};
b = 0;
a = 1;
l = 0;
c = 3 nx + 2;
For[i = 1, i <= ny, i++,
For[j = 1, j <= nx, j++,
data = {a, a + 2, 3 nx + 4 + a, 3 nx + 3 + b, a + 1,
2 nx + 3 + l, 3 nx + 4 + b, 2 nx + 2 + l};
AppendTo[meshtopology, data];
a += 2;
b += 2;
l += 1;
];
l = 3 nx + 2 + c (i - 1);
a = 3 nx + 3 + c (i - 1);
b = 3 nx + 2 + c (i - 1);
];
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, Length[meshtopology[[1]]]}];
{allcoords, meshnodes, meshtopology}
];
(*Generates graphics to visualize mesh and nodes*)
GenerateGraphics[nodes_, topology_, order_] :=
Block[{meshvis, nodevis, v},
If[order == 1, v = {1, 2, 3, 4}, v = {1, 5, 2, 6, 3, 7, 4, 8}];
meshvis =
Graphics[{FaceForm[], EdgeForm[Black],
GraphicsComplex[nodes, Polygon[topology[[All, v]]]]}];
(*nodevis=Graphics[{MapIndexed[Text[#2[[1]],#1,{-1,1}]&,
nodes],{Blue,Point[nodes]}}];*)
nodevis =
Graphics[{MapIndexed[
Style[Text[#2[[1]], #1, {-1.8, 1.8}], FontSize -> 9] &,
nodes], {PointSize[Large], Black, Point[nodes]}}];
{meshvis, nodevis}
];
L = 5;
h = 5;
nx = 2;
ny = 2;
order = 2;
{allcoords, meshnodes, meshtopology} =
GenerateGridMesh[L, h, nx, ny,
order];(*Generate finite element mesh*)
{meshvis, nodevis} =
GenerateGraphics[meshnodes, meshtopology,
order];(*Generates graphics to visualize mesh*)
Show[meshvis, nodevis, AspectRatio -> Automatic, ImageSize -> Large]
which results in the following mesh:
[![enter image description here][1]][1]
I want to build a generic mesh generator for any polynomial order. Here is an example of what I need:
L = 5;
h = 5;
x = 0;
y = 0;
nx = 2;
ny = 2;
order = 3;
meshnodes = {};
dx = L/(nx order);
dy = h/(ny order);
For[irow = 1, irow <= order nx + 1, irow++,
For[icol = 1, icol <= order ny + 1, icol++,
AppendTo[meshnodes, {x, y}];
If[OddQ[Mod[irow, 3]] == True,
x += dx ;
,
x += 3 dx ;
icol += 2;
];
];
y += dy;
x = 0;
];
meshtopology = {{1, 4, 17, 14, 2, 9, 16, 11, 3, 12, 15, 8}, {4, 7, 20,
17, 5, 10, 19, 12, 6, 13, 18, 9}, {14, 17, 30, 27, 15, 22, 29,
24, 16, 25, 28, 21}, {17, 20, 33, 30, 18, 23, 32, 25, 19, 26, 31,
22}};
{meshvis, nodevis} =
GenerateGraphics[meshnodes, meshtopology,
order];(*Generates graphics to visualize mesh*)
Show[meshvis, nodevis, AspectRatio -> Automatic, ImageSize -> Large]
[![enter image description here][2]][2]
I need this to be created automatically for any dimensions of L and h, and for any node quantity.
Below is an example of a code that generates a curved mesh (nine noded elements, not serendipity).
GenerateGridMesh[R0_, RE_, nx_, ny_, order_] :=
Block[{x = 0., y = 0., dx, dy, meshnodes, i, j, meshtopology = {},
allcoords, k, l},
meshnodes = {};
k = 0;
(*meshnodes=Flatten[Table[Table[{R Cos[\[Theta]],
R Sin[\[Theta]]},{R,R0,RE,(RE-R0)/(nx order-2)}],{\[Theta],0,Pi/2,
Pi/2 /(ny order-2)}],1]//N;*)
r = (RE/R0)^(1/(-2 + nx order));
meshnodes =
Flatten[Table[
Table[{ R0 r^(n - 1) Cos[\[Theta]],
R0 r^(n - 1) Sin[\[Theta]]}, {n, 1.,
nx order - 1}], {\[Theta], 0, Pi/2, Pi/2 /(ny order - 2)}],
1] // N;
k = 0;
For[i = 1 , i < ny, i++,
l = 1;
For[j = 1, j < nx, j++,
(*AppendTo[meshtopology,{j+k,j+2+k,4 nx+j+k,4 nx-2+j+k,j+1+k,j+1+
nx 2+k,j+nx 4-1+k,2 nx+ j-1+k,2 nx+ j+k}];*)
AppendTo[
meshtopology, {l + k, l + 2 + k, 4 nx + l + k, 4 nx - 2 + l + k,
l + 1 + k, l + 1 + nx 2 + k, l + nx 4 - 1 + k, 2 nx + l - 1 + k,
2 nx + l + k}];
l += 2;
];
k += 4 nx - 2;
];
If[order == 2,
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, 9}];
,
allcoords =
Table[meshnodes[[meshtopology[[i, j]]]], {i, 1,
Length[meshtopology]}, {j, 1, 4}];
];
{allcoords, meshnodes, meshtopology}
]
GenerateGraphics[nodes_, topology_, order_] :=
Block[{meshvis, nodevis},
If[order == 2,
meshvis =
Graphics[{FaceForm[], EdgeForm[Blue],
GraphicsComplex[nodes,
Polygon[topology[[All, {1, 5, 2, 6, 3, 7, 4, 8}]]]]}];
,
meshvis =
Graphics[{FaceForm[], EdgeForm[Blue],
GraphicsComplex[nodes,
Polygon[topology[[All, {1, 2, 3, 4}]]]]}];
];
nodevis =
Graphics[{MapIndexed[Text[#2[[1]], #1, {-1, 1}] &, nodes], {Blue,
Point[nodes]}}];
{meshvis, nodevis}
];
interpolatingQuadBezierCurve[pts_List] /; Length[pts] == 3 :=
BezierCurve[{pts[[1]], 1/2 (-pts[[1]] + 4 pts[[2]] - pts[[3]]),
pts[[3]]}];
interpolatingQuadBezierCurve[ptslist_List] :=
interpolatingQuadBezierCurve /@ ptslist;
interpolatingQuadBezierCurveComplex[coords_, indices_] :=
interpolatingQuadBezierCurve[Map[coords[[#]] &, indices]]
GenerateGraphics[nodes_, topology_] := Block[{meshvis, nodevis},
nodevis =
Graphics[{MapIndexed[
Style[Text[#2[[1]], #1, {-1.8, 1.8}], FontSize -> 12] &,
nodes], {PointSize[Large], Black, Point[nodes]}}];
Show[nodevis]]
order = 2;
serendipity = False;
{allcoords, nnodes, topol} = GenerateGridMesh[100, 200, 5, 4, order];
linestopology = Flatten[Table[
{{topol[[i]][[1]], topol[[i]][[5]], topol[[i]][[2]]},
{topol[[i]][[2]], topol[[i]][[6]], topol[[i]][[3]]},
{topol[[i]][[3]], topol[[i]][[7]], topol[[i]][[4]]},
{topol[[i]][[4]], topol[[i]][[8]], topol[[i]][[1]]}
}, {i, 1, Length[topol]}], 1];
Show[GenerateGraphics[nnodes, topol],
Graphics[interpolatingQuadBezierCurveComplex[nnodes, linestopology]],
ImageSize -> Automatic]
[![enter image description here][3]][3]
[1]: https://i.sstatic.net/UgBFK.png
[2]: https://i.sstatic.net/NKxIs.png
[3]: https://i.sstatic.net/gQfL3.png