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