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.