Geodesicfractal
Here is an example of geometric application of geodesic dynamic relaxation method presented in:
M. Miki, S. Adriaenssens, T. Igarashi, and K. Kawaguchi, The geodesic dynamic relaxation method for problems of equilibrium with equality constraint conditions, Int. J. Numer. Meth. Engng,99, pages 682–710. doi: 10.1002/nme.4713.
How it works
Settings (Grasshopper + Rhinoceros)

Geodesic Fractal
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))
Donut
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"),)
Heart
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"),)
Cube
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"),)
page revision: 6, last edited: 01 May 2015 07:45