36
37
38
39 USE elbufdef_mod
42 use element_mod , only : nixs
43
44
45
46#include "implicit_f.inc"
47
48
49
50#include "param_c.inc"
51#include "units_c.inc"
52#include "comlock.inc"
53#include "com01_c.inc"
54#include "com04_c.inc"
55#include "com08_c.inc"
56#include "scr14_c.inc"
57#include "scr17_c.inc"
58
59
60
61 INTEGER IXS(NIXS,*),IPARG(NPARG,*)
62 my_real ,
DIMENSION(3,*) :: x,a,ar,fcluster,mcluster
63 my_real ,
DIMENSION(LSKEW,*) :: skew
64 TYPE (CLUSTER_) ,DIMENSION(NCLUSTER) :: CLUSTER
65 TYPE (ELBUF_STRUCT_),DIMENSION(NGROUP) :: ELBUF_TAB
66 TYPE (H3D_DATABASE) :: H3D_DATA
67
68
69
70 INTEGER I,J,K,IL,IEL,NG,NFT,NNOD,ISKN,N,N1,N2,N3,N4,NINDX,IFAIL,IPID
71 INTEGER CLUSTERNOD(NCLUSTER),LCLUSTER(NCLUSTER),LCL(NCLUSTER),
72 . IAD(NCLUSTER)
73 INTEGER INDX(NCLUSTER)
74 my_real,
DIMENSION(NPROPG,*) :: geo
75 my_real,
DIMENSION(3) :: fbot,ftop,mbot,mtop,m1,xg,x1,x2
76 my_real,
DIMENSION(3,NCLUSTER) :: vn,vx,vy
77 my_real :: fn,ft,mr,mb,dmg,xm,ym,zm,dx1,dy1,dz1,dx2,dy2,dz2,
78 . fx,fy,fz,momx,momy,momz,
norm,critf,critm,drx,dry,drz,
79 . sx,sy,sz,tx,ty,tz
80 my_real,
DIMENSION(NCLUSTER) :: tthick
81
82 nindx = 0
83 tthick(1:ncluster) = zero
84
85 DO i = 1, ncluster
86 IF (cluster(i)%OFF == 0) cycle
87 nnod = cluster(i)%NNOD
88 iskn = cluster(i)%SKEW
89 ifail= cluster(i)%IFAIL
90
91
92
93 x1(1:3) = zero
94 DO j = 1,nnod
95 n1 = cluster(i)%NOD1(j)
96 x1(1) = x1(1) + x(1,n1)
97 x1(2) = x1(2) + x(2,n1)
98 x1(3) = x1(3) + x(3,n1)
99 END DO
100 xm = x1(1) / nnod
101 ym = x1(2) / nnod
102 zm = x1(3) / nnod
103
104
105
106 IF (ifail > 0 .and. iskn == 0) THEN
107
108
109 vn(1,i) = zero
110 vn(2,i) = zero
111 vn(3,i) = zero
112
113 IF (cluster(i)%TYPE == 1) THEN
114 DO j = 1,cluster(i)%NEL
115 ng = cluster(i)%NG(j)
116 iel = cluster(i)%ELEM(j)
117 nft = iparg(3,ng)
118 ipid = ixs(10,nft+iel)
119 n1 = ixs(2,nft+iel)
120 n2 = ixs(3,nft+iel)
121 n3 = ixs(4,nft+iel)
122 n4 = ixs(5,nft+iel)
123 tthick(i) = geo(41,ipid)
124 sx = x(1,n3) - x(1,n1)
125 sy = x(2,n3) - x(2,n1)
126 sz = x(3,n3) - x(3,n1)
127 tx = x(1,n4) - x(1,n2)
128 ty = x(2,n4) - x(2,n2)
129 tz = x(3,n4) - x(3,n2)
130 vn(1,i) = vn(1,i) + sy*tz - sz*ty
131 vn(2,i) = vn(2,i) + sz*tx - sx*tz
132 vn(3,i) = vn(3,i) + sx*ty - sy*tx
133 END DO
134
135 ELSE
136 n1 = cluster(i)%NOD1(nnod)
137 n2 = cluster(i)%NOD1(1)
138 sx = xm - x(1,n1)
139 sy = ym - x(2,n1)
140 sz = zm - x(3,n1)
141 tx = xm - x(1,n2)
142 ty = ym - x(2,n2)
143 tz = zm - x(3,n2)
144 vn(1,i) = vn(1,i) + sy*tz - sz*ty
145 vn(2,i) = vn(2,i) + sz*tx - sx*tz
146 vn(3,i) = vn(3,i) + sx*ty - sy*tx
147 DO j = 1,nnod-1
148 n1 = cluster(i)%NOD1(j)
149 n2 = cluster(i)%NOD1(j+1)
150 sx = xm - x(1,n1)
151 sy = ym - x(2,n1)
152 sz = zm - x(3,n1)
153 tx = xm - x(1,n2)
154 ty = ym - x(2,n2)
155 tz = zm - x(3,n2)
156 vn(1,i) = vn(1,i) + sy*tz - sz*ty
157 vn(2,i) = vn(2,i) + sz*tx - sx*tz
158 vn(3,i) = vn(3,i) + sx*ty - sy*tx
159 END DO
160 END IF
161
162 norm = one / sqrt(vn(1,i)**2 + vn(2,i)**2 + vn(3,i)**2)
163 vn(1,i) = vn(1,i)*
norm
164 vn(2,i) = vn(2,i)*
norm
165 vn(3,i) = vn(3,i)*
norm
166
167
168
169 n1 = cluster(i)%NOD1(1)
170 n2 = cluster(i)%NOD1(2)
171 vx(1,i) = x(1,n1) - xm
172 vx(2,i) = x(2,n1) - ym
173 vx(3,i) = x(3,n1) - zm
174 vy(1,i) = vn(2,i)*vx(3,i) - vn(3,i)*vx(2,i)
175 vy(2,i) = vn(3,i)*vx(1,i) - vn(1,i)*vx(3,i)
176 vy(3,i) = vn(1,i)*vx(2,i) - vn(2,i)*vx(1,i)
177 norm = one / sqrt(vy(1,i)**2 + vy(2,i)**2 + vy(3,i)**2)
178 vy(1,i) = vy(1,i)*
norm
179 vy(2,i) = vy(2,i)*
norm
180 vy(3,i) = vy(3,i)*
norm
181 vx(1,i) = vy(2,i)*vn(3,i) - vy(3,i)*vn(2,i)
182 vx(2,i) = vy(3,i)*vn(1,i) - vy(1,i)*vn(3,i)
183 vx(3,i) = vy(1,i)*vn(2,i) - vy(2,i)*vn(1,i)
184 norm = one / sqrt(vx(1,i)**2 + vx(2,i)**2 + vx(3,i)**2)
185 vx(1,i) = vx(1,i)*
norm
186 vx(2,i) = vx(2,i)*
norm
187 vx(3,i) = vx(3,i)*
norm
188
189 ENDIF
190
191
192
193 fbot = zero
194 ftop = zero
195 DO j = 1,nnod
196 n1 = cluster(i)%NOD1(j)
197 n2 = cluster(i)%NOD2(j)
198 fbot(1) = fbot(1) + a(1,n1)
199 fbot(2) = fbot(2) + a(2,n1)
200 fbot(3) = fbot(3) + a(3,n1)
201 ftop(1) = ftop(1) + a(1,n2)
202 ftop(2) = ftop(2) + a(2,n2)
203 ftop(3) = ftop(3) + a(3,n2)
204 END DO
205
206
207
208 mbot = zero
209 mtop = zero
210
211 IF (cluster(i)%TYPE == 1 .and. iskn == 0 .and. tthick(i) > zero) THEN
212 DO j = 1,nnod
213 n1 = cluster(i)%NOD1(j)
214 n2 = cluster(i)%NOD2(j)
215
216 drx = x(1,n2) - xm
217 dry = x(2,n2) - ym
218 drz = sign(tthick(i), x(3,n2) - zm)
219 mtop(1) = mtop(1) + dry*a(3,n2) - drz*a(2,n2)
220 mtop(2) = mtop(2) + drz*a(1,n2) - drx*a(3,n2)
221 mtop(3) = mtop(3) + drx*a(2,n2) - dry*a(1,n2)
222
223 drx = x(1,n1) - xm
224 dry = x(2,n1) - ym
225 mbot(1) = mbot(1) + dry*a(3,n1)
226 mbot(2) = mbot(2) - drx*a(3,n1)
227 mbot(3) = mbot(3) + drx*a(2,n1) - dry*a(1,n1)
228 END DO
229 ELSE
230 DO j = 1,nnod
231 n1 = cluster(i)%NOD1(j)
232 n2 = cluster(i)%NOD2(j)
233
234 drx = x(1,n2) - xm
235 dry = x(2,n2) - ym
236 drz = x(3,n2) - zm
237 mtop(1) = mtop(1) + dry*a(3,n2) - drz*a(2,n2)
238 mtop(2) = mtop(2) + drz*a(1,n2) - drx*a(3,n2)
239 mtop(3) = mtop(3) + drx*a(2,n2) - dry*a(1,n2)
240
241 drx = x(1,n1) - xm
242 dry = x(2,n1) - ym
243 drz = x(3,n1) - zm
244 mbot(1) = mbot(1) + dry*a(3,n1) - drz*a(2,n1)
245 mbot(2) = mbot(2) + drz*a(1,n1) - drx*a(3,n1)
246 mbot(3) = mbot(3) + drx*a(2,n1) - dry*a(1,n1)
247 END DO
248 END IF
249
250 IF (cluster(i)%TYPE == 1) THEN
251 fx = (ftop(1) - fbot(1))*half
252 fy = (ftop(2) - fbot(2))*half
253 fz = (ftop(3) - fbot(3))*half
254 momx = (mtop(1) - mbot(1))*half
255 momy = (mtop(2) - mbot(2))*half
256 momz = (mtop(3) - mbot(3))*half
257 ELSE
258 fx = ftop(1)
259 fy = ftop(2)
260 fz = ftop(3)
261 momx = mtop(1)
262 momy = mtop(2)
263 momz = mtop(3)
264 DO j = 1,nnod
265 n1 = cluster(i)%NOD1(j)
266 n2 = cluster(i)%NOD2(j)
267 momx = momx + ar(1,n2)
268 momy = momy + ar(2,n2)
269 momz = momz + ar(3,n2)
270 END DO
271 ENDIF
272
273 cluster(i)%FOR(1) = fx
274 cluster(i)%FOR(2) = fy
275 cluster(i)%FOR(3) = fz
276 cluster(i)%MOM(1) = momx
277 cluster(i)%MOM(2) = momy
278 cluster(i)%MOM(3) = momz
279
280 ENDDO
281
282
283
284
285 nindx = 0
286 indx = 0
287
288 DO i = 1, ncluster
289
290 IF (cluster(i)%OFF == 0) THEN
291 cluster(i)%FOR(1) = zero
292 cluster(i)%FOR(2) = zero
293 cluster(i)%FOR(3) = zero
294 cluster(i)%MOM(1) = zero
295 cluster(i)%MOM(2) = zero
296 cluster(i)%MOM(3) = zero
297 END IF
298
299 nnod = cluster(i)%NNOD
300 iskn = cluster(i)%SKEW
301 ifail= cluster(i)%IFAIL
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324 IF (ifail > 0) THEN
325 IF (iskn > 0) THEN
326 fbot(1) = cluster(i)%FOR(1)*skew(1,iskn) +
327 . cluster(i)%FOR(2)*skew(2,iskn) +
328 . cluster(i)%FOR(3)*skew(3,iskn)
329 fbot(2) = cluster(i)%FOR(1)*skew(4,iskn) +
330 . cluster(i)%FOR(2)*skew(5,iskn) +
331 . cluster(i)%FOR(3)*skew(6,iskn)
332 fbot(3) = cluster(i)%FOR(1)*skew(7,iskn) +
333 . cluster(i)%FOR(2)*skew(8,iskn) +
334 . cluster(i)%FOR(3)*skew(9,iskn)
335 m1(1) = cluster(i)%MOM(1)*skew(1,iskn) +
336 . cluster(i)%MOM(2)*skew(2,iskn) +
337 . cluster(i)%MOM(3)*skew(3,iskn)
338 m1(2) = cluster(i)%MOM(1)*skew(4,iskn) +
339 . cluster(i)%MOM(2)*skew(5,iskn) +
340 . cluster(i)%MOM(3)*skew(6,iskn)
341 m1(3) = cluster(i)%MOM(1)*skew(7,iskn) +
342 . cluster(i)%MOM(2)*skew(8,iskn) +
343 . cluster(i)%MOM(3)*skew(9,iskn)
344 ELSE
345 fbot(1) = cluster(i)%FOR(1)*vx(1,i) +
346 . cluster(i)%FOR(2)*vx(2,i) +
347 . cluster(i)%FOR(3)*vx(3,i)
348 fbot(2) = cluster(i)%FOR(1)*vy(1,i) +
349 . cluster(i)%FOR(2)*vy(2,i) +
350 . cluster(i)%FOR(3)*vy(3,i)
351 fbot(3) = cluster(i)%FOR(1)*vn(1,i) +
352 . cluster(i)%FOR(2)*vn(2,i) +
353 . cluster(i)%FOR(3)*vn(3,i)
354 m1(1) = cluster(i)%MOM(1)*vx(1,i) +
355 . cluster(i)%MOM(2)*vx(2,i) +
356 . cluster(i)%MOM(3)*vx(3,i)
357 m1(2) = cluster(i)%MOM(1)*vy(1,i) +
358 . cluster(i)%MOM(2)*vy(2,i) +
359 . cluster(i)%MOM(3)*vy(3,i)
360 m1(3) = cluster(i)%MOM(1)*vn(1,i) +
361 . cluster(i)%MOM(2)*vn(2,i) +
362 . cluster(i)%MOM(3)*vn(3,i)
363
364
365 ENDIF
366 fn = abs(fbot(3))
367 ft = sqrt(fbot(1)*fbot(1) + fbot(2)*fbot(2))
368 mr = abs(m1(3))
369 mb = sqrt(m1(1)*m1(1) + m1(2)*m1(2))
370 ENDIF
371
372
373 IF (ifail == 1) THEN
374
375 critf =
max(fn/cluster(i)%FMAX(1),ft/cluster(i)%FMAX(2))
376 critm =
max(mr/cluster(i)%MMAX(1),mb/cluster
377 dmg =
max(critf,critm)
378
379 ELSEIF (ifail == 2) THEN
380
381 dmg = fourth*(
min(one+em10, fn/cluster(i)%FMAX(1)) +
382 .
min(one+em10, ft/cluster(i)%FMAX(2)) +
383 .
min(one+em10, mr/cluster(i)%MMAX(1)) +
384 .
min(one+em10, mb/cluster(i)%MMAX(2)))
385
386 ELSEIF (ifail == 3) THEN
387
388 dmg =
389 . cluster(i)%AX(1)*(fn/cluster(i)%FMAX(1))**cluster(i)%NX(1)
390 . + cluster(i)%AX(2)*(ft/cluster(i)%FMAX(2))**cluster(i)%NX(2)
391 . + cluster(i)%AX(3)*(mr/cluster(i)%MMAX(1))**cluster
392 . + cluster(i)%AX(4)*(mb/cluster(i)%MMAX(2))**cluster(i)%NX(4)
393 ELSE
394 dmg = zero
395 ENDIF
396 cluster(i)%FAIL = dmg
397
398 IF (dmg > one) THEN
399
400 nindx = nindx+1
401 indx(nindx) = i
402 idel7nok = 1
403 cluster(i)%OFF = 0
404 cluster(i)%FOR(1) = zero
405 cluster(i)%FOR(2) = zero
406 cluster(i)%FOR(3) = zero
407 cluster(i)%MOM(1) = zero
408 cluster(i)%MOM(2) = zero
409 cluster(i)%MOM(3) = zero
410 IF (cluster(i)%TYPE == 1) THEN
411 DO j = 1,cluster(i)%NEL
412 ng = cluster(i)%NG(j)
413 iel = cluster(i)%ELEM(j)
414 elbuf_tab(ng)%GBUF%OFF(iel) = zero
415 END DO
416 ELSE
417
418 ENDIF
419 ENDIF
420
421 ENDDO
422
423 IF (anim_v(19) + h3d_data%N_VECT_CLUST_FORCE > 0) THEN
424 DO i = 1, ncluster
425 nnod = cluster(i)%NNOD
426 DO j = 1,nnod
427 n = cluster(i)%NOD1(j)
428 fcluster(1,n) = cluster(i)%FOR(1)
429 fcluster(2,n) = cluster(i)%FOR(2)
430 fcluster(3,n) = cluster(i)%FOR(3)
431 n = cluster(i)%NOD2(j)
432 fcluster(1,n) = cluster(i)%FOR(1)
433 fcluster(2,n) = cluster(i)%FOR(2)
434 fcluster(3,n) = cluster(i)%FOR(3)
435 ENDDO
436 ENDDO
437 ENDIF
438 IF (anim_v(20) + h3d_data%N_VECT_CLUST_MOM > 0) THEN
439 DO i = 1, ncluster
440 nnod = cluster(i)%NNOD
441 DO j = 1,nnod
442 n = cluster(i)%NOD1(j)
443 mcluster(1,n) = cluster(i)%MOM(1)
444 mcluster(2,n) = cluster(i)%MOM(2)
445 mcluster(3,n) = cluster(i)%MOM(3)
446 n = cluster(i)%NOD2(j)
447 mcluster(1,n) = cluster(i)%MOM(1)
448 mcluster(2,n) = cluster(i)%MOM(2)
449 mcluster(3,n) = cluster(i)%MOM(3)
450 ENDDO
451 ENDDO
452 ENDIF
453
454 IF (nindx > 0) THEN
455 DO j=1,nindx
456#include "lockon.inc"
457 WRITE(iout ,1000) cluster(indx(j))%ID
458 WRITE(istdo,1100) cluster(indx(j))%ID,tt
459#include "lockoff.inc"
460 END DO
461 ENDIF
462
463 1000 FORMAT(5x,'DELETE ELEMENT CLUSTER,ID=',i10)
464 1100 FORMAT(5x,'DELETE ELEMENT CLUSTER,ID=',i10,', AT TIME ',1pe16.9)
465
466 RETURN
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB