//phi=(1+sqrt(5))/2;
//
//points=[
//   [0,1,phi], [0,1,-phi], [0,-1,phi], [0,-1,-phi],
//   [1,phi,0], [1,-phi,0], [-1,phi,0], [-1,-phi,0],
//   [phi,0,1], [phi,0,-1], [-phi,0,1], [-phi,0,-1]
//];
//
//triangles=[
//   [0,4,8],[0,8,2],[0,2,10],[0,10,6],[0,6,4],
//   [3,11,7],[3,7,5],[3,5,9],[3,9,1],[3,1,11],
//   [7,11,10],[11,1,6],[1,9,4],[9,5,8],[5,7,2],
//   [2,8,5],[8,4,9],[4,6,1],[6,10,11],[10,2,7]
//];

//polyhedron(points=points,faces=triangles);

function z(x,y)=(x/4+y/4+0.75)*sin(180*x)*sin(180*y);

//function z(x,y)=sqrt(1+x*x-y*y);
steps=7;

ps=concat(
   [for(x=[-1:2/(steps-1):1]) for(y=[-1*2/sqrt(3)+x/2:(4/sqrt(3))/(steps-1):1*2/sqrt(3)+x/2]) [x,y,z(x,y)]]
       ,
   [for(x=[-1:2/(steps-1):1]) for(y=[-1*2/sqrt(3)+x/2:(4/sqrt(3))/(steps-1):1*2/sqrt(3)+x/2]) [x,y,z(x,y)-0.5]]);
ts=concat(
   flatten([for(x=[0:steps-2]) for(y=[0:steps-2]) [[x+y*steps,(x+1)+y*steps,x+(y+1)*steps],[(x+1)+y*steps,(x+1)+(y+1)*steps,x+(y+1)*steps]]]),
   flatten([for(x=[0:steps-2]) for(y=[0:steps-2]) [
       [steps*steps+(x+1)+y*steps,steps*steps+x+y*steps,steps*steps+x+(y+1)*steps],
       [steps*steps+(x+1)+(y+1)*steps,steps*steps+(x+1)+y*steps,steps*steps+x+(y+1)*steps]]]),
   [for(x=[0:steps-2]) [x,x+steps*steps,x+1+steps*steps]],       
   [for(x=[0:steps-2]) [x,x+1+steps*steps,x+1]],
       
   [for(x=[0:steps-2]) [x+steps*steps+steps*(steps-1),x+steps*(steps-1),x+1+steps*steps+steps*(steps-1)]],       
   [for(x=[0:steps-2]) [x+1+steps*steps+steps*(steps-1),x+steps*(steps-1),x+1+steps*(steps-1)]],
   [for(y=[0:steps-2]) [y*steps,(y+1)*steps,y*steps+steps*steps]],
   [for(y=[0:steps-2]) [(y+1)*steps,(y+1)*steps+steps*steps,y*steps+steps*steps]],
   [for(y=[0:steps-2]) [(y+1)*steps+steps-1,y*steps+steps-1,y*steps+steps*steps+steps-1]],
   [for(y=[0:steps-2]) [(y+1)*steps+steps*steps+steps-1,(y+1)*steps+steps-1,y*steps+steps*steps+steps-1]]
);

echo(len(ps));
//echo(ts);
polyhedron(points=ps,faces=ts);

points=ps;
triangles=ts;


for(facet=triangles) {
    assert(3==len(facet),"Mesh facets must be triangles");
}

function flatten(list) = [for(a=list) for (b=a) b];

function tquicksort(arr) = !(len(arr)>0) ? [] : let(
    pivot   = arr[floor(len(arr)/2)],
    lesser  = [ for (y = arr) if (y[0]  < pivot[0] || (y[0]==pivot[0] && y[1]<pivot[1])) y ],
    equal   = [ for (y = arr) if (y[0] == pivot[0] && y[1]==pivot[1]) y ],
    greater = [ for (y = arr) if (y[0]  > pivot[0] || (y[0]==pivot[0] && y[1]>pivot[1])) y ]
) concat(
    tquicksort(lesser), equal, tquicksort(greater)
);

function tbsearch(pair,list,mymin=0,mymax=-1)=
    let (
       realmax=(mymax<mymin)?len(list)-1:mymax,
       test=floor((realmax+mymin)/2)
    )
    (realmax==mymin)?
       realmax
    :((list[test][0]<pair[0])||(list[test][0]==pair[0] && list[test][1]<pair[1]))?
 //[pair,test+1,realmax]
       tbsearch(pair,list,test+1,realmax)
    :
 //[pair,mymin,test]
       tbsearch(pair,list,mymin,test)
;
    
function find_point(pair,list)=
    (pair[0]>pair[1])?
       tbsearch([pair[1],pair[0]],list)
    :
       tbsearch(pair,list)
;
    
//Butterfly interpolation is:
//(1/2)*(p1+p2)+2w*(p3+p4) -w(p5+p6+p7+p8) 
//alpha (p1+p2) +beta (p3+p4) + gamma (p5+p6+p7+p8)
//Referencing
//"A Butterfly Subdivision Scheme for Surface Interpolation with
//Tension Control"
//By Nira Dyn, David Levin, and John Gregory
//1/16 is probably the max value we want for omega
omega=1/16;
alpha=1/2;
beta=+2*omega;
gamma=-omega;

function interpolated_point(segment,full_triangle_list,points,tension=0)=
   let(
      i1=segment[0],
      i2=segment[1],
      i3=full_triangle_list[tbsearch([i1,i2],full_triangle_list)][2],
      i4=full_triangle_list[tbsearch([i2,i1],full_triangle_list)][2],
      i5=full_triangle_list[tbsearch([i3,i2],full_triangle_list)][2],
      i6=full_triangle_list[tbsearch([i1,i3],full_triangle_list)][2],
      i7=full_triangle_list[tbsearch([i2,i4],full_triangle_list)][2],
      i8=full_triangle_list[tbsearch([i4,i1],full_triangle_list)][2],
      p1=points[i1],
      p2=points[i2],
      p3=points[i3],
      p4=points[i4],
      p5=points[i5],
      p6=points[i6],
      p7=points[i7],
      p8=points[i8]
   )
   (p1+p2)*alpha+(p3+p4)*beta+(p5+p6+p7+p8)*gamma
;
    
//full_triangle_list=tquicksort(flatten([for(facet=triangles) [[facet[0],facet[1],facet[2]],[facet[1],facet[2],facet[0]],[facet[2],facet[0],facet[1]]]]));
//tension=0;

//n=[1,2,0];

//inpoints=tquicksort([for(t=triangles) for(i=[0,1,2]) if(t[i]<t[n[i]]) [t[i],t[n[i]],interpolated_point([t[i],t[n[i]]],full_triangle_list,points,tension)]]);

//new_points=concat(points, [for (s=inpoints) s[2]]);

//
//new_triangles=flatten([for(t=triangles) let(
//       a=len(points)+find_point([t[0],t[1]],inpoints),
//       b=len(points)+find_point([t[1],t[2]],inpoints),
//       c=len(points)+find_point([t[2],t[0]],inpoints)
//      ) [[t[0],a,c],[a,t[1],b],[b,t[2],c],[a,b,c]]]);

//translate([5,0,0]) { polyhedron(points=new_points,faces=new_triangles);
//}

function interpolate_mesh(points,triangles)=
   let(
       ftm=tquicksort(flatten([for(facet=triangles) [[facet[0],facet[1],facet[2]],[facet[1],facet[2],facet[0]],[facet[2],facet[0],facet[1]]]])),
       n=[1,2,0],
       tension=0,
       ip=tquicksort([for(t=triangles) for(i=[0,1,2]) if(t[i]<t[n[i]]) [t[i],t[n[i]],interpolated_point([t[i],t[n[i]]],ftm,points,tension)]]),
       np=concat(points, [for (s=ip) s[2]]),
       nt=flatten([for(t=triangles) let(
       a=len(points)+find_point([t[0],t[1]],ip),
       b=len(points)+find_point([t[1],t[2]],ip),
       c=len(points)+find_point([t[2],t[0]],ip)
      ) [[t[0],a,c],[a,t[1],b],[b,t[2],c],[a,b,c]]])
    )
    [np,nt]
 ;
 
new=interpolate_mesh(points,triangles);     
translate([3,0,0]) { polyhedron(points=new[0],faces=new[1]);
}

newnew=interpolate_mesh(new[0],new[1]);     
translate([6,0,0]) { polyhedron(points=newnew[0],faces=newnew[1]);
}

newnewnew=interpolate_mesh(newnew[0],newnew[1]);     
translate([9,0,0]) { polyhedron(points=newnewnew[0],faces=newnewnew[1]);
}
////
//newnewnewnew=interpolate_mesh(newnewnew[0],newnewnew[1]);     
//translate([12,0,0]) { polyhedron(points=newnewnewnew[0],faces=newnewnewnew[1]);
//}

ps2=concat(
   [for(x=[-1:2/((steps*8)-1):1]) for(y=[-1*2/sqrt(3)+x/2:(4/sqrt(3))/((steps*8)-1):1*2/sqrt(3)+x/2]) [x,y,z(x,y)]]
       ,
   [for(x=[-1:2/((steps*8)-1):1]) for(y=[-1*2/sqrt(3)+x/2:(4/sqrt(3))/((steps*8)-1):1*2/sqrt(3)+x/2]) [x,y,z(x,y)-0.5]]);
ts2=concat(
   flatten([for(x=[0:(steps*8)-2]) for(y=[0:(steps*8)-2]) [[x+y*(steps*8),(x+1)+y*(steps*8),x+(y+1)*(steps*8)],[(x+1)+y*(steps*8),(x+1)+(y+1)*(steps*8),x+(y+1)*(steps*8)]]]),
   flatten([for(x=[0:(steps*8)-2]) for(y=[0:(steps*8)-2]) [
       [(steps*8)*(steps*8)+(x+1)+y*(steps*8),(steps*8)*(steps*8)+x+y*(steps*8),(steps*8)*(steps*8)+x+(y+1)*(steps*8)],
       [(steps*8)*(steps*8)+(x+1)+(y+1)*(steps*8),(steps*8)*(steps*8)+(x+1)+y*(steps*8),(steps*8)*(steps*8)+x+(y+1)*(steps*8)]]]),
   [for(x=[0:(steps*8)-2]) [x,x+(steps*8)*(steps*8),x+1+(steps*8)*(steps*8)]],       
   [for(x=[0:(steps*8)-2]) [x,x+1+(steps*8)*(steps*8),x+1]],
       
   [for(x=[0:(steps*8)-2]) [x+(steps*8)*(steps*8)+(steps*8)*((steps*8)-1),x+(steps*8)*((steps*8)-1),x+1+(steps*8)*(steps*8)+(steps*8)*((steps*8)-1)]],       
   [for(x=[0:(steps*8)-2]) [x+1+(steps*8)*(steps*8)+(steps*8)*((steps*8)-1),x+(steps*8)*((steps*8)-1),x+1+(steps*8)*((steps*8)-1)]],
   [for(y=[0:(steps*8)-2]) [y*(steps*8),(y+1)*(steps*8),y*(steps*8)+(steps*8)*(steps*8)]],
   [for(y=[0:(steps*8)-2]) [(y+1)*(steps*8),(y+1)*(steps*8)+(steps*8)*(steps*8),y*(steps*8)+(steps*8)*(steps*8)]],
   [for(y=[0:(steps*8)-2]) [(y+1)*(steps*8)+(steps*8)-1,y*(steps*8)+(steps*8)-1,y*(steps*8)+(steps*8)*(steps*8)+(steps*8)-1]],
   [for(y=[0:(steps*8)-2]) [(y+1)*(steps*8)+(steps*8)*(steps*8)+(steps*8)-1,(y+1)*(steps*8)+(steps*8)-1,y*(steps*8)+(steps*8)*(steps*8)+(steps*8)-1]]
);

color("pink") translate([0,3,0]) { polyhedron(points=ps2,faces=ts2);
}

color("pink") translate([9,3,0]) { polyhedron(points=ps2,faces=ts2);
}

for(p=ps) {
    color("red") translate(p) sphere(r=0.025);
    color("red") translate(p+[9,0,0]) sphere(r=0.025);
    color("red") translate(p+[9,3,0]) sphere(r=0.025);    
}