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.

## Geodesic Fractal

```import rhinoscriptsyntax as rs
import Rhino.Geometry as rg
import random
import math
STOP=5

(_residual,_jacob,S1,S2,S3,filename)=input

for i in range(STOP+2):

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)
if(_cycle==11):
return []
elif _cycle==2 and deepness<STOP:
f.append(rg.Polyline())
f.append(rg.Polyline())
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:
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