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);

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:
//alpha*(p1+p2)/2 +beta (p3+p4) + gamma (p5+p6+p7+p8)
alpha=1/2+6/128;
beta=-4/128;
gamma=-1/256;
//...tension is something involving those values, but the paper I
//had as a reference seems to no longer be on-line.
//so I'm just leaving them as free variables.

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([i4,i2],full_triangle_list)][2],
      i8=full_triangle_list[tbsearch([i1,i4],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([5,0,0]) { polyhedron(points=new[0],faces=new[1]);
}

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

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

newnewnewnew=interpolate_mesh(newnewnew[0],newnewnew[1]);     
translate([20,0,0]) { polyhedron(points=newnewnewnew[0],faces=newnewnewnew[1]);
}

newnewnewnewnew=interpolate_mesh(newnewnewnew[0],newnewnewnew[1]);     
translate([25,0,0]) { polyhedron(points=newnewnewnewnew[0],faces=newnewnewnewnew[1]);
}

echo(newnewnewnewnew[0][0]);