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