-
Notifications
You must be signed in to change notification settings - Fork 0
/
GrowingShoot.py
535 lines (489 loc) · 23.9 KB
/
GrowingShoot.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
'''This code implements a 3D rod model for growing plant shoots. Growth is localized at the organ tip, in a growing region of constant length.
This code includes lignification (by rod stiffening), plant responses to gravity (sensed by means of statoliths), to bending (proprioception)
and an endogenous oscillator.'''
from __future__ import print_function
from __future__ import unicode_literals
from fenics import *
import time
import numpy as np
import multiprocessing as mp
import matplotlib as mpl
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
####################################### PLOT SETTINGS AND FUNCTIONS #########################################
#mpl.rcParams['text.usetex']=True # For displaying in LaTex font
height = 15 # figure height
width = 15 # figure width
ArrowL = 0.01 # lenght of arrows for directors and statoliths direction
MarkerSize = 0.65 # Marker size for tip projections
TickSize = 20
LabelSize = 24
TitleSize = 28
LabelSpace= 28
TickSpace = 12
TipPts = 12000 # Number of displayed points for tip projections
RodLineWidth = 3 # Line width of rod
OldRodColor = '#8d5524'
RodColor = '#4f8f00'
xr = 0.05 # right axis limit
xl = -xr # left axis limit
Background = 'w' # facecolor
def save_plot():
# Plot
fig = plt.figure(figsize=(height,width))
ax = fig.add_subplot(111, projection='3d',facecolor=Background)
plt.title('Shoot length: %.5f cm, Shoot time: %.4f min' % (Lt*100,t*Tau_s/60),fontsize=TitleSize)
ax.quiver(x[ArrowIndex],y[ArrowIndex],z[ArrowIndex],D1[2],D1[0],D1[1], length=ArrowL, normalize=True, color='b')
ax.quiver(x[ArrowIndex],y[ArrowIndex],z[ArrowIndex],D2[2],D2[0],D2[1], length=ArrowL, normalize=True)
ax.quiver(x[ArrowIndex],y[ArrowIndex],z[ArrowIndex],D3[2],D3[0],D3[1], length=ArrowL, normalize=True, color='r')
ax.quiver(x[ArrowIndex],y[ArrowIndex],z[ArrowIndex],H[2], H[0], H[1], length=ArrowL, normalize=True, color='k')
ax.auto_scale_xyz([xl,xr], [xl,xr], [0, (xr-xl)])
X1, X2 = ax.get_xlim3d()
Y1, Y2 = ax.get_ylim3d()
Z1, Z2 = ax.get_zlim3d()
ax.plot(Xtip[-TipPts:],Ytip[-TipPts:],[Z1]*np.size(Ztip[-TipPts:]),'.',markersize=MarkerSize)
ax.plot(Xtip[-TipPts:],[Y2]*np.size(Ytip[-TipPts:]),Ztip[-TipPts:],'.',markersize=MarkerSize)
ax.plot([X1]*np.size(Xtip[-TipPts:]),Ytip[-TipPts:],Ztip[-TipPts:],'.',markersize=MarkerSize)
ax.plot(x[young],y[young],z[young],color=RodColor,linewidth=RodLineWidth)
ax.plot(x[old],y[old],z[old],color=OldRodColor,linewidth=RodLineWidth)
# Set labels
ax.set_xlabel("$\mathbf{e}_3$ [m]",fontsize=LabelSize, labelpad=LabelSpace)
ax.set_ylabel("$\mathbf{e}_1$ [m]",fontsize=LabelSize, labelpad=LabelSpace)
ax.set_zlabel("$\mathbf{e}_2$ [m]",fontsize=LabelSize, labelpad=LabelSpace)
# Set ticks
ax.xaxis.set_tick_params(pad=TickSpace, labelsize=TickSize)
ax.yaxis.set_tick_params(pad=TickSpace, labelsize=TickSize)
ax.zaxis.set_tick_params(pad=TickSpace, labelsize=TickSize)
# Save and close
fig.savefig('Movie/movie%d.png' % frame_num )
plt.close(fig)
def save_plot2():
fig2 = plt.figure(figsize=(height,width))
plt.title('TOP VIEW - Shoot length: %.2f cm, Shoot time: %.4f min' % (Lt*100,t*Tau_s/60),fontsize=TitleSize)
plt.plot(Xtip[-TipPts:],Ytip[-TipPts:],'.')
plt.xlabel("$\mathbf{e}_3$ [m]",fontsize=LabelSize)
plt.ylabel("$\mathbf{e}_1$ [m]",fontsize=LabelSize)
plt.tick_params(axis='both', labelsize=TickSize)
fig2.savefig('Movie2/Tip%d.png' % frame_num )
plt.close(fig2)
############################################# POST-PROCESSING FUNCTIONS ##########################################
def reconstruct():
# Function to reconstruct the shape of the rod, starting from the computed angles
global x,y,z,Xtip,Ytip,Ztip,D1,D2,D3,H,old,young
ds = np.array([j-i for i, j in zip(space_ref[:-1], space_ref[1:])]) # reference space-steps
# Evaluate fields at the mesh points
space_cur, stretch_num, psi_num, phi_num, chi_num, thetaH_num, alphaH_num = np.empty(0), np.empty(0), np.empty(0), np.empty(0), np.empty(0), np.empty(0), np.empty(0)
for point in space_ref:
space_cur = np.append(space_cur,s(point,t)) # corresponding mesh in the current configuration
stretch_num = np.append(stretch_num,stretch_fun(point,t))
psi_num = np.append(psi_num,psi(point))
phi_num = np.append(phi_num,phi(point))
chi_num = np.append(chi_num,chi(point))
thetaH_num = np.append(thetaH_num,thetaH(point))
alphaH_num = np.append(alphaH_num,alphaH(point))
h1_num = np.cos(thetaH_num)
h2_num = np.sin(thetaH_num)*np.cos(alphaH_num)
h3_num = np.sin(thetaH_num)*np.sin(alphaH_num)
# Integrate the d3 director
D31 = stretch_num*np.sin(psi_num)*np.cos(phi_num)
D32 = stretch_num*np.sin(psi_num)*np.sin(phi_num)
D33 = stretch_num*np.cos(psi_num)
for n in range(1,2*nx+1):
x[n] = x[n-1]+ds[n-1]*D33[n] # integration along e3
y[n] = y[n-1]+ds[n-1]*D31[n] # integration along e1
z[n] = z[n-1]+ds[n-1]*D32[n] # integration along e2
# Update tip coordinates
Xtip = np.append(Xtip,x[-1])
Ytip = np.append(Ytip,y[-1])
Ztip = np.append(Ztip,z[-1])
# Update directors and statoliths direction, at selected points
CHI = chi_num[ArrowIndex]
PSI = psi_num[ArrowIndex]
PHI = phi_num[ArrowIndex]
D1 = (np.cos(CHI)*np.cos(PSI)*np.cos(PHI)-np.sin(CHI)*np.sin(PHI),np.cos(CHI)*np.cos(PSI)*np.sin(PHI)+np.sin(CHI)*np.cos(PHI),-np.cos(CHI)*np.sin(PSI))
D2 = (-np.sin(CHI)*np.cos(PSI)*np.cos(PHI)-np.cos(CHI)*np.sin(PHI),-np.sin(CHI)*np.cos(PSI)*np.sin(PHI)+np.cos(CHI)*np.cos(PHI),np.sin(CHI)*np.sin(PSI))
D3 = (np.sin(PSI)*np.cos(PHI),np.sin(PSI)*np.sin(PHI),np.cos(PSI))
H = h1_num[ArrowIndex]*D1+h2_num[ArrowIndex]*D2+h3_num[ArrowIndex]*D3
# Update old and young
old = space_cur < (Lt-lg)
young = [not(x) for x in old]
def save_data():
# Save solutions
for i in range(max(Nr,Nrp)):
U.assign(SolutionArray[i])
output_file = HDF5File(mesh.mpi_comm(), 'Data/Frame%d/solution%d.h5' % (frame_num,i), "w")
output_file.write(U, 'Data/Frame%d/solution%d' % (frame_num,i))
output_file.close()
################################################# SOLVER FUNCTIONS ################################################
def Picard(maxiter=250,tol=1.0E-11):
# maxiter is the max no of iterations allowed
global U_k, U, a, L, bcs, chi, phi, psi, u1S, u2S, u3S, w1, w2, wp1, wp2, alphaH, thetaH
eps = 1.0 # error measure ||u-u_k||
Nit = 0 # iteration counter
while eps > tol and Nit < maxiter:
Nit += 1
solve(a == L, U, bcs)
chi, phi, psi, u1S, u2S, u3S, w1, w2, wp1, wp2, alphaH, thetaH = U.split()
diff = U.vector().get_local() - U_k.vector().get_local()
eps = np.linalg.norm(diff, ord=np.Inf)
print('iter=%d: norm=%g' % (Nit, eps))
U_k.assign(U) # update for next iteration
def Newton(maxiter=50,tol=1E-14):
global U, F, bcs, J
problem = NonlinearVariationalProblem(F, U, bcs, J)
solver = NonlinearVariationalSolver(problem)
prm = solver.parameters
prm["newton_solver"]["absolute_tolerance"] = tol
prm["newton_solver"]["relative_tolerance"] = tol
prm["newton_solver"]["maximum_iterations"] = maxiter
solver.solve()
########################################### GLOBAL CONSTANTS ##############################################
g = 9.81 # m/s^2 - gravity acceleration
########################################### MODEL PARAMETERS ##############################################
# Time-scales
Tau_a = 2*60 # s - time-scale for statoliths avalanche dynamics
Tau_s = 12*60 # s - time-scale for nondimensionalizing the equations
Tau_m = 12*60 # s - Memory time for gravitropism
Tau_mp = 12*60 # s - Memory time for proprioception
Tau_r = 12*60 # s - Reaction time for gravitropism
Tau_rp = 12*60 # s - Reaction time for proprioception
Tau_g = 20*60*60 # s - Growth time for the upper part
Tau_e = 24*60 # s - Time period for endogenous oscillator
Tau_l = 6*24*60*60 # s - Maturation/Lignification time
# Geometry
L0 = 6.8E-2 # m - initial length
lg = 6E-2 # m - growing zone
R = 5E-4 # m - radius
A = np.pi*R**2 # m^2 - cross sectional area
# BioMechanical properties of plant tissue
Lig = 200 # rod stiffening ratio
EY = 10E+6 # N/m^2 - Young Modulus
I = np.pi*R**4/4 # m^4 - Second moment of inertia
B0 = EY*I # Bending stiffness
Bmax = Lig*EY*I
nu = 0.5 # Poisson's ratio
mu = EY/(2*(1+nu))# shear modulus
J = 2*I # parameter depending on the cross-sectional shape
Jmax = Lig*J
rho = 1000 # Kg/m^3 - density
alpha= 0 # internal oscillator sensitivity
beta = 0.8 # gravitropic sensitivity
eta = 20 # proprioceptic sensitivity
# Loads
Q1 = 1E-6 # N/m - distributed load along e1 per unit length
Q2 = -rho*g*A # N/m - distributed load along e2 per unit length
Q3 = 0 # N/m - distributed load along e3 per unit length
al1 = 1E-7 # N - apical load along e1
al2 = 0 # N - apical load along e2
al3 = 1E-7 # N - apical load along e3
########################################### REFERENCE MAPS ##############################################
# Critical time of a material point: when a material point exists the growth zone
def Tc(S0):
if S0==L0:
return np.Infinity
else:
return np.log(lg/(L0-S0))*Tau_g/Tau_s
Tc0 = Tc(0) # critical time of the left end point
# Shoot length at dimensionless time t
def length_fun(t):
if Tau_g == np.Infinity:
return L0
elif t <= Tc0:
return L0*np.exp(t*Tau_s/Tau_g)
else:
return max(L0,lg)+lg*(t-max(Tc0,0))*Tau_s/Tau_g
# Critical dimensionless time of a spatial point: when a point at position s at time t exited the growth zone
def Tc_(s,t):
return t+(s+lg-length_fun(t))*Tau_g/(Tau_s*lg)
# Reference map: material point corresponding to a point s at dimensionless time t
def S0(s,t):
if Tau_g == np.Infinity:
return s
elif 0<=s<(L0-lg):
return s
elif (L0-lg)<=s<(length_fun(t)-lg):
return (s-length_fun(Tc_(s,t)))*np.exp(-Tc_(s,t)*Tau_s/Tau_g)+L0
else:
return (s-length_fun(t))*np.exp(-t*Tau_s/Tau_g)+L0
# Inverse of the reference map - dimensionless time t
def s(S0,t):
if S0<=(L0-lg):
return S0
elif (S0>(L0-lg) and t>=Tc(S0)):
return length_fun(Tc(S0))-lg
else:
return length_fun(t)-(L0-S0)*np.exp(t*Tau_s/Tau_g)
# Stretch function in terms of the material reference configuration - dimensionless time t
def stretch_fun(S0,t):
if Tau_g == np.Infinity:
return 1
elif S0<=(L0-lg):
return 1
elif (S0>(L0-lg) and t>=Tc(S0)):
return lg/(L0-S0)
else:
return np.exp(t*Tau_s/Tau_g)
########################################### FEM IMPLEMENTATION ##############################################
# Create mesh
nx = 2**8
mesh = IntervalMesh(2*nx,0,L0)
space_ref= mesh.coordinates()
# Time-step
Nr = 2**6 # number of time-steps for gravitropic delay
dt = (Tau_r/Tau_s)/Nr # time-step
Nrp = int((Tau_rp/Tau_s)/dt) # number of time-steps for proprioceptic delay
# Function space
P1 = FiniteElement('P',interval, 1)
element = MixedElement([P1, P1, P1, P1, P1, P1, P1, P1, P1, P1, P1, P1])
V = FunctionSpace(mesh, element)
# Boundary conditions
def boundary_D_l(x, on_boundary):
return on_boundary and near(x[0],0,1E-14)
bc_chi = DirichletBC(V.sub(0), Constant(0), boundary_D_l)
bc_phi = DirichletBC(V.sub(1), Constant(np.pi/2), boundary_D_l)
bc_psi = DirichletBC(V.sub(2), Constant(np.pi/2), boundary_D_l)
bcs = [bc_chi,bc_phi,bc_psi]
# Initial conditions for statolith direction
alphaH0 = Constant(np.pi/2)
thetaH0 = Constant(np.pi/2)
# Define expressions used in variational forms
h = Expression('h', h=dt, degree=1)
alpha = Constant(alpha)
beta = Constant(beta)
eta = Constant(eta)
q1 = Constant(Q1)
q2 = Constant(Q2)
q3 = Constant(Q3)
p1 = Expression('al',al=al1,t=0,degree=1)
p2 = Expression('al',al=al2,t=0,degree=1)
p3 = Expression('al',al=al3,t=0,degree=1)
B = Expression('x[0]>=(L0-lg*exp(-t*ts/tg)) ? B0:(Bmax-(Bmax-B0)*exp(-(x[0]<(L0-lg) ? (t*ts+(L0-lg-x[0])*tg/lg):(t*ts-tg*log10(lg/(L0-x[0]))/log10(exp(1))))/tl))',
tl=Tau_l, B0=B0, Bmax=Bmax, L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=0., degree=1)
muJ = Expression('x[0]>=(L0-lg*exp(-t*ts/tg)) ? J0:(Jmax-(Jmax-J0)*exp(-(x[0]<(L0-lg) ? (t*ts+(L0-lg-x[0])*tg/lg):(t*ts-tg*log10(lg/(L0-x[0]))/log10(exp(1))))/tl))',
tl=Tau_l, J0=mu*J, Jmax=mu*Jmax, L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=0., degree=1)
r = Constant(R)
Tg = Constant(Tau_g)
Tm = Constant(Tau_m)
Tmp = Constant(Tau_mp)
Ts = Constant(Tau_s)
Ta = Constant(Tau_a)
stretch= Expression('x[0]<(L0-lg) ? 1 : (x[0]>=(L0-lg*exp(-t*ts/tg)) ? exp(t*ts/tg) : lg/(L0-x[0]))',
L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=0., degree=1) # stretch
stretch_rp = Expression('x[0]<(L0-lg) ? 1 : (x[0]>=(L0-lg*exp(-t*ts/tg)) ? exp(t*ts/tg) : lg/(L0-x[0]))',
L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=-Tau_rp/Tau_s, degree=1) # stretch evaluated at previous times, according to the proprioceptive delay
s = Expression('x[0]<(L0-lg) ? x[0] : (x[0]>=(L0-lg*exp(-t*ts/tg)) ? (Lt-(L0-x[0])*exp(t*ts/tg)):((x[0]<=(L0-lg*exp(-tc0*ts/tg)) ? (L0*exp(log10(lg/(L0-x[0]))/log10(exp(1)))):((L0>lg ? L0:lg)+lg*(log10(lg/(L0-x[0]))/log10(exp(1))-(tc0*ts>0 ? tc0*ts/tg:0))))-lg) )',
L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=0., Lt=L0, tc0=Tc0, degree=1) # motion
rg = Expression('x[0]>=(L0-lg*exp(-t*ts/tg)) ? 1 : 0', L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=0., degree=1) # 1 in the growing region, 0 otherwise
length = Expression('t<=tc0 ? (L0*exp(t*ts/tg)):((L0>lg ? L0:lg)+lg*(t*ts-(tc0*ts>0 ? tc0*ts:0))/tg)',
L0=L0, lg=lg, ts=Tau_s, tg=Tau_g, t=0., tc0=Tc0, degree=1) # current total length
v1 = Expression('cos(omega*t)', omega=2*np.pi*Tau_s/Tau_e, t=0, degree=1) # internal oscillator component (along d1)
v2 = Expression('sin(omega*t)', omega=2*np.pi*Tau_s/Tau_e, t=0, degree=1) # internal oscillator component (along d2)
# Define variational problem in the reference domain [0,L0]
test_chi, test_phi, test_psi, test_u1S, test_u2S, test_u3S, test_w1, test_w2, test_wp1, test_wp2, test_alphaH, test_thetaH = TestFunctions(V)
# Retarded unknowns
U_r = Function(V)
chi_r, phi_r, psi_r, u1S_r, u2S_r, u3S_r, w1_r, w2_r, wp1_r, wp2_r, alphaH_r, thetaH_r = split(U_r) # retarded functions for gravitropic delay
U_rp = Function(V)
chi_rp, phi_rp, psi_rp, u1S_rp, u2S_rp, u3S_rp, w1_rp, w2_rp, wp1_rp, wp2_rp, alphaH_rp, thetaH_rp = split(U_rp) # retarded functions for gravitropic delay
# Previous step unknowns
U_n = Function(V)
chi_n, phi_n, psi_n, u1S_n, u2S_n, u3S_n, w1_n, w2_n, wp1_n, wp2_n, alphaH_n, thetaH_n = split(U_n)
# Unknowns
U = Function(V)
chi, phi, psi, u1S, u2S, u3S, w1, w2, wp1, wp2, alphaH, thetaH = split(U)
# Functional defined by using backward Euler in time
n1 = p1+q1*(length-s) # resultant contact force - e1 component
n2 = p2+q2*(length-s) # resultant contact force - e2 component
n3 = p3+q3*(length-s) # resultant contact force - e3 component
u1 = (psi.dx(0)*sin(chi)-phi.dx(0)*cos(chi)*sin(psi))/stretch # flexural strain 1
U1 = u1-u1S
u2 = (psi.dx(0)*cos(chi)+phi.dx(0)*sin(chi)*sin(psi))/stretch # flexural strain 2
U2 = u2-u2S
u3 = (chi.dx(0)+phi.dx(0)*cos(psi))/stretch # torsional strain
U3 = u3-u3S
m1 = B*(cos(psi)*cos(phi)*(U1*cos(chi)-U2*sin(chi))-sin(phi)*(U1*sin(chi)+U2*cos(chi))) \
+ muJ*U3*sin(psi)*cos(phi) # resultant contact couple - e1 component
m2 = B*(cos(psi)*sin(phi)*(U1*cos(chi)-U2*sin(chi))+cos(phi)*(U1*sin(chi)+U2*cos(chi))) \
+ muJ*U3*sin(psi)*sin(phi) # resultant contact couple - e2 component
m3 = -B*sin(psi)*(U1*cos(chi)-U2*sin(chi)) \
+ muJ*U3*cos(psi) # resultant contact couple - e3 component
u1_rp = (psi_rp.dx(0)*sin(chi_rp)-phi_rp.dx(0)*cos(chi_rp)*sin(psi_rp))/stretch_rp
u2_rp = (psi_rp.dx(0)*cos(chi_rp)+phi_rp.dx(0)*sin(chi_rp)*sin(psi_rp))/stretch_rp
H1_r = cos(thetaH_r)
H2_r = sin(thetaH_r)*cos(alphaH_r)
RHS_alpha = cos(alphaH)*sin(psi)*sin(phi)+(cos(psi)*sin(chi)*sin(phi)-cos(chi)*cos(phi))*sin(alphaH)
RHS_theta = cos(thetaH)*( cos(alphaH)*(cos(chi)*cos(phi)-cos(psi)*sin(chi)*sin(phi))+sin(alphaH)*sin(psi)*sin(phi) )-sin(thetaH)*( cos(phi)*sin(chi)+cos(chi)*cos(psi)*sin(phi) )
F = - m1*test_chi.dx(0)*dx + stretch*(n3*sin(psi)*sin(phi)-n2*cos(psi))*test_chi*dx \
- m2*test_psi.dx(0)*dx + stretch*(n1*cos(psi)-n3*sin(psi)*cos(phi))*test_psi*dx \
- m3*test_phi.dx(0)*dx + stretch*(n2*sin(psi)*cos(phi)-n1*sin(psi)*sin(phi))*test_phi*dx \
+ ((u1S-u1S_n)*r*Tg/Ts-h*rg*alpha*v1-h*rg*beta*w1-h*rg*eta*r*wp1)*test_u1S*dx \
+ ((u2S-u2S_n)*r*Tg/Ts-h*rg*alpha*v2-h*rg*beta*w2-h*rg*eta*r*wp2)*test_u2S*dx \
+ u3S*test_u3S*dx \
+ ((w1-w1_n)*Tm/Ts+h*(w1+H2_r))*test_w1*dx \
+ ((w2-w2_n)*Tm/Ts+h*(w2-H1_r))*test_w2*dx \
+ ((wp1-wp1_n)*Tmp/Ts+h*(wp1+u1_rp))*test_wp1*dx \
+ ((wp2-wp2_n)*Tmp/Ts+h*(wp2+u2_rp))*test_wp2*dx \
+ ((alphaH-alphaH_n)*sin(thetaH)-RHS_alpha*h*Ts/Ta)*test_alphaH*dx \
+ ((thetaH-thetaH_n)-RHS_theta*h*Ts/Ta)*test_thetaH*dx
du = TrialFunction(V)
J = derivative(F,U,du)
# LINEAR PROBLEM FOR PICARD
chi_t, phi_t, psi_t, u1S_t, u2S_t, u3S_t, w1_t, w2_t, wp1_t, wp2_t, alphaH_t, thetaH_t = TrialFunctions(V)
U_k = Function(V)
chi_k, phi_k, psi_k, u1S_k, u2S_k, u3S_k, w1_k, w2_k, wp1_k, wp2_k, alphaH_k, thetaH_k = split(U_k)
u1_k = psi_t.dx(0)*sin(chi_k)-phi_t.dx(0)*cos(chi_k)*sin(psi_k)/stretch
U1_k = u1_k-u1S_t
u2_k = psi_t.dx(0)*cos(chi_k)+phi_t.dx(0)*sin(chi_k)*sin(psi_k)/stretch
U2_k = u2_k-u2S_t
u3_k = chi_t.dx(0)+phi_t.dx(0)*cos(psi_k)/stretch
U3_k = u3_k-u3S_t
m1_k = B*(cos(psi_k)*cos(phi_k)*(U1_k*cos(chi_k)-U2_k*sin(chi_k)) -sin(phi_k)*(U1_k*sin(chi_k)+U2_k*cos(chi_k))) \
+ muJ*U3_k*sin(psi_k)*cos(phi_k)
m2_k = B*(cos(psi_k)*sin(phi_k)*(U1_k*cos(chi_k)-U2_k*sin(chi_k))+cos(phi_k)*(U1_k*sin(chi_k)+U2_k*cos(chi_k))) \
+ muJ*U3_k*sin(psi_k)*sin(phi_k)
m3_k = -B*sin(psi_k)*(U1_k*cos(chi_k)-U2_k*sin(chi_k)) \
+ muJ*U3_k*cos(psi_k)
RHS_alpha_k = cos(alphaH_k)*sin(psi_k)*sin(phi_k)+(cos(psi_k)*sin(chi_k)*sin(phi_k)-cos(chi_k)*cos(phi_k))*sin(alphaH_k)
RHS_theta_k = cos(thetaH_k)*( cos(alphaH_k)*(cos(chi_k)*cos(phi_k)-cos(psi_k)*sin(chi_k)*sin(phi_k))+sin(alphaH_k)*sin(psi_k)*sin(phi_k) ) -sin(thetaH_k)*( cos(phi_k)*sin(chi_k)+cos(chi_k)*cos(psi_k)*sin(phi_k) )
a = - m1_k*test_chi.dx(0)*dx \
- m2_k*test_psi.dx(0)*dx \
- m3_k*test_phi.dx(0)*dx \
+ (u1S_t*r*Tg/Ts-h*rg*beta*w1_t-h*rg*eta*r*wp1_t)*test_u1S*dx \
+ (u2S_t*r*Tg/Ts-h*rg*beta*w2_t-h*rg*eta*r*wp2_t)*test_u2S*dx \
+ u3S_t*test_u3S*dx \
+ (Tm/Ts+h)*w1_t*test_w1*dx \
+ (Tm/Ts+h)*w2_t*test_w2*dx \
+ (Tmp/Ts+h)*wp1_t*test_wp1*dx \
+ (Tmp/Ts+h)*wp2_t*test_wp2*dx \
+ alphaH_t*sin(thetaH_k)*test_alphaH*dx \
+ thetaH_t*test_thetaH*dx
L = - stretch*(n3*sin(psi_k)*sin(phi_k)-n2*cos(psi_k))*test_chi*dx \
- stretch*(n1*cos(psi_k)-n3*sin(psi_k)*cos(phi_k))*test_psi*dx \
- stretch*(n2*sin(psi_k)*cos(phi_k)-n1*sin(psi_k)*sin(phi_k))*test_phi*dx \
+ ((r*Tg/Ts)*u1S_n+h*rg*alpha*v1)*test_u1S*dx \
+ ((r*Tg/Ts)*u2S_n+h*rg*alpha*v2)*test_u2S*dx \
+ (w1_n*Tm/Ts-h*H2_r)*test_w1*dx \
+ (w2_n*Tm/Ts+h*H1_r)*test_w2*dx \
+ (wp1_n*Tmp/Ts-h*u1_rp)*test_wp1*dx \
+ (wp2_n*Tmp/Ts-h*u2_rp)*test_wp2*dx \
+ (alphaH_n*sin(thetaH_k)+RHS_alpha_k*h*Ts/Ta)*test_alphaH*dx \
+ (thetaH_n+RHS_theta_k*h*Ts/Ta)*test_thetaH*dx
################################################## INITIALIZE STORAGE ##################################################
t = 0
Lt = length_fun(t) # current length
ArrowIndex = np.array([int(i) for i in np.linspace(0,2*nx,5)]) # material points where to plot directors
# Initialize rod coordinates
x = np.zeros(2*nx+1)
y = np.zeros(2*nx+1)
z = np.zeros(2*nx+1)
# Initialize tip Coordinates
Xtip, Ytip, Ztip = np.empty(0), np.empty(0), np.empty(0)
################################################## INITIAL ELASTIC EQUILIBRIUM ###########################################
# SOLVE INITIAL GUESS: STRAIGHT CONFIGURATION
F0 = chi.dx(0)*test_chi*dx \
+ psi.dx(0)*test_psi*dx \
+ phi.dx(0)*test_phi*dx \
+ u1S*test_u1S*dx \
+ u2S*test_u2S*dx \
+ u3S*test_u3S*dx \
+ w1*test_w1*dx \
+ w2*test_w2*dx \
+ wp1*test_wp1*dx \
+ wp2*test_wp2*dx \
+ (alphaH-alphaH0)*test_alphaH*dx \
+ (thetaH-thetaH0)*test_thetaH*dx
J0 = derivative(F0,U,du)
problem0 = NonlinearVariationalProblem(F0, U, bcs, J0)
solver0 = NonlinearVariationalSolver(problem0)
solver0.solve()
# SOLVE INITIAL ELASTIC EQUILIBRIUM
F_n = - m1*test_chi.dx(0)*dx + stretch*(n3*sin(psi)*sin(phi)-n2*cos(psi))*test_chi*dx \
- m2*test_psi.dx(0)*dx + stretch*(n1*cos(psi)-n3*sin(psi)*cos(phi))*test_psi*dx \
- m3*test_phi.dx(0)*dx+stretch*(n2*sin(psi)*cos(phi)-n1*sin(psi)*sin(phi))*test_phi*dx \
+ u1S*test_u1S*dx \
+ u2S*test_u2S*dx \
+ u3S*test_u3S*dx \
+ w1*test_w1*dx \
+ w2*test_w2*dx \
+ wp1*test_wp1*dx \
+ wp2*test_wp2*dx \
+ (alphaH-alphaH0)*test_alphaH*dx \
+ (thetaH-thetaH0)*test_thetaH*dx
J_n = derivative(F_n,U,du)
problem_n = NonlinearVariationalProblem(F_n, U, bcs, J_n)
solver_n = NonlinearVariationalSolver(problem_n)
solver_n.solve()
# Initialization array of solutions at previous steps
SolutionArray = []
for n in range(max(Nr,Nrp)):
SolutionArray.append(U.copy(deepcopy=True))
chi, phi, psi, u1S, u2S, u3S, w1, w2, wp1, wp2, alphaH, thetaH = U.split()
# Curve reconstruction from the angles
frame_num = 0 # Initialize frame counter
reconstruct()
# Plot and Save
proc = mp.Process(target=save_plot)
proc.daemon = True
proc.start()
proc.join()
proc2 = mp.Process(target=save_plot2)
proc2.daemon = True
proc2.start()
proc2.join()
#################################################### SOLVE ########################################################
while t<10000:
# Update frame counter
frame_num += 1
# Update current time and all functions depending on time
t += dt
stretch.t = t
stretch_rp.t = t-Tau_rp/Tau_s
s.t = t
B.t = t
muJ.t = t
Lt = length_fun(t) # current length
s.Lt = Lt
rg.t = t
length.t = t
p1.t = t
p2.t = t
p3.t = t
v1.t = t # internal oscillator component along d1
v2.t = t # internal oscillator component along d2
# Update variables
U_n.assign(SolutionArray[-1]) # update the previous step
U_r.assign(SolutionArray[-Nr]) # update the retarded solution for gravitropism
U_rp.assign(SolutionArray[-Nrp]) # update the retarded solution for proprioception
start = time.time()
try:
# Solve the problem through nonlinear FEniCS solver
Newton()
except:
# Use Picard's method
print('Newton did not converge: trying Picard first...')
U.assign(SolutionArray[-1])
U_k.assign(SolutionArray[-1])
Picard()
# Solve the problem through nonlinear FEniCS solver
Newton(50,1E-10)
end = time.time()
print('Frame: %d, Shoot Time: %g min, Shoot Length: %g cm, Computation Time: %g sec' % (frame_num, t*Tau_s/60, Lt*100, end-start))
# Update
chi, phi, psi, u1S, u2S, u3S, w1, w2, wp1, wp2, alphaH, thetaH = U.split() # update variables
SolutionArray = SolutionArray[1:] # remove the first element
SolutionArray.append(U.copy(deepcopy=True)) # add the last solution
##### PLOT AND SAVE ####
# Curve reconstruction
reconstruct()
# Plot and Save
proc = mp.Process(target=save_plot)
proc.daemon = True
proc.start()
proc.join()
proc2 = mp.Process(target=save_plot2)
proc2.daemon = True
proc2.start()
proc2.join()
if np.mod(frame_num,max(Nr,Nrp))==0:
save_data()
print('Saving data...')