39
40
41
43 USE elbufdef_mod
44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "parit_c.inc"
52#include "com20_c.inc"
53#include "com08_c.inc"
54#include "scr02_c.inc"
55#include "scr18_c.inc"
56
57
58
59 INTEGER, INTENT(IN) :: NFT
60 INTEGER :: NEL,IMAT,NDDL,ITASK
61 INTEGER, DIMENSION(NEL) :: NC1,NC2,NC3,NC4
62 my_real,
DIMENSION(NEL,NDDL),
INTENT(INOUT)::
63 . var_reg
64 my_real,
DIMENSION(NEL),
INTENT(IN) ::
65 .
area,off,px1,py1,px2,py2,thk,le,thk0,area0
67 . dt2t
68 TYPE(NLOCAL_STR_),TARGET :: NLOC_DMG
69 TYPE(BUF_NLOC_) ,TARGET :: BUFNL
70
71
72
73 INTEGER I,II,K,N1,N2,N3,N4,L_NLOC,J,NDNOD
75 . dx, dy, dz, l2,ntn,ntn_unl,ntn_vnl,xi,ntvar,a,zeta,
76 . b1,b2,b3,b4,sspnl,nth1,nth2,bth1,bth2,k1,k2,k12,
77 . dtnl_th,dtnl,le_max,maxstif,dtnod,dt2p
78 my_real,
DIMENSION(:,:),
ALLOCATABLE ::
79 . f1,f2,f3,f4,sti1,sti2,sti3,sti4
80 my_real,
DIMENSION(:) ,
ALLOCATABLE ::
81 . btb11,btb12,btb22,vol
82 INTEGER, DIMENSION(:), ALLOCATABLE ::
83 . POS1,POS2,POS3,POS4
84 my_real,
POINTER,
DIMENSION(:) ::
85 . vnl,fnl,unl,stifnl,mass,mass0,vnl0
86 my_real,
POINTER,
DIMENSION(:,:) ::
87 . massth,unlth,vnlth,fnlth
88 my_real,
DIMENSION(:,:),
ALLOCATABLE ::
89 . stifnlth,dtn
90
91
92 my_real,
PARAMETER :: csta = 40.0d0
93
94 my_real,
PARAMETER :: cdamp = 0.7d0
95
96
97
98
99
100 l2 = nloc_dmg%LEN(imat)**2
101 xi = nloc_dmg%DAMP(imat)
102 zeta = nloc_dmg%DENS(imat)
103 sspnl = nloc_dmg%SSPNL(imat)
104 l_nloc = nloc_dmg%L_NLOC ! length of non-local tables
105 le_max = nloc_dmg%LE_MAX(imat)
106
107 ALLOCATE(f1(nel,nddl),f2(nel,nddl),f3(nel,nddl),f4(nel,nddl))
108
109 IF (nodadt > 0) THEN
110
111 ALLOCATE(sti1(nel,nddl),sti2(nel,nddl),sti3(nel,nddl),sti4(nel,nddl))
112
113 mass => nloc_dmg%MASS(1:l_nloc)
114
115 mass0 => nloc_dmg%MASS0(1:l_nloc)
116 ENDIF
117 ALLOCATE(btb11(nel),btb12(nel),btb22(nel),vol(nel),
118 . pos1(nel),pos2(nel),pos3(nel),pos4(nel))
119
120 vnl => nloc_dmg%VNL(1:l_nloc)
121 vnl0 => nloc_dmg%VNL_OLD(1:l_nloc)
122 unl => nloc_dmg%UNL(1:l_nloc)
123
124 ntn = four*four
125
126
127
128
129
130# include "vectorize.inc"
131 DO i=1,nel
132
133
134 n1 = nloc_dmg%IDXI(nc1(i))
135 n2 = nloc_dmg%IDXI(nc2(i))
136 n3 = nloc_dmg%IDXI(nc3(i))
137 n4 = nloc_dmg%IDXI(nc4(i))
138
139
140 pos1(i) = nloc_dmg%POSI(n1)
141 pos2(i) = nloc_dmg%POSI(n2)
142 pos3(i) = nloc_dmg%POSI(n3)
143 pos4(i) = nloc_dmg%POSI(n4)
144
145
146 vol(i) =
area(i)*thk(i)
147
148
149 btb11(i) = px1(i)**2 + py1(i)**2
150 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i)
151 btb22(i) = px2(i)**2 + py2(i)**2
152
153 ENDDO
154
155
156
157
158
159 IF ((nddl > 1).AND.(l2>zero)) THEN
160
161
162 IF (nddl > 2) THEN
163 IF (nodadt > 0) THEN
164 ALLOCATE(stifnlth(nel,nddl+1))
165 ALLOCATE(dtn(nel,nddl+1))
166 ENDIF
167 ndnod = nddl+1
168 ELSE
169 IF (nodadt > 0) THEN
170 ALLOCATE(stifnlth(nel,nddl))
171 ALLOCATE(dtn(nel,nddl))
172 ENDIF
173 ndnod = nddl
174 ENDIF
175
176
177 massth => bufnl%MASSTH(1:nel,1:ndnod)
178 unlth => bufnl%UNLTH(1:nel,1:ndnod)
179 vnlth => bufnl%VNLTH(1:nel,1:ndnod)
180 fnlth => bufnl%FNLTH(1:nel,1:ndnod)
181
182 DO k = 1,ndnod
183 DO i = 1,nel
184
185 fnlth(i,k) = zero
186
187 IF (nodadt > 0) THEN
188 stifnlth(i,k) = em20
189 ENDIF
190 ENDDO
191 ENDDO
192
193
194 DO k = 1, nddl
195
196
197 IF ((nddl==2).AND.(k==2)) THEN
198 nth1 = (z0(k,nddl) - zth(k,nddl)) / (zth(k-1,nddl) - zth(k,nddl))
199 nth2 = (z0(k,nddl) - zth(k-1,nddl)) / (zth(k,nddl) - zth(k-1,nddl))
200 ELSE
201 nth1 = (z0(k,nddl) - zth(k+1,nddl)) / (zth(k,nddl) - zth(k+1,nddl))
202 nth2 = (z0(k,nddl) - zth(k,nddl)) / (zth(k+1,nddl) - zth(k,nddl))
203 ENDIF
204
205
206 DO i = 1,nel
207
208 IF ((nddl==2).AND.(k==2)) THEN
209 bth1 = (one/(zth(k-1,nddl) - zth(k,nddl)))*(one/thk(i))
210 bth2 = (one/(zth(k,nddl) - zth(k-1,nddl)))*(one/thk(i))
211 ELSE
212 bth1 = (one/(zth(k,nddl) - zth(k+1,nddl)))*(one/thk(i))
213 bth2 = (one/(zth(k+1,nddl) - zth(k,nddl)))*(one/thk(i))
214 ENDIF
215
216
217 k1 = l2*(bth1**2) + nth1**2
218 k12 = l2*(bth1*bth2)+ (nth1*nth2)
219 k2 = l2*(bth2**2) + nth2**2
220
221
222 IF ((nddl==2).AND.(k==2)) THEN
223 fnlth(i,k-1) = fnlth(i,k-1) + (k1*unlth(i,k-1) + k12*unlth(i,k)
224 . + xi*((nth1**2)*vnlth(i,k-1)
225 . + (nth1*nth2)*vnlth(i,k
226 . - (nth1*var_reg(i,k)))*vol(i)*wf(k,nddl)
227 fnlth(i,k) = fnlth(i,k) + (k12*unlth(i,k-1) + k2*unlth(i,k)
228 . + xi*(nth1*nth2*vnlth(i,k-1)
229 . + (nth2**2)*vnlth(i,k))
230 . - nth2*var_reg(i,k))*vol(i)*wf(k,nddl)
231 ELSE
232 fnlth(i,k) = fnlth(i,k) + (k1*unlth(i,k) + k12*unlth(i,k+1)
233 . + xi*((nth1**2)*vnlth(i,k)
234 . + (nth1*nth2)*vnlth(i,k+1))
235 . - (nth1*var_reg(i,k)))*vol(i)*wf(k,nddl
236 fnlth(i,k+1) = fnlth(i,k+1) + (k12*unlth(i,k)
237 . + xi*(nth1*nth2*vnlth(i,k)
238 . + (nth2**2)*vnlth(i,k+1))
239 . - nth2*var_reg(i,k))*vol(i)*wf(k,nddl)
240 ENDIF
241
242
243 IF (nodadt > 0) THEN
244 IF ((nddl==2).AND.(k==2)) THEN
245 stifnlth(i,k-1) = stifnlth(i,k-1) +
max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
246 stifnlth(i,k) = stifnlth(i,k) +
max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
247 ELSE
248 stifnlth(i,k) = stifnlth(i,k) +
max(abs(k1)+abs(k12),abs(k12)+abs(k2))*vol(i)*wf(k,nddl)
249 stifnlth(i,k+1) = stifnlth(i,k+1) +
max
250 ENDIF
251 ENDIF
252
253 ENDDO
254 ENDDO
255
256
257 IF (nodadt > 0) THEN
258
259
260 dtnod = ep20
261 DO k = 1,ndnod
262 DO i = 1,nel
263 dtn(i,k) = dtfac1(11)*cdamp*sqrt(two*massth(i,k)/
max(stifnlth(i,k),em20))
264 dtnod =
min(dtn(i,k),dtnod)
265 ENDDO
266 ENDDO
267
268
269 IF ((idtmin(11)==3).OR.(idtmin(11)==4).OR.(idtmin(11)==8)) THEN
270
271 IF (dtnod < dtmin1(11)*(sqrtTHEN
272 DO k = 1,ndnod
273 DO i = 1,nel
274 IF (dtn(i,k) < dtmin1(11)) THEN
275 dt2p = dtmin1(11)/(dtfac1(11)*cdamp
276 massth(i,k) =
max(massth(i,k),csta*half*stifnlth(i,k)*dt2p*dt2p*onep00001)
277 ENDIF
278 ENDDO
279 ENDDO
280 ENDIF
281 dtnod = dtmin1(11)*(sqrt(csta))
282 ENDIF
283
284
285 IF (dtnod < dt2t) THEN
286 dt2t =
min(dt2t,dtnod)
287 ENDIF
288 ENDIF
289
290 DO k = 1,ndnod
291 DO i = 1,nel
292
293 vnlth(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*dt12
294 ENDDO
295 ENDDO
296
297 DO k = 1,ndnod
298 DO i = 1,nel
299
300 unlth(i,k) = unlth(i,k) + vnlth(i,k)*dt1
301 ENDDO
302 ENDDO
303
304
305 DO k = 1, nddl
306
307 IF ((nddl==2).AND.(k==2)) THEN
308 nth1 = (z0(k,nddl) - zth(k,nddl))/(zth(k-1,nddl) - zth(k,nddl))
309 nth2 = (z0(k,nddl) - zth(k-1,nddl)) /(zth(k,nddl) - zth(k-1,nddl))
310 ELSE
311 nth1 = (z0(k,nddl) - zth(k+1,nddl))/(zth(k,nddl) - zth(k+1,nddl))
312 nth2 = (z0(k,nddl) - zth(k,nddl)) /(zth(k+1,nddl) - zth(k,nddl))
313 ENDIF
314
315 DO i = 1,nel
316
317 IF ((nddl==2).AND.(k==2)) THEN
318 var_reg(i,k) = nth1*unlth(i,k-1) + nth2*unlth(i,k)
319 ELSE
320 var_reg(i,k) = nth1*unlth(i,k) + nth2*unlth(i
321 ENDIF
322 ENDDO
323 ENDDO
324 ENDIF
325
326
327
328
329
330 DO k = 1, nddl
331
332
333# include "vectorize.inc"
334 DO i = 1, nel
335
336
337 IF (off(i) /= zero) THEN
338 ! computation of
the product len**2 * btxb
339
340
341 b1 = (l2 / vol(i))*wf(k,nddl)*(btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1)
342 . - btb11(i)*unl(pos3(i)+k-1) - btb12(i)*unl(pos4(i)+k-1))
343
344 b2 = (l2 / vol(i))*wf(k,nddl)*(btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1)
345 . - btb12(i)*unl(pos3(i)+k-1) - btb22(i)*unl(pos4(i)+k-1))
346
347 b3 = (l2 / vol(i))*wf(k,nddl)*(-btb11(i)*unl(pos1(i)+k-1) - btb12(i)*unl(pos2(i)+k-1)
348 . + btb11(i)*unl(pos3(i)+k-1) + btb12(i)*unl(pos4(i)+k-1))
349
350 b4 = (l2 / vol(i))*wf(k,nddl)*(-btb12(i)*unl(pos1(i)+k-1) - btb22(i)*unl(pos2(i)+k-1)
351 . + btb12(i)*unl(pos3(i)+k-1) + btb22(i)*unl(pos4(i)+k-1))
352
353
354 ntn_unl = (unl(pos1(i)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) + unl(pos4(i)+k-1))/ntn
355
356
357 ntn_vnl = (vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) + vnl(pos4(i)+k-1))/ntn
358 IF (nodadt > 0) THEN
359 ntn_vnl
360 . sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1)),
361 . sqrt(mass(pos3(i
362 . sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1)))*ntn_vnl
363 ENDIF
364
365
366 ntn_unl = ntn_unl * vol(i) * wf(k,nddl)
367 ntn_vnl = ntn_vnl * xi * vol(i) * wf(k,nddl)
368
369
370 ntvar = var_reg(i,k)*fourth*vol(i)*wf(k,nddl)
371
372
373 a = ntn_unl + ntn_vnl - ntvar
374 f1(i,k) = a + b1
375 f2(i,k) = a + b2
376 f3(i,k) = a + b3
377 f4(i,k) = a + b4
378
379
380 IF (nodadt > 0) THEN
381 sti1(i,k) = wf(k,nddl)*(abs((l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
382 . abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
383 . abs(-(l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
384 . abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)))
385 sti2(i,k) = wf(k,nddl)*(abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
386 . abs((l2/vol(i))*btb22(i) + one/ntn*vol(i)) +
387 . abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
388 . abs(-(l2/vol(i))*btb22(i) + one
389 sti3(i,k) = wf(k,nddl)*(abs(-(l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
390 . abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
391 . abs((l2/vol(i))*btb11(i) + one/ntn*vol(i)) +
392 . abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)))
393 sti4(i,k) = wf(k,nddl)*(abs(-(l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
394 . abs(-(l2/vol(i))*btb22(i) + one/ntn*vol(i)) +
395 . abs((l2/vol(i))*btb12(i) + one/ntn*vol(i)) +
396 . abs((l2/vol(i))*btb22(i) + one/ntn*vol(i)))
397 ENDIF
398
399
400 ELSE
401 IF (nodadt > 0) THEN
402
403 f1(i,k) = wf(k,nddl)*sqrt(mass(pos1(i)+k-1)/mass0(pos1(i)+k-1))*zeta*sspnl*
404 . half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*sqrt(area0(i))*thk0(i)
405 f2(i,k) = wf(k,nddl)*sqrt(mass(pos2(i)+k-1)/mass0(pos2(i)+k-1))*zeta*sspnl*
406 . half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*sqrt(area0(i))*thk0(i)
407 f3(i,k) = wf(k,nddl)*sqrt
408 . half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*sqrt(area0
409 f4(i,k) = wf(k,nddl)*sqrt(mass(pos4(i)+k-1)/mass0(pos4(i)+k-1))*zeta*sspnl*
410 . half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*sqrt(area0(i))*thk0(i)
411
412 sti1(i,k) = em20
413 sti2(i,k) = em20
414 sti3(i,k) = em20
415 sti4(i,k) = em20
416 ELSE
417
418 f1(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos1(i)+k-1)+vnl0(pos1(i)+k-1))*sqrt(area0(i))*thk0(i)
419 f2(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos2(i)+k-1)+vnl0(pos2(i)+k-1))*sqrt(area0(i))*thk0(i)
420 f3(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos3(i)+k-1)+vnl0(pos3(i)+k-1))*sqrt(area0(i))*thk0(i)
421 f4(i,k) = wf(k,nddl)*zeta*sspnl*half*(vnl(pos4(i)+k-1)+vnl0(pos4(i)+k-1))*sqrt(area0(i))*thk0(i)
422 ENDIF
423 ENDIF
424 ENDDO
425 ENDDO
426
427
428
429
430
431
432 IF (iparit == 0) THEN
433
434 fnl => nloc_dmg%FNL(1:l_nloc,itask+1)
435 IF (nodadt > 0) stifnl => nloc_dmg%STIFNL(1:l_nloc,itask+1)
436
437 DO i=1,nel
438
439# include "vectorize.inc"
440 DO k = 1,nddl
441
442 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - f1(i,k)
443 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - f2(i,k)
444 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - f3(i,k)
445 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - f4(i,k)
446 IF (nodadt > 0) THEN
447
448 maxstif =
max(sti1(i,k),sti2(i,k),sti3(i,k),sti4(i,k))
449
450 stifnl(pos1(i)+k-1) = stifnl(pos1(i)+k-1) + maxstif
451 stifnl(pos2(i)+k-1) = stifnl(pos2(i)+k-1) + maxstif
452 stifnl(pos3(i)+k-1) = stifnl(pos3(i)+k-1) + maxstif
453 stifnl(pos4(i)+k-1) = stifnl(pos4(i)+k-1) + maxstif
454 ENDIF
455 ENDDO
456 ENDDO
457
458
459 ELSE
460
461 DO j = 1,nddl
462
463
464 DO i=1,nel
465 ii = i + nft
466
467
468 IF (nodadt > 0) THEN
469 maxstif =
max(sti1(i,j),sti2(i,j),sti3(i,j),sti4(i,j))
470 ENDIF
471
472 k = nloc_dmg%IADC(1,ii)
473 nloc_dmg%FSKY(k,j) = -f1(i,j)
474 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
475
476 k = nloc_dmg%IADC(2,ii)
477 nloc_dmg%FSKY(k,j) = -f2(i,j)
478 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
479
480 k = nloc_dmg%IADC(3,ii)
481 nloc_dmg%FSKY(k,j) = -f3(i,j)
482 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
483
484 k = nloc_dmg%IADC(4,ii)
485 nloc_dmg%FSKY(k,j) = -f4(i,j)
486 IF (nodadt > 0) nloc_dmg%STSKY(k,j) = maxstif
487
488 ENDDO
489
490 ENDDO
491
492 ENDIF
493
494
495
496
497 IF (nodadt == 0) THEN
498 DO i = 1,nel
499
500 IF (off(i) /= zero) THEN
501
502 dtnl = (two*(
min(le(i),le_max))*sqrt(three*zeta))/
503 . sqrt(twelve*l2 + (
min(le(i),le_max))**2)
504
505 IF (nddl>1) THEN
506 IF (nddl > 2) THEN
507 dtnl_th = (two*(
min(thk(i)/nddl,le_max))*sqrt(three*zeta))/
508 . sqrt(twelve*l2 + (
min(thk(i)/nddl,le_max))**2)
509 ELSE
510 dtnl_th = (two*(
min(thk(i),le_max))*sqrt(three*zeta))/
511 . sqrt(twelve*l2 + (
min(thk(i),le_max))**2)
512 ENDIF
513 ELSE
514 dtnl_th = ep20
515 ENDIF
516
517 dt2t =
min(dt2t,dtfac1(1)*cdamp*dtnl_th,dtfac1(1)*cdamp*dtnl)
518 ENDIF
519 ENDDO
520 ENDIF
521
522
523 IF (ALLOCATED(f1)) DEALLOCATE(f1)
524 IF (ALLOCATED(f2)) DEALLOCATE(f2)
525 IF (ALLOCATED(f3)) DEALLOCATE(f3)
526 IF (ALLOCATED(f4)) DEALLOCATE(f4)
527 IF (ALLOCATED(sti1)) DEALLOCATE(sti1)
528 IF (ALLOCATED(sti2)) DEALLOCATE(sti2)
529 IF (ALLOCATED(sti3)) DEALLOCATE(sti3)
530 IF (ALLOCATED(sti4)) DEALLOCATE(sti4)
531 IF (ALLOCATED(stifnlth)) DEALLOCATE(stifnlth)
532 IF (ALLOCATED(btb11)) DEALLOCATE(btb11)
533 IF (ALLOCATED(btb12)) DEALLOCATE(btb12)
534 IF (ALLOCATED(btb22)) DEALLOCATE(btb22)
535 IF (ALLOCATED(pos1)) DEALLOCATE(pos1)
536 IF (ALLOCATED(pos2)) DEALLOCATE(pos2)
537 IF (ALLOCATED(pos3)) DEALLOCATE(pos3)
538 IF (ALLOCATED(pos4)) DEALLOCATE(pos4)
539 IF (ALLOCATED(vol)) DEALLOCATE(vol)
540
end diagonal values have been computed in the(sparse) matrix id.SOL
subroutine area(d1, x, x2, y, y2, eint, stif0)