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