#========== Lieu des  points se projetant sur les faces(intérieur) d'un polyèdre ============

restart:
with(plots):

Vdot:=proc(U,V)  #produit scalaire
local u,v; u:=expand(U): v:=expand(V):
       add(u[i]*v[i],i=1..nops(u)):
end:

Vcross:=proc(u,v)# produit vectoriel
   local U,V; U:=expand(u): V:=expand(v):
[U[2]*V[3]-U[3]*V[2],U[3]*V[1]-U[1]*V[3],U[1]*V[2]-U[2]*V[1]]:
end proc:
#---- équation du plan passant par 3 points [x1,y1,z1], [x2,y2,z2],[x3,y3,z3] -------------
EQplan:= proc(x1,y1,z1,x2,y2,z2,x3,y3,z3)
     ((-z2+z3)*y1+(z1-z3)*y2-y3*(z1 - z2))*x+
     ((z2-z3)*x1+(-z1+z3)*x2+x3*(z1-z2))*y+
     ((-y2+y3)*x1+(y1-y3)*x2-x3*(y1-y2))*z+
     (x1*y2-x2*y1)*z3+x2*y3*z1+x3*(y1*z2-y2*z1)-x1*y3*z2:
end:
#====== aretes, sommets d'un polyèdre à partir de la liste des faces ====================
calcul_aretes:= proc(Lfaces)
    convert(map(u ->convert(u,list),{seq(seq({Lfaces[i,j],Lfaces[i,j mod nops(Lfaces[i]) +1]},j=1..nops(Lfaces[i])),i=1..nops(Lfaces))}),set): end:
calcul_sommets:=proc(Lfaces) convert(convert(map(op,map(op,Lfaces)),set),list): end:

 with(geom3d):
#----- choix du polyèdre --------
cuboctahedron(gon,point(O1,0,0,0),sqrt(3)):#dodecahedron(gon,point(O1,0,0,0),sqrt(3)):
Lfaces:=faces(gon): 
Lsom:=vertices(gon):#calcul_sommets(Lfaces)
Laretes:=calcul_aretes(Lfaces):
Laretes:={seq(seq({Lfaces[i,j],Lfaces[i,j mod nops(Lfaces[i]) +1]},j=1..nops(Lfaces[i])),i=1..nops(Lfaces))}:    #nops(%); 

#==== Equation d'un prisme droit sur une face (polygone gauche)
#Equation du plan passant par une arete et  perpendiculaire à la face
EQplan_ar:=proc(A,B,N)   #---- N est un vecteur normal à la face -----
   simplify(EQplan(op(A),op(B),op(expand((A+B)/2+N) )));
end:
EQ_prismeFace:=proc(face) local i,j,LEQ_ar,N,n, G , LEQ_sign:
    n:=nops(face); 
    N:=Vcross(face[2]-face[1],face[2]-face[3]):
    LEQ_ar:=[seq( EQplan_ar(face[i],face[(i mod n) +1],N ),i=1..n)]:
    G:=convert(face,`+`)/n:
    LEQ_sign:=map(u -> sign(evalf(subs(x=G[1],y=G[2],z=G[3],u)))*u,LEQ_ar):
 end:
LEQ_prism_gon:=(Lfaces)->  map(EQ_prismeFace,Lfaces): #--- liste des équations pour les faces du polyèdre ------


EQ1:=op(evalf( map(u ->min(op(u)),LEQ_prism_gon(Lfaces)),4)):
EQ:='min'(EQ1);  #` min( op( map(u ->min(op(u)),LEQ_prism_gon(Lfaces))))`;

#===== le tracé du lieu =======
ro:=0.03:  K:=1: np:=120:
display([#polygonplot3d(Lfaces,style=line,thickness=5,color=gold),
         seq(tubeplot(expand(Laretes[i,1]+t*(Laretes[i,2]-Laretes[i,1])),radius=ro,t=0..1,color=gold),i=1..nops(Laretes)),
         implicitplot3d(EQ,x=-K..K,y=-K..K,z=-K..K,grid=[np,np,np], color=cyan,transparency=0.)
         ],style=patchnogrid,lightmodel=light3,scaling=constrained,axes=none,orientation=[-56,80,5]);