Masaaki MIKI

import rhinoscriptsyntax as rs
import Rhino.Geometry as rg
import random
import math
STOP=5
Radius=list()
 
(_residual,_jacob,S1,S2,S3,filename)=input
 
for i in range(STOP+2):
    Radius.append(math.pow(0.8,i)*0.015)
 
def residual(_x):
    global _residual
    global _jacob
    r=_residual(_x)
    return r
 
def jacob(_x):
    global _residual
    global _jacob
    j=_jacob(_x)
    if j.Length<0.001:
        dd=list()
        S=0.00001
        j=rg.Vector3d(0,0,0)
        dd.append(_x+rg.Vector3d(S,S,0))
        dd.append(_x+rg.Vector3d(-S,S,0))
        dd.append(_x+rg.Vector3d(-S,-S,0))
        dd.append(_x+rg.Vector3d(-S,S,0))
        dd.append(_x+rg.Vector3d(0,S,S))
        dd.append(_x+rg.Vector3d(0,-S,S))
        dd.append(_x+rg.Vector3d(0,-S,-S))
        dd.append(_x+rg.Vector3d(0,-S,S))
        for d in dd:
            j+=_jacob(d)/8
    return j
 
def pinv(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    N=x**2+y**2+z**2
    return rg.Vector3d(x/N,y/N,z/N)
 
def varphi(_x,_X):
    J=jacob(_X)
    iJ=pinv(J)
    N=_x.Length
    lam=_x*J
    _x=_x-lam*iJ
    _x.Unitize()
    _x*=N
    return _x
 
def psi(_x):
    for i in range(100):
        J=jacob(_x)
        iJ=pinv(J)
        Resid=residual(_x)
        dx=iJ*Resid
        _x-=dx
        if(math.fabs(Resid)<0.0001):
            break
    return _x
 
def tic(job):
    (_X,_V,_R,_f,_cycle,deepness)=job #unpacking
    global a
    global f
    global dd
    for t in range(20):
        _V=varphi(_V,_X)
        if(t==0 and _R!=0):
            N=jacob(_X)
            N.Unitize()
            _V*=0.8
            D=rg.Vector3d.CrossProduct(N,_V)
            D.Unitize()
            D*=_V.Length
            _V=_V*math.cos(_R/360*2*math.pi)+D*math.sin(_R/360*2*math.pi)
        _X+=_V*dt
        _X=psi(_X)
        _f.Add(_X)
    if(_cycle==11):
        #createPipe(_f,Radius[deepness])
        dd.append((_f,Radius[deepness]))
        return []
    elif _cycle==2 and deepness<STOP:
        f.append(rg.Polyline())
        f[len(f)-1].Add(_X)
        f.append(rg.Polyline())
        f[len(f)-1].Add(_X)
        return [(_X,_V,0,_f,_cycle+1,deepness),
        (_X,_V,30,f[len(f)-2],0,deepness+1),
        (_X,_V,-30,f[len(f)-1],0,deepness+1)]
    else:
        return [(_X,_V,0,_f,_cycle+1,deepness)]
dt=0.02
if x:
    for s in range(10):
        if(jobNum<len(jobList)):
            ret=tic(jobList[jobNum])
            jobList.extend(ret)
            jobNum=jobNum+1
    a=list()
    for l in f:
        a.append(rs.AddPolyline(l))
    if(jobNum==len(jobList)):
        if output==False:
            file=open(filename,'w')
            for d in dd:
                (p,R)=d
                file.write('P'+' , ' +str(len(p))+'\n')
                file.write('R , '+str(R)+'\n')
                for s in p:
                    file.write('C , '+str(s.X)+' , '+str(s.Y)+' , '+str(s.Z)+"\n")
            file.close()
            output=True
        print "Finished!"
else:
    output=False
    cycle=0;
    jobNum=0
    jobList=list()
    f=list()
    g=list()
    dd=list()
    f.append(rg.Polyline())
    f.append(rg.Polyline())
    f.append(rg.Polyline())
    jobList.append((S1[0],S1[1],0,f[0],cycle,0))
    jobList.append((S2[0],S2[1],0,f[1],cycle,0))
    jobList.append((S3[0],S3[1],0,f[2],cycle,0))

import Rhino.Geometry as rg
import math
R=0.9/2
r=0.6/2
def residual(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    L=math.sqrt(x**2+z**2)
    return ((x-R*x/L)**2)+((z-R*z/L)**2)+(y**2)-r**2
def jacob(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    L=math.sqrt(x**2+z**2)
    fx=2*(x-R*x/L)*(1-R/L+R*x/(L**3)*x)+2*(z-R*z/L)*(R*z/(L**3)*x)
    fz=2*(z-R*z/L)*(1-R/L+R*z/(L**3)*z)+2*(x-R*x/L)*(R*x/(L**3)*z)
    fy=2*y
    return rg.Vector3d(fx,fy,fz)
X1=rg.Point3d(0.2,0.4,-0.76)
V1=rg.Vector3d(0.12,0.24,0.32)*0.75
X2=rg.Point3d(-0.35,0.4,-0.84)
V2=rg.Vector3d(-0.18,0.24,0.3)*0.75
X3=rg.Point3d(-0.11,-0.4,-0.76)
V3=rg.Vector3d(-0.12,-0.4,0.4)*0.75
 
a=((residual,jacob,(X1,V1),(X2,V2),(X3,V3),"c:/out/donut.txt"),)

import Rhino.Geometry as rg
S=1.3
def residual(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    return (2*(S*x**2)+2*(S*y**2)+(S*z**2)-1)**3-0.1*(S*x**2)*(S*z**3)-(S*y**2)*(S*z**3)
def jacob(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    fx=3*((2*(S*x**2)+2*(S*y**2)+(S*z**2)-1)**2)*4*S*x-0.2*S*x*(S*z**3)
    fy=3*((2*(S*x**2)+2*(S*y**2)+(S*z**2)-1)**2)*4*S*y-2*S*y*(S*z**3)
    fz=3*((2*(S*x**2)+2*(S*y**2)+(S*z**2)-1)**2)*2*S*z-0.3*(S*x**2)*(S*z**2)-3*(S*y**2)*(S*z**2)
    return rg.Vector3d(fx,fy,fz)
X1=rg.Point3d(0.2,0.4,-0.68)
V1=rg.Vector3d(0.18,0.14,0.5)*0.75
X2=rg.Point3d(-0.2,0,-0.4)
V2=rg.Vector3d(-0.28,0.14,0.5)*0.75
X3=rg.Point3d(0.2,-0.4,-0.69)
V3=rg.Vector3d(-0.06,-0.2,0.5)*0.75
 
a=((residual,jacob,(X1,V1),(X2,V2),(X3,V3),"c:/out/heart.txt"),)

import Rhino.Geometry as rg
x0=1
y0=0
z0=0
def residual(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    return ((2*x-2*x0)**4+(2*y-2*y0)**4+(2*z-2*z0)**4-((2*x-2*x0)**2+(2*y-2*y0)**2+(2*z-2*z0)**2))
def jacob(_x):
    x=_x.X
    y=_x.Y
    z=_x.Z
    fx=4*((2*x-2*x0)**3)*2-2*2*(2*x-2*x0)
    fy=4*((2*y-2*y0)**3)*2-2*2*(2*y-2*y0)
    fz=4*((2*z-2*z0)**3)*2-2*2*(2*z-2*z0)
    return rg.Vector3d(fx,fy,fz)
X1=rg.Point3d(1.1,0.4,-0.55)
V1=rg.Vector3d(-0.07,0.3,0.6)*0.55
X2=rg.Point3d(1.6,0.25,-0.6)
V2=rg.Vector3d(0.06,-0.12,0.6)*0.6
X3=rg.Point3d(1.6,-0.4,-0.6)
V3=rg.Vector3d(0.2,-0.1,0.6)*0.6
 
a=((residual,jacob,(X1,V1),(X2,V2),(X3,V3),"c:/out/cube.txt"),)