Jim Propp once asked: "Is there a way to relax an approximation to a space-filling curve in continuous time so that it works out its kinks and regresses to simpler approximations? No interim self-intersections please!" Julian Ziegler Hunts developed a Fourier matrix product that can produce the analogous animation.

Julian wrote a nifty Fourier expander for recursive Koch polygons. E.g.,

ptsgnlst2Fouriermat[{0, 1, I^(2/3), 1 + I^(2/3), 2}, {1, -1, 1, -1}]

spacefills the triangle joining 0 to 2 via 1+i?3. Actually, it makes an infinite 3x3 matrix product for the coefficients a(k). Then Sum a(k+1/m) exp(2 i ? (k+1/m)) repeats the fractal on the sides of an m-gon. The GIF just accumulates consecutive harmonics. These are the definitions of functions:

ClearAll[ptsgnlst2Fouriermat]; 
ptsgnlst2Fouriermat[points_List, signs_List] := 
 Block[{Amat, Bmat, M, 
   klists = {(3 - signs)/2, Reverse[(signs + 3)/2]}, 
   pt0 = points[[1]], pt1 = points[[-1]], n = Length[signs], m = 2, c,x}, 
Amat = {#, Reverse[#]} &[
    signs*(Rest[points] - Most[points])/(pt1 - pt0)]; 
  Bmat = {#, Reverse[#]} &[-signs*pt0/(pt1 - pt0) + 
     points[[Range[n] + (1 - signs)/2]]]; 
  Function[x, c] /. c -> Append[Table[Append[
       Table[Exp[-I*x*(i - 1)/n]*Amat[[j, i]], {i, 1, n}].IdentityMatrix[m][[klists[[j]]]]/n, 
       Sum[Exp[-I*x*(i - 1)/n]*Bmat[[j, i]], {i, 1, n}]*I*(Exp[-I*x/n] - 1)/x], {j, 1, m}], 
     Table[Boole[i == m], {i, 0, m}]]];

    ClearAll[mat2series]; 
    mat2series[var_, mat_, n_, m_: 1, terms_, prodterms_, prec_: 60] := 
     Block[{A, x}, A[0] = Limit[mat[x], x -> 0]; A[x_] = mat[x]; 
      Exp[2*I*\[Pi]*Range[-terms + 1/m, terms + 1/m]*
         var].Table[(Dot @@ 
           Table[A[N[2*\[Pi]*(k + 1/m)/n^r, prec]], {r, 0, prodterms}])[[1 ;; -2, -1]], {k, -terms, terms}]]

Here is an example of 3-fold rotational symmetry:

Block[{start = 0, end = 1, m = -3., 
  f = mat2series[t, %, 7, -3., 33, MachinePrecision][[1]]}, 
 ParametricPlot[ReIm[f + Sum[
     E^(2 \[Pi] (k + 1/m) I t) I*(start - Exp[-2*I*\[Pi]/m]*end)/
         2/\[Pi]/(k + 1/m), {k, -33, 33}]], {t, 0, m}, Axes -> False]]

And here it is with higher details:

Block[{start = 0, end = 1, m = -3., 
  f = mat2series[t, %%, 7, -3., 69, MachinePrecision][[1]]}, 
 ParametricPlot[ReIm[f + Sum[
     E^(2 \[Pi] (k + 1/m) I t) I*(start - Exp[-2*I*\[Pi]/m]*end)/
         2/\[Pi]/(k + 1/m), {k, -69, 69}]], {t, 0, m}, Axes -> False]]

And given enough time the code below will produce the animation at the top of this post:

Table[Block[{start = 0, end = 2, m = 6., 
    f = mat2series[t, ptsgnlst2Fouriermat[{0, 1, I^(2/3), 1 + I^(2/3), 2}, 
    {1, -1, 1, -1}], 4, m, n, MachinePrecision][[1]]}, 
   ParametricPlot[ReIm[f + 
      Sum[E^(2 \[Pi] (k + 1/m) I t) I*(start - Exp[-2*I*\[Pi]/m]*end)/
           2/\[Pi]/(k + 1/m), {k, -n, n}]], {t, 0, 6}, Axes -> False]], {n, 99}];

Export["/Users/billgosper/Movies/hellodoily.gif", 
 Join[ConstantArray[%[[1]], 6], %, ConstantArray[%[[-1]], 9]]]