37
38
39
40 USE elbufdef_mod
43 use element_mod , only : nixc
44
45
46
47#include "implicit_f.inc"
48
49
50
51#include "param_c.inc"
52#include "com01_c.inc"
53#include "com04_c.inc"
54#include "scr03_c.inc"
55
56
57
58 TYPE(ELBUF_STRUCT_), TARGET, DIMENSION(NGROUP) :: ELBUF_TAB
59 TYPE (NLOCAL_STR_) , TARGET :: NLOC_DMG
60 INTEGER IXC(NIXC,*),NFT,NEL,NG,IPM(NPROPMI,*)
61 my_real ,
DIMENSION(NUMELC+NUMELTG),
INTENT(IN) ::
64 . x(3,*),xrefc(4,3,*),dt_nl,bufmat(*),time
65 LOGICAL :: FAILURE
66
67
68
69 INTEGER :: IMAT,NDOF,L_NLOC,N1,N2,N3,N4,
70 . K,I,NPTR,NPTS,IR,IS,NDNOD
71 INTEGER, DIMENSION(NEL) :: POS1,POS2,POS3,POS4
73 . le_min,len,damp,dens,ntn_unl,ntn_vnl,
74 . ntvar,z01(11,11),wf1(11,11),zn1(12,11),b1,b2,
75 . b3,b4,nth1,nth2,bth1,bth2,k1,k12,k2,sspnl,le_max
76 my_real,
DIMENSION(:,:),
ALLOCATABLE :: vpred,var_reg
78 . DIMENSION(:), POINTER :: fnl,unl,vnl,dnl,mnl,thck
79 my_real,
DIMENSION(NEL) :: x1,x2,x3,x4,
80 . y1,y2,y3,y4,px1,px2,py1,py2,e1x,e2x,e3x,
81 . e1y,e2y,e3y,e1z,e2z,e3z,x2l,y2l,x3l,y3l,
82 . x4l,y4l,z1,z2,z3,z4,surf,offg,vols,btb11,
83 . btb12,btb22
84 TYPE(BUF_NLOC_),POINTER :: BUFNL
85 my_real,
DIMENSION(:,:),
POINTER ::
86 . massth,fnlth,vnlth,unlth
87
88 DATA z01/
89 1 0. ,0. ,0. ,0. ,0. ,
90 1 0. ,0. ,0. ,0. ,0. ,0. ,
91 2 -.5 ,0.5 ,0. ,0. ,0. ,
92 2 0. ,0. ,0. ,0. ,0. ,0. ,
93 3 -.5 ,0. ,0.5 ,0. ,0. ,
94 3 0. ,0. ,0. ,0. ,0. ,0. ,
95 4 -.5 ,-.1666667,0.1666667,0.5 ,0. ,
96 4 0. ,0. ,0. ,0. ,0. ,0. ,
97 5 -.5 ,-.25 ,0. ,0.25 ,0.5 ,
98 5 0. ,0. ,0. ,0. ,0. ,0. ,
99 6 -.5 ,-.3 ,-.1 ,0.1 ,0.3 ,
100 6 0.5 ,0. ,0. ,0. ,0. ,0. ,
101 7 -.5 ,-.3333333,-.1666667,0.0 ,0.1666667,
102 7 0.3333333,0.5 ,0. ,0. ,0. ,0. ,
103 8 -.5 ,-.3571429,-.2142857,-.0714286,0.0714286,
104 8 0.2142857,0.3571429,0.5 ,0. ,0. ,0. ,
105 9 -.5 ,-.375 ,-.25 ,-.125 ,0.0 ,
106 9 0.125 ,0.25 ,0.375 ,0.5 ,0. ,0. ,
107 a -.5 ,-.3888889,-.2777778,-.1666667,-.0555555,
108 a 0.0555555,0.1666667,0.2777778,0.3888889,0.5 ,0. ,
109 b -.5 ,-.4 ,-.3 ,-.2 ,-.1 ,
110 b 0. ,0.1 ,0.2 ,0.3 ,0.4 ,0.5 /
111
112 DATA wf1/
113 1 1. ,0. ,0. ,0. ,0. ,
114 1 0. ,0. ,0. ,0. ,0. ,0. ,
115 2 0.5 ,0.5 ,0. ,0. ,0. ,
116 2 0. ,0. ,0. ,0. ,0. ,0. ,
117 3 0.25 ,0.5 ,0.25 ,0. ,0. ,
118 3 0. ,0. ,0. ,0. ,0. ,0. ,
119 4 0.1666667,0.3333333,0.3333333,0.1666667,0. ,
120 4 0. ,0. ,0. ,0. ,0. ,0. ,
121 5 0.125 ,0.25 ,0.25 ,0.25 ,0.125 ,
122 5 0. ,0. ,0. ,0. ,0. ,0. ,
123 6 0.1 ,0.2 ,0.2 ,0.2 ,0.2 ,
124 6 0.1 ,0. ,0. ,0. ,0. ,0. ,
125 7 0.0833333,0.1666667,0.1666667,0.1666667,0.1666667,
126 7 0.1666667,0.0833333,0. ,0. ,0. ,0. ,
127 8 0.0714286,0.1428571,0.1428571,0.1428571,0.1428571,
128 8 0.1428571,0.1428571,0.0714286,0. ,0. ,0. ,
129 9 0.0625 ,0.125 ,0.125 ,0.125 ,0.125 ,
130 9 0.125 ,0.125 ,0.125 ,0.0625 ,0. ,0. ,
131 a 0.0555556,0.1111111,0.1111111,0.1111111,0.1111111,
132 a 0.1111111,0.1111111,0.1111111,0.1111111,0.0555556,0. ,
133 b 0.05 ,0.1 ,0.1 ,0.1 ,0.1 ,
134 b 0.1 ,0.1 ,0.1 ,0.1 ,0.1 ,0.05 /
135
136 DATA zn1/
137 1 0. ,0. ,0. ,0. ,0. ,0. ,
138 1 0. ,0. ,0. ,0. ,0. ,0. ,
139 2 -.5 ,0.5 ,0. ,0. ,0. ,0. ,
140 2 0. ,0. ,0. ,0. ,0. ,0. ,
141 3 -.5 ,-.25 ,0.25 ,0.5 ,0. ,0. ,
142 3 0. ,0. ,0. ,0. ,0. ,0. ,
143 4 -.5 ,-.3333333,0. ,0.3333333,0.5 ,0. ,
144 4 0. ,0. ,0. ,0. ,0. ,0. ,
145 5 -.5 ,-.375 ,-0.125 ,0.125 ,0.375 ,0.5 ,
146 5 0. ,0. ,0. ,0. ,0. ,0. ,
147 6 -.5 ,-.4 ,-.2 ,0.0 ,0.2 ,0.4 ,
148 6 0.5 ,0. ,0. ,0. ,0. ,0. ,
149 7 -.5 ,-.4166667,-.25 ,-.0833333,0.0833333,0.25 ,
150 7 0.4166667,0.5 ,0. ,0. ,0. ,0. ,
151 8 -.5 ,-.4285715,-.2857143,-.1428572,0.0 ,0.1428572,
152 8 0.2857143,0.4285715,0.5 ,0. ,0. ,0. ,
153 9 -.5 ,-.4375 ,-.3125 ,-.1875 ,-.0625 ,0.0625 ,
154 9 0.1875 ,0.3125 ,0.4375 ,0.5 ,0. ,0. ,
155 a -.5 ,-.4444444,-.3333333,-.2222222,-.1111111,0. ,
156 a 0.1111111,0.2222222,0.3333333,0.4444444,0.5 ,0. ,
157 b -.5 ,-.45 ,-.35 ,-.25 ,-.15 ,-.05 ,
158 b 0.05 ,0.15 ,0.25 ,0.35 ,0.45 ,0.5 /
159
160
161 l_nloc = nloc_dmg%L_NLOC
162
163 fnl => nloc_dmg%FNL(1:l_nloc,1)
164 vnl => nloc_dmg%VNL(1:l_nloc)
165 dnl => nloc_dmg%DNL(1:l_nloc)
166 unl => nloc_dmg%UNL(1:l_nloc)
167 mnl => nloc_dmg%MASS(1:l_nloc)
168
169 imat = ixc(1,1+nft)
170
171 le_min = sqrt(minval(
area(nft+1:nft+nel)))
172
173 ndof = elbuf_tab(ng)%BUFLY(1)%NPTT
174
175 thck => elbuf_tab(ng)%GBUF%THK(1:nel)
176
177 IF (ndof>1) THEN
178 le_min =
min(le_min,minval(thck(1:nel)
179 ENDIF
180
181 len = nloc_dmg%LEN(imat)
182
183 le_max = nloc_dmg%LE_MAX(imat)
184
185 damp = nloc_dmg%DAMP(imat)
186
187 dens = nloc_dmg%DENS(imat)
188
189 sspnl = nloc_dmg%SSPNL(imat)
190
191 dt_nl =
min(dt_nl,0.5d0*((two*
min(le_min
192 . (sqrt(twelve*(len**2)+(
min(le_min,le_max)**2)))))
193
194 IF (ndof>1) THEN
195 IF (ndof > 2) THEN
196 ALLOCATE(vpred(nel,ndof+1))
197 ndnod = ndof+1
198 ELSE
199 ALLOCATE(vpred(nel,ndof))
200 ndnod = ndof
201 ENDIF
202 ENDIF
203
204 nptr = elbuf_tab(ng)%NPTR
205 npts = elbuf_tab(ng)%NPTS
206
207 IF (.NOT.ALLOCATED(var_reg)) ALLOCATE(var_reg(nel,ndof))
208
209
210# include "vectorize.inc"
211 DO i = 1,nel
212
213 IF (nxref == 0) THEN
214 x1(i)=x(1,ixc(2,nft+i))
215 y1(i)=x(2,ixc(2,nft+i))
216 z1(i)=x(3,ixc(2,nft+i))
217 x2(i)=x(1,ixc(3,nft+i))
218 y2(i)=x(2,ixc(3,nft+i))
219 z2(i)=x(3,ixc(3,nft+i))
220 x3(i)=x(1,ixc(4,nft+i))
221 y3(i)=x(2,ixc(4,nft+i))
222 z3(i)=x(3,ixc(4,nft+i))
223 x4(i)=x(1,ixc(5,nft+i))
224 y4(i)=x(2,ixc(5,nft+i))
225 z4(i)=x(3,ixc(5,nft+i))
226 ELSE
227 x1(i)=xrefc(1,1,nft+i)
228 y1(i)=xrefc(1,2,nft+i)
229 z1(i)=xrefc(1,3,nft+i)
230 x2(i)=xrefc(2,1,nft+i)
231 y2(i)=xrefc(2,2,nft+i)
232 z2(i)=xrefc(2,3,nft+i)
233 x3(i)=xrefc(3,1,nft+i)
234 y3(i)=xrefc(3,2,nft+i)
235 z3(i)=xrefc(3,3,nft+i)
236 x4(i)=xrefc(4,1,nft+i)
237 y4(i)=xrefc(4,2,nft+i)
238 z4(i)=xrefc(4,3,nft+i)
239 ENDIF
240
241
242 n1 = nloc_dmg%IDXI(ixc(2,nft+i))
243 n2 = nloc_dmg%IDXI(ixc(3,nft+i))
244 n3 = nloc_dmg%IDXI(ixc(4,nft+i))
245 n4 = nloc_dmg%IDXI(ixc(5,nft+i))
246
247 pos1(i) = nloc_dmg%POSI(n1)
248 pos2(i) = nloc_dmg%POSI(n2)
249 pos3(i) = nloc_dmg%POSI(n3)
250 pos4(i) = nloc_dmg%POSI(n4)
251 ENDDO
252
253
254
255 DO k = 1,ndof
256
257# include "vectorize.inc"
258 DO i = 1,nel
259 var_reg(i,k) = fourth*(dnl(pos1(i)+k-1) + dnl(pos2(i)+k-1)
260 . + dnl(pos3(i)+k-1) + dnl(pos4(i)+k-1))
261 ENDDO
262 ENDDO
263
265 . x1 ,x2 ,x3 ,x4 ,y1 ,y2 ,
266 . y3 ,y4 ,z1 ,z2 ,z3 ,z4 ,
267 . e1x ,e2x ,e3x ,e1y ,e2y ,e3y ,
268 . e1z ,e2z ,e3z )
269
270
272 . bufmat ,time ,var_reg ,
273 . failure )
274
275
276
277
278
279
280# include "vectorize.inc"
281 DO i=1,nel
282
283
284 x2l(i)=e1x(i)*(x2(i)-x1(i))+e1y(i)*(y2(i)-y1(i))+e1z(i)*(z2(i)-z1(i))
285 y2l(i)=e2x(i)*(x2(i)-x1(i))+e2y(i)*(y2(i)-y1(i))+e2z(i)*(z2(i)-z1(i))
286 x3l(i)=e1x(i)*(x3(i)-x1(i))+e1y(i)*(y3(i)-y1(i))+e1z(i)*(z3(i)-z1(i))
287 y3l(i)=e2x(i)*(x3(i)-x1(i))+e2y(i)*(y3(i)-y1(i))+e2z(i)*(z3(i)-z1(i))
288 x4l(i)=e1x(i)*(x4(i)-x1(i))+e1y(i)*(y4(i)-y1(i))+e1z(i)*(z4(i)-z1(i)
289 y4l(i)=e2x(i)*(x4(i)-x1(i))+e2y(i)*(y4(i)-y1(i))+e2z(i)*(z4(i)-z1(i))
290
291 px1(i) = half *(y2l(i)-y4l(i))
292 px2(i) = half * y3l(i)
293 py1(i) = -half *(x2l(i)-x4l(i))
294 py2(i) = -half * x3l(i)
295
296
297 btb11(i) = px1(i)**2 + py1(i)**2
298 btb12(i) = px1(i)*px2(i) + py1(i)*py2(i)
299 btb22(i) = px2(i)**2 + py2(i)**2
300
301
302 vols(i) =
area(nft+i)*thck(i)
303
304
305 offg(i) = elbuf_tab(ng)%GBUF%OFF(i)
306
307 ENDDO
308
309
310
311
312
313 IF ((ndof > 1).AND.(len>zero)) THEN
314
315
316 bufnl => elbuf_tab(ng)%NLOC(1,1)
317
318
319 massth => bufnl%MASSTH(1:nel,1:ndnod)
320 unlth => bufnl%UNLTH(1:nel ,1:ndnod)
321 vnlth => bufnl%VNLTH(1:nel ,1:ndnod)
322 fnlth => bufnl%FNLTH(1:nel ,1:ndnod)
323
324 DO k = 1,ndnod
325 DO i = 1,nel
326
327 vpred(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*(dt_nl/two)
328 ENDDO
329 ENDDO
330 DO k = 1,ndnod
331 DO i = 1,nel
332
333 fnlth(i,k) = zero
334 ENDDO
335 ENDDO
336
337
338 DO k = 1, ndof
339
340
341 IF ((ndof==2).AND.(k==2)) THEN
342 nth1 = (z01(k,ndof) - zn1(k,ndof))/(zn1(k-1,ndof) - zn1(k,ndof))
343 nth2 = (z01(k,ndof) - zn1(k-1,ndof)) /(zn1(k,ndof) - zn1(k-1,ndof
344 ELSE
345 nth1 = (z01(k,ndof) - zn1(k+1,ndof))/(zn1(k,ndof) - zn1(k+1,ndof))
346 nth2 = (z01(k,ndof) - zn1(k,ndof)) /(zn1(k+1,ndof) - zn1(k,ndof))
347 ENDIF
348
349
350 DO i = 1,nel
351
352 IF ((ndof==2).AND.(k==2)) THEN
353 bth1 = (one/(zn1(k-1,ndof) - zn1(k,ndof)))*(one/thck(i))
354 bth2 = (one/(zn1(k,ndof) - zn1(k-1,ndof)))*(one/thck(i))
355 ELSE
356 bth1 = (one/(zn1(k,ndof) - zn1(k+1,ndof)))*(one/thck(i))
357 bth2 = (one/(zn1(k+1,ndof) - zn1(k,ndof)))*(one/thck(i))
358 ENDIF
359
360
361 k1 = (len**2)*(bth1**2) + nth1**2
362 k12 = (len**2)*(bth1*bth2)+ (nth1*nth2)
363 k2 = (len**2)*(bth2**2) + nth2**2
364
365
366 IF ((ndof==2).AND.(k==2)) THEN
367 fnlth(i,k-1) = fnlth(i,k-1) + (k1*unlth(i,k-1) + k12*unlth(i,k)
368 . + damp*((nth1**2)*vpred(i,k-1)
369 . + (nth1*nth2)*vpred(i,k))
370 . - (nth1*var_reg(i,k)))*vols(i)*wf1(k,ndof)
371 fnlth(i,k) = fnlth(i,k) + (k12*unlth(i,k-1) + k2*unlth(i,k)
372 . + damp*(nth1*nth2*vpred(i,k-1)
373 . + (nth2**2)*vpred(i,k))
374 . - nth2*var_reg(i,k))*vols(i)*wf1(k,ndof)
375 ELSE
376 fnlth(i,k) = fnlth(i,k) + (k1*unlth(i,k) + k12*unlth(i,k+1)
377 . + damp*((nth1**2)*vpred
378 . + (nth1*nth2)*vpred(i,k+1))
379 . - (nth1*var_reg(i,k)))*vols(i)*wf1(k,ndof)
380 fnlth(i,k+1) = fnlth(i,k+1) + (k12*unlth(i
381 . + damp*(nth1
382 . + (nth2**2)*vpred(i,k+1))
383 . - nth2*var_reg(i,k))*vols(i)*wf1(k,ndof)
384 ENDIF
385 ENDDO
386 ENDDO
387
388 DO k = 1,ndnod
389 DO i = 1,nel
390
391 vnlth(i,k) = vnlth(i,k) - (fnlth(i,k)/massth(i,k))*dt_nl
392 ENDDO
393 ENDDO
394
395 DO k = 1,ndnod
396 DO i = 1,nel
397
398 unlth(i,k) = unlth(i,k) + vnlth(i,k)*dt_nl
399 ENDDO
400 ENDDO
401
402
403 DO k = 1, ndof
404
405 IF ((ndof==2).AND.(k==2)) THEN
406 nth1 = (z01(k,ndof) - zn1(k,ndof))/(zn1(k-1,ndof) - zn1(k,ndof))
407 nth2 = (z01(k,ndof) - zn1(k-1,ndof)) /(zn1(k,ndof) - zn1(k-1,ndof))
408 ELSE
409 nth1 = (z01(k,ndof) - zn1(k+1,ndof))/(zn1(k,ndof) - zn1(k+1,ndof))
410 nth2 = (z01(k,ndof) - zn1(k,ndof)) /(zn1(k+1,ndof) - zn1(k,ndof))
411 ENDIF
412
413 DO i = 1,nel
414
415 IF ((ndof==2).AND.(k==2)) THEN
416 var_reg(i,k) = nth1*unlth(i,k-1) + nth2*unlth(i,k)
417 ELSE
418 var_reg(i,k) = nth1*unlth(i,k) + nth2*unlth(i,k+1)
419 ENDIF
420 ENDDO
421 ENDDO
422
423 IF ((nptr>1).OR.(npts>1)) THEN
424 DO ir = 1,nptr
425 DO is = 1,npts
426 bufnl => elbuf_tab(ng)%NLOC(ir,is)
427 bufnl%MASSTH(1:nel,1:ndnod) = elbuf_tab(ng)%NLOC(1,1)%MASSTH(1:nel,1:ndnod)
428 bufnl%MASSTH(1:nel,1:ndnod) = elbuf_tab(ng)%NLOC(1,1)%UNLTH(1:nel,1:ndnod)
429 bufnl%MASSTH(1:nel,1:ndnod) = elbuf_tab(ng)%NLOC(1,1)%VNLTH(1:nel,1:ndnod)
430 bufnl%MASSTH(1:nel,1:ndnod) = elbuf_tab(ng)%NLOC(1,1)%FNLTH(1:nel,1:ndnod)
431 ENDDO
432 ENDDO
433 ENDIF
434 ENDIF
435
436
437
438
439
440 DO k = 1,ndof
441
442
443# include "vectorize.inc"
444 DO i = 1,nel
445
446
447 IF (offg(i) > zero) THEN
448
449 b1 = ((len**2)/vols(i))*wf1(k,ndof)*(btb11(i)*unl(pos1(i)+k-1) + btb12(i)*unl(pos2(i)+k-1)
450 . - btb11(i)*unl(pos3(i)+k-1) - btb12(i)*unl(pos4(i)+k-1))
451 b2 = ((len**2)/vols(i))*wf1(k,ndof)*(btb12(i)*unl(pos1(i)+k-1) + btb22(i)*unl(pos2(i)+k-1)
452 . - btb12(i)*unl(pos3(i)+k-1) - btb22(i)*unl(pos4(i)+k-1))
453 b3 = ((len**2)/vols(i))*wf1(k,ndof)*(-btb11(i)*unl(pos1(i)+k-1) - btb12(i)*unl(pos2(i)+k-1)
454 . + btb11(i)*unl(pos3(i)+k-1) + btb12(i)*unl(pos4(i)+k-1))
455 b4 = ((len**2)/vols(i))*wf1(k,ndof)*(-btb12(i)*unl(pos1(i)+k-1) - btb22(i)*unl(pos2(i)+k-1)
456 . + btb12(i)*unl(pos3(i)+k-1) + btb22(i)*unl(pos4(i)+k-1))
457
458 ntn_unl = ((unl(pos1(i)+k-1) + unl(pos2(i)+k-1) + unl(pos3(i)+k-1) +
459 . unl(pos4(i)+k-1))*fourth*fourth)*vols(i)*wf1(k,ndof)
460
461 ntn_vnl = ((vnl(pos1(i)+k-1) + vnl(pos2(i)+k-1) + vnl(pos3(i)+k-1) +
462 . vnl(pos4(i)+k-1))*fourth*fourth)*damp*vols(i)*wf1(k,ndof)
463
464 ntvar = var_reg(i,k)*fourth*vols(i)*wf1(k,ndof)
465
466 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - (ntn_unl + ntn_vnl - ntvar + b1)
467 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - (ntn_unl + ntn_vnl - ntvar + b2)
468 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - (ntn_unl + ntn_vnl - ntvar + b3)
469 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - (ntn_unl + ntn_vnl - ntvar + b4)
470
471 ELSE
472
473 fnl(pos1(i)+k-1) = fnl(pos1(i)+k-1) - wf1(k,ndof)*dens*sspnl*vnl(pos1(i)+k-1)*le_max*thck(i)
474 fnl(pos2(i)+k-1) = fnl(pos2(i)+k-1) - wf1(k,ndof)*dens*sspnl*vnl(pos2(i)+k-1)*le_max*thck(i)
475 fnl(pos3(i)+k-1) = fnl(pos3(i)+k-1) - wf1(k,ndof)*dens*sspnl*vnl(pos3(i)+k-1)*le_max*thck(i)
476 fnl(pos4(i)+k-1) = fnl(pos4(i)+k-1) - wf1(k,ndof)*dens*sspnl*vnl(pos4(i)+k-1)*le_max*thck(i)
477 ENDIF
478 ENDDO
479 ENDDO
480
481 IF (ALLOCATED(var_reg)) DEALLOCATE(var_reg)
482 IF (ALLOCATED(vpred)) DEALLOCATE(vpred)
subroutine cneveci(jft, jlt, area, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, e1x, e2x, e3x, e1y, e2y, e3y, e1z, e2z, e3z)
subroutine cnloc_matini(elbuf_str, nel, ipm, bufmat, time, varnl, failure)
subroutine area(d1, x, x2, y, y2, eint, stif0)