Help to create a 2D mesh generator (FEM)

Question

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:

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]

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]

Answer 1

I am not sure, if this fully answers your question, but you should be able to work from here. I tried to explain my code with the comments above each for-loop and I maintained your overall structure.

(*Generate Grid Mesh of dimensions axb with nx divisions in x and ny \
divisions in y*)
GenerateGridMesh[aa_, bb_, nx_, ny_, p_] := 
  Block[{x = 0., y = 0., dx, dy, meshnodes = {}, i, j, 
    meshtopology = {}, allcoords, k, topolsz, l, data, c, a, b}, k = 0;
   meshnodes = {};
   (*determine the distance between each node*)
   dx = aa/(p nx);
   dy = bb/(p ny);
   (*Generate node coordinates, 
   meshnodes should contain (p nx+1)(ny+1)+(p \
ny+1)(nx+1)-(nx+1)(ny+1) nodes*)
   For[hl = 0, hl < ny, hl++, (*loop over ny horizontal lines, 
    the last one is done below, after the loop *)
    For[hln = 0, hln < p nx + 1, 
     hln++, (*loop over p nx+1 nodes on the horizontal line*)
     AppendTo[meshnodes, {hln dx, dy p hl}];
     ];
    For[vl = 1, vl <= p - 1, 
     vl++, (*loop over the p-1 horizontal "lines" that are not part \
of the mesh so we can assign the coordinates to the nodes on the \
vertical lines*)
     For[vln = 0, vln < nx + 1, 
       vln++,(*loop over the nx+1 nodes on the vertical lines*)
       AppendTo[meshnodes, {vln p dx , (p hl + vl) dy}];
       ];
     ];
    ];
   (*Now for the last horizontal line, 
   note that we need to set hl to ny, 
   since we started counting from 0:*)
   hl = ny;
   For[hln = 0, hln < p nx + 1, 
    hln++, (*loop over p nx+1 nodes on the horizontal line*)
    AppendTo[meshnodes, {hln dx, dy p hl}];
    ];
   (*generate the list of cells - each cell is a list of node-
   IDs that are on its border.
   We have nx ny cells, each cell has 4p nodes*)
   meshtopology = Table[{}, nx ny];
   (*label the cells (cx,cy), i.e. (0,0),(1,0,1),...,(nx-1,0),(0,
   1),...,(nx-1,ny-1) etc.*)
   For[cy = 0, cy < ny, cy++,
    For[cx = 0, cx < nx, cx++,
      (*bottom edge*)
      For[i = 0, i < p + 1, i++,
       AppendTo[meshtopology[[cx + cy nx + 1]], 
         i + cx p + cy ((nx + 1) (p - 1) + p nx + 1) + 1];
       ];
      (*right edge*)
      For[i = 0, i < p - 1, i++,
       AppendTo[meshtopology[[cx + cy nx + 1]], 
         cx + 1 + i (nx + 1) + (cy + 1) (p nx + 1) + 
          cy (p - 1) (nx + 1) + 1];
       ];
      (*top edge, from right to left*)
      For[i = p, i >= 0, i--,
       AppendTo[meshtopology[[cx + cy nx + 1]], 
         i + cx p + (cy + 1) ((nx + 1) (p - 1) + p nx + 1) + 1];
       ];
      (*left edge,from top to bottom*)
      For[i = p - 2, i >= 0, i--,
       AppendTo[meshtopology[[cx + cy nx + 1]], 
         cx + i (nx + 1) + (cy + 1) (p nx + 1) + 
          cy (p - 1) (nx + 1) + 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_, p_] := Block[{meshvis, nodevis, v},
   If[order == 1,
    v = {1, 2, 3, 4},
    v = Table[i, {i, 1, 4 p}];
    ];
   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}
   ];

There is really not that much complicated stuff in my solution, I just used some basic for loop indexing wizardry.

When trying to understand my solution, you should note, that I started counting from 0 for all my indices and added 1 at the very end where necessary. This is in part because I am very experienced in C++ (which starts counting from 0 while Mathematica counts from 1), but also because we have many multiplications where it is usefull to have the zeroth row. I recommend making a general sketch of the mesh you want of order $p$ and count the nodes for each horizontal line, for each cell, etc. That way you should arrive at the same equations that I found.

Note, that I changed the ordering of the nodes in each cell. I am going counter-clockwise around the cell, starting at the bottom left corner. I achieved this in part by letting the For-loops for the top and left edge "run backwards".

I am sure that one could work more with Table and similar commands, but since you used For-loops, I did too.

Using the functions like so:

L = 12;
h = 10;
nx = 4;
ny = 5;
order = 6;
{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]

yields the following image

I am also not sure, what you meant with the dimensionality of L and h, but since your examples are all 2-D meshes, I assumed that they are only lengths (and hence have basically no influence on the mesh apart from the distance between the nodes).