Unprotect[Power]; 0. ^0 = 1; 0^0 = 1;

a = {0, 0}; b = {2, 0}; c = {0, 2};
phi1 = 
      If[b[[1]] == a[[1]], Pi/2, 
        ArcTan[(b[[2]] - a[[2]])/(b[[1]] - 
                a[[1]])]];  
    phi2 = Pi/2 - 
        If[c[[2]] == a[[2]], Pi/2, 
          ArcTan[(c[[1]] - a[[1]])/(c[[2]] - 
                  a[[2]])]]; 
    alpha = Pi/2 - 
        If[b[[1]] == c[[1]], Pi/2, 
          ArcTan[(b[[2]] - c[[2]])/(b[[1]] - 
                  c[[1]])]] - 
        If[b[[2]] == a[[2]], Pi/2, 
          ArcTan[(b[[1]] - a[[1]])/(b[[2]] - 
                  a[[2]])]]; 
    ab = Sqrt[(b[[1]] - a[[1]])^2 + (b[[2]] - 
                a[[2]])^2];

area[x_, y_, z_] := 
      Abs[Det[{{x[[1]], y[[1]], z[[1]]}, {x[[2]], y[[2]], 
                z[[2]]}, {1, 1, 1}}]/2]
barycentric[a_, b_, c_][
        p_] := {area[p, b, c], area[a, p, c], area[a, b, p]}/
        area[a, b, c]
euklidic[a_, b_, c_][{x_, y_, z_}] := 
      x a + y b + z c

Bernstein[j_][u_] := 
        Multinomial[Sequence @@ j]  u[[1]]^j[[1]] u[[2]]^
            j[[2]] u[[3]]^j[[3]];

bb[f_, a_, b_, c_][j_] := 
      Block[{n = Plus @@ j, u = j[[1]]/n, v = j[[2]]/n}, 
        f[euklidic[a, b, c][{u, v, 1 - u - v}]]]

casteljau[f_, n_, a_, b_, c_][u_] := 
      Block[{j = {0, 0, 0}, val = 0}, 
        For[,j[[1]] ≤ n,j[[1]]++,For[
            j[[2]] = 0, 
            j[[2]] ≤ 
              n - j[[1]], j[[2]]++, 
{j[[3]] = n - j[[2]] - j[[1]], 
              val += bb[f, a, b, c][j] Bernstein[j][u]}]]; 
        val]

f[x_, 
        y_] := Bernstein[{2, 2, 0}][barycentric[a, b, c][{x, 
            y}]]

f[{x_, y_}] := f[x, y]

q[f_, n_, a_, b_, c_][phi_, r_] := 
      Block[{s = Sin[alpha] ab/Sin[alpha + phi - phi1], 
          x = a[[1]] + r s Cos[phi], 
          y = a[[2]] + r s Sin[phi]}, {x, 
          y, casteljau[f, n, a, b, 
              c][barycentric[a, b, c][{x, y}]]}]

ParametricPlot3D[q[f, 10, a, b, c][phi, r], {r, 0, 1}, {phi, 
          phi1, phi2}, PlotPoints -> 15, PlotRange -> All, 
        Compiled -> False];

f[x_, 
          y_] := Bernstein[{3, 0, 0}][barycentric[a, b, c][{x, 
              y}]];

ParametricPlot3D[q[f, 10, a, b, c][phi, r], {r, 0, 1}, {phi, 
          phi1, phi2}, PlotPoints -> 15, PlotRange -> All, 
        Compiled -> False];

f[x_, y_] := Sin[8  x y]

ParametricPlot3D[q[f, 10, a, b, c][phi, r], {r, 0, 1}, {phi, 
          phi1, phi2}, PlotPoints -> 20, PlotRange -> All, 
        Compiled -> False];

For[n = 1, n ≤ 20, n++, 
        ParametricPlot3D[q[f, n][phi, r], {r, 0, 1}, {phi, phi1, phi2}, 
          PlotPoints -> 30, PlotRange -> All, 
          Compiled -> False, Boxed -> False, 
          Axes -> False]];

r[f_, n_][x_, y_] := 
      Block[{
left = Floor[x], 
          down = Floor[y], 
a = 
            If[y - down ≤ 1 - (x - left), {left, 
                down}, {left + 1, down + 1}], 
b = {left + 
                1, down}, 
c = {left, 
              down + 1}
}, casteljau[f, n, a, b, 
            c][barycentric[a, b, c][{x, y}]]]

f[x_, y_] := Sin[x y]

ParametricPlot3D[{x, y, r[f, 13][x, y]}, {x, -3, 
          3}, {y, -3, 3}, PlotPoints -> 40];

Plot3D[f[x, y], {x, -3, 3}, {y, -3, 3}, 
        PlotPoints -> 40];