SetOptions[$FrontEndSession,AutoStyleOptions->{"HighlightFormattingErrors"->False}];
ClassicLights = {{"Ambient",RGBColor[0.312,0.188,0.4]},{"Directional",RGBColor[0.8,0,0],ImageScaled[{2,0,2}]},{"Directional",RGBColor[0,0.8,0],ImageScaled[{2,2,2}]},{"Directional",RGBColor[0,0,0.8],ImageScaled[{0,2,2}]}}

M={{1,Cos[Pi/3]},{0,Sin[Pi/3]}};
Mi=Inverse[M];

(*map a point in the triangle T to the corresponding point on the standard 2-simplex S*)
TtoS[xy_]:=Module[{p1,p2},{p1,p2}=Mi.xy;{p1,p2,1-p1-p2}];

(*calculate the probability density as a function on T*)
f[xy_,a_]:=Module[{a1,a2,a3,x,y,z,B},{a1,a2,a3}=a;
{x,y,z}=TtoS[xy];
B=Gamma[a1+a2+a3]/(Gamma[a1]*Gamma[a2]*Gamma[a3]);
B*(x^(a1-1))*(y^(a2-1))*(z^(a3-1))];

(*recursively tesselate triangle (x1,x2,x3) by bisecting its edges and return a list of triangles as in input suitable for Graphics3D*)
Tesselate[x1_,x2_,x3_,alpha_]:=Module[{d1=(x2-x1),d2=(x3-x2),d3=(x1-x3),m,g},m=Max[{d1.d1,d2.d2,d3.d3}];
If[m<0.0002,(*stop when the length of every edge is below threshold*){{EdgeForm[],(*don't display edges of triangles*)Polygon[{{x1[[1]],x1[[2]],f[x1,alpha]},{x2[[1]],x2[[2]],f[x2,alpha]},{x3[[1]],x3[[2]],f[x3,alpha]}}]}},Module[{m1=.5*(x1+x2),m2=.5*(x2+x3),m3=.5*(x3+x1)},Join[Tesselate[x1,m1,m3,alpha],Tesselate[m1,x2,m2,alpha],Tesselate[m3,m1,m2,alpha],Tesselate[m3,m2,x3,alpha]]]]]

Dirichlet[alpha_, zmax_:6] :=
Show[Graphics3D[Tesselate[{0,0},{1,0},{Cos[Pi/3],Sin[Pi/3]},alpha],
Lighting->ClassicLights
],
Graphics3D[triangle],
Boxed->False,PlotRange->{{0,1},{0,1},{0,zmax}},
Axes->True, AxesEdge->{None,None,{1,-1}},
BoxRatios->{1,1,1/2},
ViewPoint->{-1,-2,1.7}]

Dirichlet[{2,4,2},10]

Dirichlet[{0.1,0.1,0.1},20]

Dirichlet[{1,1,1},4]