56
57
58
60
61
62
63#include "implicit_f.inc"
64
65
66
67#include "mvsiz_p.inc"
68#include "impl1_c.inc"
69#include "comlock.inc"
70
71
72
73 INTEGER JLT, JLT_NEW,NIN,NSN,INACTI,IEDGE, NFT,
74 . CAND_N(*),CN_LOC(MVSIZ),
75 . CAND_E(*),CE_LOC(MVSIZ), NSVG(MVSIZ), KINI(*), CAND_T(*)
76 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
77 . INTTH,MSEGLO(*),IRTLM(2,NSN) ,SUBTRIA(MVSIZ),
78 . NSMS(*),IXX(MVSIZ,13),MVOISIN(4,*),ITRIV(4,MVSIZ),
79 . INTPLY ,NSNE,ISPT2(*),IZERO
81 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pmax_gap,
82 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz
83 . kb1(mvsiz), kb2(mvsiz), kb3(mvsiz), kb4(mvsiz),
84 . kc1(mvsiz), kc2(mvsiz), kc3(mvsiz), kc4(mvsiz),
85 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
86 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
87 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
88 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
89 . vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),vy1(mvsiz),
90 . vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),vz1(mvsiz),vz2(mvsiz),
91 . vz3(mvsiz),vz4(mvsiz),vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),
92 . time_s(nsn),pene(mvsiz),secnd_fr(6,nsn),
93 . xx0(mvsiz,17),yy0(mvsiz,17),zz0(mvsiz,17),gap_m(*),
94 . vx(mvsiz,17),vy(mvsiz,17),vz(mvsiz,17),gaps(*),gap_nm(12,*),
95 . pene_old(5,nsn),stif_old(2,nsn),penmin,eps0,
96 . nm1(mvsiz), nm2(mvsiz), nm3(mvsiz),dgap_nm(4,*),marge
97 INTEGER IRECT(4,*),ITAB(*),ICONT_I(*)
98 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IKNON
99 my_real,
DIMENSION(MVSIZ),
INTENT(OUT) :: penref
100
101
102
103#include "scr05_c.inc"
104#include "com08_c.inc"
105
106
107
108 INTEGER I, IG, J, K, L, INTERSECT,MG,N,II,NSLID,ITR,ITS,ITT,ITQ,NS
110 . x5(mvsiz), y5(mvsiz), z5(mvsiz),
111 . vx5(mvsiz), vy5(mvsiz), vz5(mvsiz),
112 . al1(mvsiz), al2(mvsiz), al3(mvsiz), al4(mvsiz),
113 . nxx(mvsiz,17),nyy(mvsiz,17),nzz(mvsiz,17),
114 . x51, x52, x53, x54,
115 . y51, y52, y53, y54,
116 . z51, z52, z53, z54,
117 . xo16, xo17, xo15, xo14, x163, x153, x152, x142, x141,
118 . yo16, yo17, yo15, yo14, y163, y153, y152, y142, y141,
119 . zo16, zo17, zo15, zo14, z163, z153, z152, z142, z141,
120 . x164, x174, x171,
121 . y164, y174, y171,
122 . z164, z174, z171,
123 . xi1, xi2, xi3, xi4, xi5,
124 . yi1, yi2, yi3, yi4, yi5,
125 . zi1, zi2, zi3, zi4, zi5,
126 . xia, xib, xic,
127 . yia, yib, yic,
128 . zia, zib, zic,
129 . xo1, xo2, xo3, xo4, xo5, xoi,
130 . yo1, yo2, yo3, yo4, yo5, yoi,
131 . zo1, zo2, zo3, zo4, zo5, zoi,xs,ys,zs,
132 . xm1, xm2, xm3, xm4, xm5,
133 . ym1, ym2, ym3, ym4, ym5,
134 . zm1, zm2, zm3, zm4, zm5
136 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
137 . sx125(mvsiz),sx235(mvsiz),sx345(mvsiz),sx415(mvsiz),
138 . sy125(mvsiz),sy235(mvsiz),sy345(mvsiz),sy415(mvsiz),
139 . sz125(mvsiz),sz235(mvsiz),sz345(mvsiz),sz415(mvsiz),
140 . sx1250(mvsiz),sx2350(mvsiz),sx3450(mvsiz),sx4150(mvsiz),
141 . sy1250(mvsiz),sy2350(mvsiz),sy3450(mvsiz),sy4150(mvsiz),
142 . sz1250(mvsiz),sz2350(mvsiz),sz3450(mvsiz),sz4150(mvsiz),
143 . sx2114(mvsiz),sx3215(mvsiz),sx4316(mvsiz),sx1417(mvsiz),
144 . sy2114(mvsiz),sy3215(mvsiz),sy4316(mvsiz),sy1417(mvsiz),
145 . sz2114(mvsiz),sz3215(mvsiz),sz4316(mvsiz),sz1417(mvsiz),
146 . sx21140(mvsiz),sx32150(mvsiz),sx43160(mvsiz),sx14170(mvsiz),
147 . sy21140(mvsiz),sy32150(mvsiz),sy43160(mvsiz),sy14170(mvsiz),
148 . sz21140(mvsiz),sz32150(mvsiz),sz43160(mvsiz),sz14170(mvsiz),
149 . la(mvsiz),lb(mvsiz), lc(mvsiz),
150 . xa(mvsiz),ya(mvsiz),za(mvsiz),
151 . xb(mvsiz),yb(mvsiz),zb(mvsiz),
152 . xc(mvsiz),yc(mvsiz),zc(mvsiz),
153 . xd(mvsiz),yd(mvsiz),zd(mvsiz),
154 . xe(mvsiz),ye(mvsiz),ze(mvsiz),
155 . xf(mvsiz),yf(mvsiz),zf(mvsiz),adt0(mvsiz),
156 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),
157 . xxl(17),yyl(17),zzl(17)
159 . s2,d1,d2,d3,d4,unp,zerom,h0,nqx,nqy,nqz,la2,lb2,lc2,
160 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
161 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
162 . sz,sz1,sz2,sz3,sz4,sz5,sz0,
163 . v12,v23,v34,v41,
164 . nx1, nx2, nx3, nx4,
165 . ny1, ny2, ny3, ny4,
166 . nz1, nz2, nz3, nz4,
167 . nx14, nx15, nx16, nx17,
168 . ny14, ny15, ny16, ny17,
169 . nz14, nz15, nz16, nz17,
170 . sox125,sox235,sox345,sox415,
171 . soy125,soy235,soy345,soy415,
172 . soz125,soz235,soz345,soz415,
173 . sox2114,sox3215,sox4316,sox1417,
174 . soy2114,soy3215,soy4316,soy1417,
175 . soz2114,soz3215,soz4316,soz1417,
176 . xo12,xo23,xo34,xo41,sxo1,sxo2,sxo3,sxo4,sxo5,
177 . yo12,yo23,yo34,yo41,syo1,syo2,syo3,syo4,syo5,
178 . zo12,zo23,zo34,zo41,szo1,szo2,szo3,szo4,szo5,
179 . sxo12,sxo23,sxo34,sxo41,vo12,vo23,vo34,vo41,
180 . syo12,syo23,syo34,syo41,
181 . szo12,szo23,szo34,szo41,
182 . x13_2,y13_2,z13_2,x7_3 ,y7_3 ,z7_3 ,
183 . x9_4 ,y9_4 ,z9_4 ,x11_1,y11_1,z11_1,
184 . x6_4 ,y6_4 ,z6_4 ,x8_1 ,y8_1 ,z8_1 ,
185 . x10_2,y10_2,z10_2,x12_3,y12_3,z12_3,
186 . prx,pry,prz,psx,psy,psz,ptx,pty,ptz,
187 . nqrx,nqry,nqrz,nqsx,nqsy,nqsz,nqtx,nqty,nqtz,
188 . nrx,nry,nrz,nsx,nsy,nsz,ntx,nty,ntz,
189 . nax,nay,naz,nbx,nby,nbz,ncx,ncy,ncz,
190 . sax,say,saz,sbx,sby,sbz,scx,scy,scz,
191 . trx,try,trz,tsx,tsy,tsz,ttx,tty,ttz,
192 . tr2,ts2,tt2,aaa,bbb,vr,vs,vt,nnx,nny,nnz,
193 . xab,xbc,xca,yab,ybc,yca,zab,zbc,zca,
194 . xce,yce,zce,xaf,yaf,zaf,
195 . xpa,ypa,zpa,xpb,ypb,zpb,xpc,ypc,zpc,xpi,ypi,zpi,
196 . abx,bcx,cax,aby,bcy,cay,abz,bcz,caz,xbd,ybd,zbd
198 . gap2, ds2,t1,t2,t3,xxx,yyy,zzz,nni,ni2,
199 . al1num,al2num,al3num,al4num
200 . x23d,y23d,z23d,x34d,y34d,z34d,x41d,y41d,z41d,
201 . x12d,y12d,z12d,gap2d,xi5d,yi5d,zi5d,s2d, la0,lb0,lc0,
202 . hla, adt,overw,eps,pene_sh,overw0,vol,eps2,dgap,rx,ry,rz
203 . pene_tm(mvsiz),la_tm(mvsiz),lb_tm(mvsiz),n_tm(3,mvsiz),
204 .
area(mvsiz),f_pene,fac_k,fac_p
205 INTEGER ICONTACT(MVSIZ),IT0(3,20),IT(3,20,MVSIZ),IC(10,20),
206 . ISTEP(MVSIZ),IRTLM_OLD(MVSIZ),SUBTRIA_N(MVSIZ),
207 . NSS,MGLOB,ibug,JG,IZLIM,ip,IKEEP,LARGEP,IPEN0(MVSIZ),IER,
208 . NOD1,NOD2,NOR_OLD,ISHEL(MVSIZ),ICONT_R(MVSIZ),
209 . SUBTRIA_TM(MVSIZ),IFIRST
211 . unpt,zeromt,eps02,tole,epseg,marge025,tol_ints,tol_b,epsext
212
213
214
215
216
217
218
219
220
221
222
223! | 6 | 15 3 2 | 5 8 9 | 8 9 3 2 | |/ 7 \|
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242 DATA ic /
243 1 5, 1, 2 , 14, 3, 4 , 3, 4, 1, 2,
244 2 5, 2, 3 , 15, 4, 1 , 4, 1, 2, 3,
245 3 5, 3, 4 , 16, 1, 2 , 1, 2, 3, 4,
246 4 5, 4, 1 , 17, 2, 3 , 2, 3, 4, 1,
247 5 14, 2, 1 , 5, 6, 7 , 6, 7, 2, 1,
248 6 15, 3, 2 , 5, 8, 9 , 8, 9, 3, 2,
249 7 16, 4, 3 , 5,10,11 , 10,11, 4, 3,
250 8 17, 1, 4 , 5,12,13 , 12,13, 1, 4,
251 9 14, 1, 6 , 0, 7, 2 , 7, 2, 1, 6,
252 . 15, 2, 8 , 0, 9, 3 , 9, 3, 2, 8,
253 1 16, 3,10 , 0,11, 4 , 11, 4, 3,10,
254 2 17, 4,12 , 0,13, 1 , 13, 1, 4,12,
255 3 14, 7, 2 , 0, 1, 6 , 1, 6, 7, 2,
256 4 15, 9, 3 , 0, 2, 8 , 2, 8, 9, 3,
257 5 16,11, 4 , 0, 3,10 , 3,10,11, 4,
258 6 17,13, 1 , 0, 4,12 , 4,12,13, 1,
259 7 14, 6, 7 , 0, 2, 1 , 2, 1, 6, 7,
260 8 15, 8, 9 , 0, 3, 2 , 3, 2, 8, 9,
261 9 16,10,11 , 0, 4, 3 , 4, 3,10,11,
262 . 17,12,13 , 0, 1, 4 , 1, 4,12,13/
263
264
265
266
267
268! | +----------+----------+----------+----------+
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296 DATA it0 /
297 1 5, 2, 4,
298 2 6, 3, 1,
299 3 7, 4, 2,
300 4 8, 1, 3,
301 5 1, 9,13,
302 6 2,10,14,
303 7 3,11,15,
304 8 4,12,16,
305 9 -1,17, 5,
306 . -1,18, 6,
307 1 -1,19, 7,
308 2 -1,20, 8,
309 3 -1, 5,17,
310 4 -1, 6,18,
311 5 -1, 7,19,
312 6 -1, 8,20,
313 7 -1, 9,13,
314 8 -1,10,14,
315 9 -1,11,15,
316 . -1,12,16/
317
318
319
320
321 nor_old = 0
322 unp = one + em4
323 zerom = zero - em4
324 eps = (two+half)/hundred
325 unpt = one + eps
326 zeromt = zero - eps
327
328 eps2=eps*eps
329
330 eps = em3
331 eps02=eps*eps
332 eps = eps0
333 overw0=one
334 marge025=fourth*marge
335
336
337 IF (iresp==1) THEN
338 tol_ints = two*em3
339 tol_b = two*em6
340 ELSE
341 tol_ints = em08
342 tol_b = em08
343 END IF
344 nod1=0
345 nod2=0
346
347
348
349 DO i=1,jlt
350 n=cand_n(i)
351
352
353
354
355
356 ipen0(i)=0
357 icont_r(i) = 0
358 IF(n <= nsn)THEN
359 icontact(i) = irtlm(1,n)
360 subtria(i)=irtlm(2,n)
361 ELSE
362 icontact(i) =
irtlm_fi(nin)%P(1,n-nsn
363 subtria(i) =
irtlm_fi(nin)%P(2,n-nsn)
364 ENDIF
365 irtlm_old(i) = icontact(i)
366
367 IF(icontact(i) <= zero)THEN
368 istep(i) = 2
369 icontact(i) = 0
370 ELSEIF(icontact(i) == mseglo(cand_e(i)))THEN
371 istep(i) = 1
372 ELSE
373 icontact(i) = 0
374 istep(i) = 0
375 ENDIF
376 IF(stif(i) == zero)THEN
377 icontact(i) = 0
378 istep(i) = 0
379 ENDIF
380 IF(n <= nsn)THEN
381 IF (pene_old(5,n)> zero.AND.istep(i) == 1) ipen0(i)=1
382 ELSE
383 IF (
pene_oldfi(nin)%P(5,n-nsn)> zero.AND.istep(i) == 1) ipen0
384 ENDIF
385
386 IF(istep(i) == 1.AND.subtria(i) > 20) THEN
387 subtria(i) = subtria(i) - 20
388 istep(i) = 6
389 ENDIF
390 IF(iedge/=0)THEN
391 IF(cand_t(i)==2)THEN
392
393 icontact(i) = 0
394 istep(i) = 0
395 ENDIF
396 ENDIF
397 nm1(i) = zero
398 nm2(i) = zero
399 nm3(i) = zero
400 ENDDO
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427 DO i=1,jlt
428 ishel(i)=0
429 IF(istep(i) == 0)cycle
430
431
432 xo14 = xx0(i,14)
433 yo14 = yy0(i,14)
434 zo14 = zz0(i,14)
435
436 xo15 = xx0(i,15)
437 yo15 = yy0(i,15)
438 zo15 = zz0(i,15)
439
440 xo16 = xx0(i,16)
441 yo16 = yy0(i,16)
442 zo16 = zz0(i,16)
443
444 xo17 = xx0(i,17)
445 yo17 = yy0(i,17)
446 zo17 = zz0(i,17)
447
448 x51 = xx0(i,1) - xx0(i,5)
449 y51 = yy0(i,1) - yy0(i,5)
450 z51 = zz0(i,1) - zz0(i,5)
451
452 x52 = xx0(i,2) - xx0(i,5)
453 y52 = yy0(i,2) - yy0(i,5)
454 z52 = zz0(i,2) - zz0(i,5)
455
456 x53 = xx0(i,3) - xx0(i,5)
457 y53 = yy0(i,3) - yy0(i,5)
458 z53 = zz0(i,3) - zz0(i,5)
459
460 x54 = xx0(i,4) - xx0(i,5)
461 y54 = yy0(i,4) - yy0(i,5)
462 z54 = zz0(i,4) - zz0(i,5)
463
464
465 sx1250(i) = y51*z52 - z51*y52
466 sy1250(i) = z51*x52 - x51*z52
467 sz1250(i) = x51*y52 - y51*x52
468
469 sx2350(i) = y52*z53 - z52*y53
470 sy2350(i) = z52*x53 - x52*z53
471 sz2350(i) = x52*y53 - y52*x53
472
473 sx3450(i) = y53*z54 - z53*y54
474 sy3450(i) = z53*x54 - x53*z54
475 sz3450(i) = x53*y54 - y53*x54
476
477 sx4150(i) = y54*z51 - z54*y51
478 sy4150(i) = z54*x51 - x54*z51
479 sz4150(i) = x54*y51 - y54*x51
480
481
482 x141 = xx0(i,1) - xx0(i,14)
483 y141 = yy0(i,1) - yy0(i,14)
484 z141 = zz0(i,1) - zz0(i,14)
485
486 x142 = xx0(i,2) - xx0(i,14)
487 y142 = yy0(i,2) - yy0(i,14)
488 z142 = zz0(i,2) - zz0(i,14)
489
490 x152 = xx0(i,2) - xx0(i,15)
491 y152 = yy0(i,2) - yy0(i,15)
492 z152 = zz0(i,2) - zz0(i,15)
493
494 x153 = xx0(i,3) - xx0(i,15)
495 y153 = yy0(i,3) - yy0(i,15)
496 z153 = zz0(i,3) - zz0(i,15)
497
498 x163 = xx0(i,3) - xx0(i,16)
499 y163 = yy0(i,3) - yy0(i,16)
500 z163 = zz0(i,3) - zz0(i,16)
501
502 x164 = xx0(i,4) - xx0(i,16)
503 y164 = yy0(i,4) - yy0(i,16)
504 z164 = zz0(i,4) - zz0(i,16)
505
506 x174 = xx0(i,4) - xx0(i,17)
507 y174 = yy0(i,4) - yy0(i,17)
508 z174 = zz0(i,4) - zz0(i,17)
509
510 x171 = xx0(i,1) - xx0(i,17)
511 y171 = yy0(i,1) - yy0(i,17)
512 z171 = zz0(i,1) - zz0(i,17)
513
514
515 x13_2 = xx0(i,2) - xx0(i,13)
516 y13_2 = yy0(i,2) - yy0(i,13)
517 z13_2 = zz0(i,2) - zz0(i,13)
518
519 x7_3 = xx0(i,3) - xx0(i,7)
520 y7_3 = yy0(i,3) - yy0(i,7)
521 z7_3 = zz0(i,3) - zz0(i,7)
522
523 x9_4 = xx0(i,4) - xx0(i,9)
524 y9_4 = yy0(i,4) - yy0(i,9)
525 z9_4 = zz0(i,4) - zz0(i,9)
526
527 x11_1 = xx0(i,1) - xx0(i,11)
528 y11_1 = yy0(i,1) - yy0(i,11)
529 z11_1 = zz0(i,1) - zz0(i,11)
530
531 x6_4 = xx0(i,4) - xx0(i,6)
532 y6_4 = yy0(i,4) - yy0(i,6)
533 z6_4 = zz0(i,4) - zz0(i,6)
534
535 x8_1 = xx0(i,1) - xx0(i,8)
536 y8_1 = yy0(i,1) - yy0(i,8)
537 z8_1 = zz0(i,1) - zz0(i,8)
538
539 x10_2 = xx0(i,2) - xx0(i,10)
540 y10_2 = yy0(i,2) - yy0(i,10)
541 z10_2 = zz0(i,2) - zz0(i,10)
542
543 x12_3 = xx0(i,3) - xx0(i,12)
544 y12_3 = yy0(i,3) - yy0(i,12)
545 z12_3 = zz0(i,3) - zz0(i,12)
546
547
548
549 IF(mvoisin(1,cand_e(i))/=0)THEN
550 sx21140(i) = y142*z141 - z142*y141
551 sy21140(i) = z142*x141 - x142*z141
552 sz21140(i) = x142*y141 - y142*x141
553 ELSE
554 sx21140(i) = sx1250(i)
555 sy21140(i) = sy1250(i)
556 sz21140(i) = sz1250(i)
557 ENDIF
558
559
560 IF(ixx(i,3) /= ixx(i,4))THEN
561 IF(mvoisin(2,cand_e(i))/=0)THEN
562 sx32150(i) = y153*z152 - z153*y152
563 sy32150(i) = z153*x152 - x153*z152
564 sz32150(i) = x153*y152 - y153*x152
565 ELSE
566 sx32150(i) = sx2350(i)
567 sy32150(i) = sy2350(i)
568 sz32150(i) = sz2350(i)
569 ENDIF
570
571 IF(mvoisin(3,cand_e(i))/=0)THEN
572 sx43160(i) = y164*z163 - z164*y163
573 sy43160(i) = z164*x163 - x164*z163
574 sz43160(i) = x164*y163 - y164*x163
575 ELSE
576 sx43160(i) = sx3450(i)
577 sy43160(i) = sy3450(i)
578 sz43160(i) = sz3450(i)
579 ENDIF
580
581 IF(mvoisin(4,cand_e(i))/=0)THEN
582 sx14170(i) = y171*z174 - z171*y174
583 sy14170(i) = z171*x174 - x171*z174
584 sz14170(i) = x171*y174 - y171*x174
585 ELSE
586 sx14170(i) = sx4150(i)
587 sy14170(i) = sy4150(i)
588 sz14170(i) = sz4150(i)
589 ENDIF
590
591
592 nxx(i,1) = y13_2 * z6_4 - z13_2 * y6_4
593 nyy(i,1) = z13_2 * x6_4 - x13_2 * z6_4
594 nzz(i,1) = x13_2 * y6_4 - y13_2 * x6_4
595
596 nxx(i,2) = y7_3 * z8_1 - z7_3 * y8_1
597 nyy(i,2) = z7_3 * x8_1 - x7_3 * z8_1
598 nzz(i,2) = x7_3 * y8_1 - y7_3 * x8_1
599
600 nxx(i,3) = y9_4 * z10_2 - z9_4 * y10_2
601 nyy(i,3) = z9_4 * x10_2 - x9_4 * z10_2
602 nzz(i,3) = x9_4 * y10_2 - y9_4 * x10_2
603
604 nxx(i,4) = y11_1 * z12_3 - z11_1 * y12_3
605 nyy(i,4) = z11_1 * x12_3 - x11_1 * z12_3
606 nzz(i,4) = x11_1 * y12_3 - y11_1 * x12_3
607
608 nxx(i,5) = fourth*(nxx(i,1)+nxx(i,2)+nxx(i,3)+nxx(i,4))
609 nyy(i,5) = fourth*(nyy(i,1)+nyy(i,2)+nyy(i,3)+nyy(i,4))
610 nzz(i,5) = fourth*(nzz(i,1)+nzz(i,2)+nzz(i,3)+nzz(i,4))
611
612 ELSE
613
614 IF(mvoisin(2,cand_e(i))/=0)THEN
615 sx32150(i) = y153*z152 - z153*y152
616 sy32150(i) = z153*x152 - x153*z152
617 sz32150(i) = x153*y152 - y153*x152
618 ELSE
619 sx32150(i) = sx1250(i)
620 sy32150(i) = sy1250(i)
621 sz32150(i) = sz1250(i)
622 ENDIF
623
624 IF(mvoisinTHEN
625 sx43160(i) = y164*z163
626 sy43160(i) = z164*x163 - x164*z163
627 sz43160(i) = x164*y163 - y164*x163
628 ELSE
629 sx43160(i) = sx1250(i)
630 sy43160(i) = sy1250(i)
631 sz43160(i) = sz1250(i)
632 ENDIF
633
634 IF(mvoisin(4,cand_e(i))/=0)THEN
635 sx14170(i) = y171*z174 - z171*y174
636 sy14170(i) = z171*x174 - x171*z174
637 sz14170(i) = x171*y174 - y171*x174
638 ELSE
639 sx14170(i) = sx1250(i)
640 sy14170(i) = sy1250(i)
641 sz14170(i) = sz1250(i)
642 ENDIF
643
644 nxx(i,1) =sx1250(i)+sx14170(i)+sx21140(i)
645 nyy(i,1) =sy1250(i)+sy14170(i)+sy21140(i)
646 nzz(i,1) =sz1250(i)+sz14170(i)+sz21140(i)
647
648 nxx(i,2) =sx1250(i)+sx21140(i)+sx32150(i)
649 nyy(i,2) =sy1250(i)+sy21140(i)+sy32150(i)
650 nzz(i,2) =sz1250(i)+sz21140(i)+sz32150(i)
651
652 nxx(i,3) =sx1250(i)+sx32150(i)+sx14170(i)
653 nyy(i,3) =sy1250(i)+sy32150(i)+sy14170(i)
654 nzz(i,3) =sz1250(i)+sz32150(i)+sz14170(i)
655
656 nxx(i,4) =nxx(i,3)
657 nyy(i,4) =nyy(i,3)
658 nzz(i,4) =nzz(i,3)
659
660 nxx(i,5) =nxx(i,3)
661 nyy(i,5) =nyy(i,3)
662 nzz(i,5) =nzz(i,3)
663
664 ENDIF
665
666
667
668 n=cand_n(i)
669 l=cand_e(i)
670
671
672
673
674
675 IF(gaps(i)>zero.OR.gap_m(l)>zero) THEN
676 ishel(i)=1
677 ELSEIF(gap_nm(5,l)>zero.OR.gap_nm(6,l)>zero.OR.gap_nm(7,l)>zero) THEN
678 ishel(i)=1
679 ELSEIF(gap_nm(8,l)>zero.OR.gap_nm(9,l)>zero.OR.gap_nm(10,l)>zero) THEN
680 ishel(i)=1
681 ELSEIF(gap_nm(11,l)>zero.OR.gap_nm(12,l)>zero) THEN
682 ishel(i)=1
683 END IF
684
685 IF(gaps(i)>zero.AND.gap_m(l)>zero) ishel(i)=2
686 IF(ishel(i)>0) THEN
687 nx1 = nxx(i,1)
688 ny1 = nyy(i,1)
689 nz1 = nzz(i,1)
690 aaa = gaps(i)+gap_nm(1,l)
691 aaa = aaa/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
692 nx1 = nx1 * aaa
693 ny1 = ny1 * aaa
694 nz1 = nz1 * aaa
695
696 nx2 = nxx(i,2)
697 ny2 = nyy(i,2)
698 nz2 = nzz(i,2)
699 aaa = gaps(i)+gap_nm(2,l)
700 aaa = aaa/sqrt(nx2*nx2+ny2*ny2+nz2*nz2)
701 nx2 = nx2 * aaa
702 ny2 = ny2 * aaa
703 nz2 = nz2 * aaa
704
705 IF(ixx(i,3) /= ixx(i,4))THEN
706 nx3 = nxx(i,3)
707 ny3 = nyy(i,3)
708 nz3 = nzz(i,3)
709 aaa = gaps(i)+gap_nm(3,l)
710 aaa = aaa/sqrt(nx3*nx3+ny3*ny3+nz3*nz3)
711 nx3 = nx3 * aaa
712 ny3 = ny3 * aaa
713 nz3 = nz3 * aaa
714
715 nx4 = nxx(i,4)
716 ny4 = nyy(i,4)
717 nz4 = nzz(i,4)
718 aaa = gaps(i)+gap_nm(4,l)
719 aaa = aaa/sqrt(nx4*nx4+ny4*ny4+nz4*nz4)
720 nx4 = nx4 * aaa
721 ny4 = ny4 * aaa
722 nz4 = nz4 * aaa
723
724 xx(i,1) = xx0(i,1) + nx1
725 yy(i,1) = yy0(i,1) + ny1
726 zz(i,1) = zz0(i,1) + nz1
727
728 xx(i,2) = xx0(i,2) + nx2
729 yy(i,2) = yy0(i,2) + ny2
730 zz(i,2) = zz0(i,2) + nz2
731
732 xx(i,3) = xx0(i,3) + nx3
733 yy(i,3) = yy0(i,3) + ny3
734 zz(i,3) = zz0(i,3) + nz3
735
736 xx(i,4) = xx0(i,4) + nx4
737 yy(i,4) = yy0(i,4) + ny4
738 zz(i,4) = zz0(i,4) + nz4
739
740 xx(i,5) = fourth*(xx(i,1)+xx(i,2)+xx(i,3)+xx(i,4))
741 yy(i,5) = fourth*(yy(i,1)+yy(i,2)+yy(i,3)+yy(i,4))
742 zz(i,5) = fourth*(zz(i,1)+zz(i,2)+zz(i,3)+zz(i,4))
743 IF(ixx(i,10)==ixx(i,11))THEN
744 xx(i,16) = xx0(i,10)
745 yy(i,16) = yy0(i,10)
746 zz(i,16) = zz0(i,10)
747 ELSE
748 xx(i,16) = fourth*(xx0(i,4)+xx0(i,3)+xx0(i,10)+xx0(i,11))
749 yy(i,16) = fourth*(yy0(i,4)+yy0(i,3)+yy0(i,10)+yy0(i,11))
750 zz(i,16) = fourth*(zz0(i,4)+zz0(i,3)+zz0(i,10)+zz0(i,11))
751 ENDIF
752 nx16 = sx43160(i)
753 ny16 = sy43160(i)
754 nz16 = sz43160(i)
755 aaa = gaps(i)+fourth*(
756 . gap_nm(3,l)+gap_nm(4,l)
757 . + gap_nm(9,l)+gap_nm(10,l) )
758 aaa = aaa/sqrt(nx16*nx16+ny16*ny16+nz16*nz16)
759 xx(i,16) = xx(i,16) + nx16 * aaa
760 yy(i,16) = yy(i,16) + ny16 * aaa
761 zz(i,16) = zz(i,16) + nz16 * aaa
762 ELSE
763 nx3 = nxx(i,3)
764 ny3 = nyy(i,3)
765 nz3 = nzz(i,3)
766 aaa = gaps(i)+gap_nm(3,l)
767 aaa = aaa/sqrt(nx3*nx3+ny3*ny3+nz3*nz3)
768 nx3 = nx3 * aaa
769 ny3 = ny3 * aaa
770 nz3 = nz3 * aaa
771
772 nx4 = nxx(i,4)
773 ny4 = nyy(i,4)
774 nz4 = nzz(i,4)
775
776 xx(i,1) = xx0(i,1) + nx1
777 yy(i,1) = yy0(i,1) + ny1
778 zz(i,1) = zz0(i,1) + nz1
779
780 xx(i,2) = xx0(i,2) + nx2
781 yy(i,2) = yy0(i,2) + ny2
782 zz(i,2) = zz0(i,2) + nz2
783
784 xx(i,3) = xx0(i,3) + nx3
785 yy(i,3) = yy0(i,3) + ny3
786 zz(i,3) = zz0(i,3) + nz3
787
788 xx(i,4) = xx(i,3)
789 yy(i,4) = yy(i,3)
790 zz(i,4) = zz(i,3)
791
792 xx(i,5) = xx(i,3)
793 yy(i,5) = yy(i,3)
794 zz(i,5) = zz(i,3)
795 ENDIF
796 IF(ixx(i,6) == ixx(i,7))THEN
797 xx(i,14) = xx0(i,6)
798 yy(i,14) = yy0(i,6)
799 zz(i,14) = zz0(i,6)
800 ELSE
801 xx(i,14) = fourth*(xx0(i,2)+xx0(i,1)+xx0(i,6)+xx0(i,7))
802 yy(i,14) = fourth*(yy0(i,2)+yy0(i,1)+yy0(i,6)+yy0(i,7))
803 zz(i,14) = fourth*(zz0(i,2)+zz0(i,1)+zz0(i,6)+zz0(i,7))
804 ENDIF
805 nx14 = sx21140(i)
806 ny14 = sy21140(i)
807 nz14 = sz21140(i)
808 aaa = gaps(i)+fourth*(
809 . gap_nm(1,l)+gap_nm(2,l)
810 . + gap_nm(5,l)+gap_nm(6,l) )
811 aaa = aaa/sqrt(nx14*nx14+ny14*ny14+nz14*nz14)
812 xx(i,14) = xx(i,14) + nx14 * aaa
813 yy(i,14) = yy(i,14) + ny14 * aaa
814 zz(i,14) = zz(i,14) + nz14 * aaa
815 IF(ixx(i,8)==ixx(i,9))THEN
816 xx(i,15) = xx0(i,8)
817 yy(i,15) = yy0(i,8)
818 zz(i,15) = zz0(i,8)
819 ELSE
820 xx(i,15) = fourth*(xx0(i,3)+xx0(i,2)+xx0(i,8)+xx0(i,9))
821 yy(i,15) = fourth*(yy0(i,3)+yy0(i,2)+yy0(i,8)+yy0(i,9))
822 zz(i,15) = fourth*(zz0(i,3)+zz0(i,2)+zz0(i,8)+zz0(i,9))
823 ENDIF
824 nx15 = sx32150(i)
825 ny15 = sy32150(i)
826 nz15 = sz32150(i)
827 aaa = gaps(i)+fourth*(
828 . gap_nm(2,l)+gap_nm(3,l)
829 . + gap_nm(7,l)+gap_nm(8,l) )
830 aaa = aaa/sqrt(nx15*nx15+ny15*ny15+nz15*nz15)
831 xx(i,15) = xx(i,15) + nx15 * aaa
832 yy(i,15) = yy(i,15) + ny15 * aaa
833 zz(i,15) = zz(i,15) + nz15 * aaa
834 IF(ixx(i,12)==ixx(i,13))THEN
835 xx(i,17) = xx0(i,12)
836 yy(i,17) = yy0(i,12)
837 zz(i,17) = zz0(i,12)
838 ELSE
839 xx(i,17) = fourth*(xx0(i,1)+xx0(i,4)+xx0(i,12)+xx0(i,13))
840 yy(i,17) = fourth*(yy0(i,1)+yy0(i,4)+yy0(i,12)+yy0(i,13))
841 zz(i,17) = fourth*(zz0(i,1)+zz0(i,4)+zz0(i,12)+zz0(i,13))
842 ENDIF
843 nx17 = sx14170(i)
844 ny17 = sy14170(i)
845 nz17 = sz14170(i)
846 aaa = gaps(i)+fourth*(
847 . gap_nm(4,l)+gap_nm(1,l)
848 . + gap_nm(11,l)+gap_nm(12,l) )
849 aaa = aaa/sqrt(nx17*nx17+ny17*ny17+nz17*nz17)
850 xx(i,17) = xx(i,17) + nx17 * aaa
851 yy(i,17) = yy(i,17) + ny17 * aaa
852 zz(i,17) = zz(i,17) + nz17 * aaa
853
854 x51 = xx(i,1) - xx(i,5)
855 y51 = yy(i,1) - yy(i,5)
856 z51 = zz(i,1) - zz(i,5)
857
858 x52 = xx(i,2) - xx(i,5)
859 y52 = yy(i,2) - yy(i,5)
860 z52 = zz(i,2) - zz(i,5)
861
862 x53 = xx(i,3) - xx(i,5)
863 y53 = yy(i,3) - yy(i,5)
864 z53 = zz(i,3) - zz(i,5)
865
866 x54 = xx(i,4) - xx(i,5)
867 y54 = yy(i,4) - yy(i,5)
868 z54 = zz(i,4) - zz(i,5)
869
870
871 sx125(i) = y51*z52 - z51*y52
872 sy125(i) = z51*x52 - x51*z52
873 sz125(i) = x51*y52 - y51*x52
874
875 sx235(i) = y52*z53 - z52*y53
876 sy235(i) = z52*x53 - x52*z53
877 sz235(i) = x52*y53 - y52*x53
878
879 sx345(i) = y53*z54 - z53*y54
880 sy345(i) = z53*x54 - x53*z54
881 sz345(i) = x53*y54 - y53*x54
882
883 sx415(i) = y54*z51 - z54*y51
884 sy415(i) = z54*x51 - x54*z51
885 sz415(i) = x54*y51 - y54*x51
886
887 ELSE
888
889 xx(i,1:17) = xx0(i,1:17)
890 yy(i,1:17) = yy0(i,1:17)
891 zz(i,1:17) = zz0(i,1:17)
892
893 sx125(i)=sx1250(i)
894 sy125(i)=sy1250(i)
895 sz125(i)=sz1250(i)
896
897 sx235(i) = sx2350(i)
898 sy235(i) = sy2350(i)
899 sz235(i) = sz2350(i)
900
901 sx345(i) = sx3450(i)
902 sy345(i) = sy3450(i)
903 sz345(i) = sz3450(i)
904
905 sx415(i) = sx4150(i)
906 sy415(i) = sy4150(i)
907 sz415(i) = sz4150(i)
908 ENDIF
909
910
911 ENDDO
912
913
914
915
916
917 DO i=1,jlt
918 pene(i) = zero
919 pene_tm(i) = zero
920 IF(istep(i) == 0)cycle
921 IF(istep(i) /= 1.AND.tt /= zero)cycle
922
923
924 IF (tt==zero) THEN
925 itq = 1
926 ELSE
927 itq = subtria(i)
928 END IF
929 k = ic(1,itq)
930 xa(i) = xx(i,k)
931 ya(i) = yy(i,k)
932 za(i) = zz(i,k)
933
934 k = ic(2,itq)
935 xb(i) = xx(i,k)
936 yb(i) = yy(i,k)
937 zb(i) = zz(i,k)
938
939 k = ic(3,itq)
940 xc(i) = xx(i,k)
941 yc(i) = yy(i,k)
942 zc(i) = zz(i,k)
943
944 k = ic(4,itq)
945 xd(i) = xx(i,k)
946 yd(i) = yy(i,k)
947 zd(i) = zz(i,k)
948
949 k = 15
950 xe(i) = xx(i,k)
951 ye(i) = yy(i,k)
952 ze(i) = zz(i,k)
953
954 k = 17
955 xf(i) = xx(i,k)
956 yf(i) = yy(i,k)
957 zf(i) = zz(i,k)
958
959 IF(itq ==1)THEN
960 nqx = sx125(i)
961 nqy = sy125(i)
962 nqz = sz125(i)
963 ELSEIF(itq ==2)THEN
964 nqx = sx235(i)
965 nqy = sy235(i)
966 nqz = sz235(i)
967 ELSEIF(itq ==3)THEN
968 nqx = sx345(i)
969 nqy = sy345(i)
970 nqz = sz345(i)
971 ELSEIF(itq ==4)THEN
972 nqx = sx415(i)
973 nqy = sy415(i)
974 nqz = sz415(i)
975 ENDIF
976
977 s2 =
max(em30,sqrt(nqx*nqx + nqy*nqy + nqz*nqz))
979 aaa = one/s2
980 nqx = nqx*aaa
981 nqy = nqy*aaa
982 nqz = nqz*aaa
983
984 bbb = (xi(i)-xa(i))*nqx+(yi(i)-ya(i))*nqy+(zi(i)-za(i))*nqz
985
986 IF (tt == zero .AND. istep(i)/= 1) THEN
987 IF (abs(bbb) < penmin.AND.bbb/=zero) THEN
988
989 CALL s_in_m4(xx(i,1),yy(i,1),zz(i,1),xx(i,2),yy(i,2),zz(i,2),
990 1 xx(i,3),yy(i,3),zz(i,3),xx(i,4),yy(i,4),zz(i,4),
991 1 xi(i),yi(i),zi(i),ier)
992 IF (ier ==0) THEN
993 pene(i) = abs(bbb)
994 n=cand_n(i)
995#include "lockon.inc"
996 IF(n <= nsn)THEN
997 IF(irtlm(1,n)==0)THEN
998 irtlm(1,n)=-mseglo(cand_e(i))
999 irtlm(2,n)=1
1000 time_s(n) = pene(i)
1001 ELSE
1002 IF( time_s(n)<pene(i) .OR.
1003 * (time_s(n)==pene(i).AND.-irtlm(1,n) < mseglo(cand_e(i)) ) )THEN
1004 irtlm(1,n)=-mseglo(cand_e(i))
1005 irtlm(2,n)=1
1006 time_s(n) = pene(i)
1007 ENDIF
1008 ENDIF
1009 ELSE
1010 IF(
irtlm_fi(nin)%P(1,n-nsn) == 0)
THEN
1011 irtlm_fi(nin)%P(1,n-nsn)=-mseglo(cand_e(i))
1014 ELSE
1015 IF(
time_sfi(nin)%P(n-nsn)<pene(i) .OR.
1016 * (
time_sfi(nin)%P(n-nsn)==pene(i).AND.-
irtlm_fi(nin)%P(1,n-nsn) < mseglo(cand_e(i)) ) )
THEN
1017 irtlm_fi(nin)%P(1,n-nsn)=-mseglo(cand_e(i))
1020 ENDIF
1021 ENDIF
1022 ENDIF
1023#include "lockoff.inc"
1024 istep(i)= 1
1025 subtria(i) =1
1026 ELSE
1027 istep(i)= 0
1028 cycle
1029 END IF
1030 ELSE
1031 cycle
1032 END IF
1033 ELSE
1034 pene(i) = -
min(zero,bbb)
1035 END IF
1036 largep = 0
1037 IF(pene(i)*pene(i) > eps2*s2.AND.ipen0(i)==0) largep = 1
1038 pene_tm(i) =pene(i)
1039
1040
1041 n1(i) = nqx
1042 n2(i) = nqy
1043 n3(i) = nqz
1044 n_tm(1,i) = n1(i)
1045 n_tm(2,i) = n2(i)
1046 n_tm(3,i) = n3(i)
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067 xab = xb(i)-xa(i)
1068 yab = yb(i)-ya(i)
1069 zab = zb(i)-za(i)
1070 xbc = xc(i)-xb(i)
1071 ybc = yc(i)-yb(i)
1072 zbc = zc(i)-zb(i)
1073 xca = xa(i)-xc(i)
1074 yca = ya(i)-yc(i)
1075 zca = za(i)-zc(i)
1076
1077 xbd = xd(i)-xb(i)
1078 ybd = yd(i)-yb(i)
1079 zbd = zd(i)-zb(i)
1080 xce = xe(i)-xc(i)
1081 yce = ye(i)-yc(i)
1082 zce = ze(i)-zc(i)
1083 xaf = xf(i)-xa(i)
1084 yaf = yf(i)-ya(i)
1085 zaf = zf(i)-za(i)
1086
1087 xia = xa(i)-xi(i)
1088 yia = ya(i)-yi(i)
1089 zia = za(i)-zi(i)
1090 xib = xb(i)-xi(i)
1091 yib = yb(i)-yi(i)
1092 zib = zb(i)-zi(i)
1093 xic = xc(i)-xi(i)
1094 yic = yc(i)-yi(i)
1095 zic = zc(i)-zi(i)
1096 sx = - yab*zca + zab*yca
1097 sy = - zab*xca + xab*zca
1098 sz = - xab*yca + yab*xca
1099 s2 = sx*sx+sy*sy+sz*sz
1100 sax = yib*zic - zib*yic
1101 say = zib*xic - xib*zic
1102 saz = xib*yic - yib*xic
1103 la(i) = (sx*sax+sy*say+sz*saz)/s2
1104 sbx = yic*zia - zic*yia
1105 sby = zic*xia - xic*zia
1106 sbz = xic*yia - yic*xia
1107 lb(i) = (sx*sbx+sy*sby+sz*sbz)/s2
1108 lc(i) = one - la(i) - lb(i)
1109 izlim=0
1110 subtria_tm(i) = subtria(i)
1111 la_tm(i) = la(i)
1112 lb_tm(i) = lb(i)
1113 epseg =eps
1114 epsext = eps
1115
1116 IF (nsne >0) epsext=em04
1117 IF (ispt2(i) >0) epseg=em04
1118 IF(ixx(i,3) /= ixx(i,4))THEN
1119
1120
1121
1122 aaa = zerom
1124 istep(i)=3
1125 subtria_n(i)=it0(2,itq)
1126 subtria(i)=it0(2,itq)
1127 xxl(1:17)=xx(i,1:17)
1128 yyl(1:17)=yy(i,1:17)
1129 zzl(1:17)=zz(i,1:17)
1130 izlim=-1
1131 gap2 = gaps(i)+gap_m(cand_e(i))
1133 1 izlim ,istep(i),subtria_n(i),subtria(i),
1134 2 largep ,pene(i) ,xxl ,yyl ,zzl ,
1135 3 sx125(i),sy125(i),sz125(i),sx235(i),sy235(i),
1136 4 sz235(i),sx345(i),sy345(i),sz345(i),sx415(i),
1137 5 sy415(i),sz415(i),xi(i) ,yi(i) ,zi(i) ,
1138 6 n1(i) ,n2(i) ,n3(i) ,la(i) ,lb(i) ,
1139 7 lc(i) ,gap2 ,epseg )
1140 subtria_tm(i) = subtria(i)
1141
1142 ELSEIF(lc(i)<aaa)
THEN
1143 istep(i)=3
1144 subtria_n(i)=it0(3,itq)
1145 subtria(i)=it0(3,itq)
1146 xxl(1:17)=xx
1147 yyl(1:17)=yy(i,1:17)
1148 zzl(1:17)=zz(i,1:17)
1149 izlim=-1
1150 gap2 = gaps(i)+gap_m(cand_e(i))
1152 1 izlim ,istep(i),subtria_n(i),subtria(i),
1153 2 largep ,pene(i) ,xxl ,yyl ,zzl ,
1154 3 sx125(i),sy125(i),sz125(i),sx235(i),sy235(i),
1155 4 sz235(i),sx345(i),sy345(i),sz345(i),sx415(i),
1156 5 sy415(i),sz415(i),xi(i) ,yi(i) ,zi(i) ,
1157 6 n1(i) ,n2(i) ,n3(i) ,la(i) ,lb(i) ,
1158 7 lc(i) ,gap2 ,epseg )
1159 subtria_tm(i) = subtria(i)
1160
1161 ELSEIF(la(i)<-epsext.OR.(pene(i)==zero.AND.la(i)<zero))THEN
1162
1163 istep(i)=3
1164 subtria_n(i)=it0(1,itq)
1165 izlim=-1
1166 ELSEIF(la(i)<epseg)THEN
1167
1168 vol = n1(i)*xbd + n2(i)*ybd + n3(i)*zbd
1169 gap2 = gaps(i)+gap_m(cand_e(i))
1170 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1171 IF(largep == 1)THEN
1172
1173 istep(i)=3
1174 subtria_n(i)=it0(1,itq)
1175 izlim = -1
1176 ELSEIF(vol < zero)THEN
1177
1178 izlim = 1
1179 ELSEIF(la(i)<zerom)THEN
1180 istep(i)=3
1181 subtria_n(i)=it0(1,itq)
1182 izlim = -1
1183 ENDIF
1184 ENDIF
1185 IF(izlim == -1)cycle
1186
1187 ELSE
1188
1189
1190 IF(la(i)<-epsext.OR.(pene(i)==zero.AND.la(i)<zero))THEN
1191 istep(i)=3
1192 subtria_n(i)=5
1193 izlim = -1
1194 ELSEIF(la(i)<epseg)THEN
1195
1196 vol = n1(i)*xbd + n2(i)*ybd + n3(i)*zbd
1197 gap2 = gaps(i)+gap_m(cand_e(i))
1198 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1199 IF(largep == 1)THEN
1200
1201 istep(i)=3
1202 subtria_n(i)=5
1203 izlim = -1
1204 ELSEIF(vol < zero)THEN
1205
1206 izlim=1
1207 ELSEIF(la(i)<zerom)THEN
1208 istep(i)=3
1209 subtria_n(i)=5
1210 izlim = -1
1211 ENDIF
1212 ENDIF
1213 IF(lb(i)<-epsext.OR.(pene(i)==zero.AND.lb(i)<zero))THEN
1214 istep(i)=3
1215 IF(lb(i) < la(i)) subtria_n(i)=6
1216 izlim = -1
1217 ELSEIF(lb(i)<epseg)THEN
1218
1219 vol = n1(i)*xce + n2(i)*yce + n3(i)*zce
1220 gap2 = gaps(i)+gap_m(cand_e(i))
1221 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1222 IF(largep == 1)THEN
1223
1224 istep(i)=3
1225 IF(lb(i) < la(i).OR.izlim==1) subtria_n(i)=6
1226 izlim = -1
1227 ELSEIF(vol < zero)THEN
1228
1229 izlim=1
1230 ELSEIF(lb(i)<zerom)THEN
1231 istep(i)=3
1232 IF(lb(i) < la(i).OR.izlim==1) subtria_n(i)=6
1233 izlim = -1
1234 ENDIF
1235 ENDIF
1236 IF(lc(i)<-epsext.OR.(pene(i)==zero.AND.lc(i)<zero))THEN
1237 istep(i)=3
1238 IF(lc(i) < la(i).AND.lc(i) < lb(i)) subtria_n(i)=8
1239 izlim = -1
1240 ELSEIF(lc(i)<epseg)THEN
1241
1242 vol = n1(i)*xaf + n2(i)*yaf + n3(i)*zaf
1243 gap2 = gaps(i)+gap_m(cand_e(i))
1244 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1245 IF(largep == 1)THEN
1246
1247 istep(i)=3
1248 IF((lc(i) < la(i).AND.lc(i) < lb(i)).OR.izlim==1) subtria_n(i)=8
1249 izlim = -1
1250 ELSEIF(vol < zero)THEN
1251
1252 izlim=1
1253 ELSEIF(lc(i)<zerom)THEN
1254 istep(i)=3
1255 IF((lc(i) < la(i).AND.lc(i) < lb(i)).OR.izlim==1) subtria_n(i)=8
1256 izlim = -1
1257 ENDIF
1258 ENDIF
1259 IF(izlim == -1)cycle
1260 ENDIF
1261
1262 la_tm(i) = la(i)
1263 lb_tm(i) = lb(i)
1264 IF(la(i)<zero)THEN
1265 IF(lb(i)<zero)THEN
1266 la(i) = zero
1267 lb(i) = zero
1268 lc(i) = one
1269 ELSEIF(lc(i)<zero)THEN
1270 lc(i) = zero
1271 la(i) = zero
1272 lb(i) = one
1273 ELSE
1274 la(i) = zero
1275 aaa = lb(i) + lc(i)
1276 lb(i) = lb(i)/aaa
1277 lc(i) = lc(i)/aaa
1278 ENDIF
1279 ELSEIF(lb(i)<zero)THEN
1280 IF(lc(i)<zero)THEN
1281 lb(i) = zero
1282 lc(i) = zero
1283 la(i) = one
1284 ELSE
1285 lb(i) = zero
1286 aaa = lc(i) + la(i)
1287 lc(i) = lc(i)/aaa
1288 la(i) = la(i)/aaa
1289 ENDIF
1290 ELSEIF(lc(i)<zero)THEN
1291 lc(i) = zero
1292 aaa = la(i) + lb(i)
1293 la(i) = la(i)/aaa
1294 lb(i) = lb(i)/aaa
1295 ENDIF
1296
1297
1298
1299 IF(izlim==1 .OR. ispt2(i)>0)THEN
1300
1301 nm1(i) = n1(i)
1302 nm2(i) = n2(i)
1303 nm3(i) = n3(i)
1304 k = ic(1,itq)
1305 nax = nxx(i,k)
1306 nay = nyy(i,k)
1307 naz = nzz(i,k)
1308
1309 k = ic(2,itq)
1310 nbx = nxx(i,k)
1311 nby = nyy(i,k)
1312 nbz = nzz(i,k)
1313
1314 k = ic(3,itq)
1315 ncx = nxx(i,k)
1316 ncy = nyy(i,k)
1317 ncz = nzz(i,k)
1318 n1(i) = la(i)*nax + lb(i)*nbx + lc(i)*ncx
1319 n2(i) = la(i)*nay + lb(i)*nby + lc(i)*ncy
1320 n3(i) = la(i)*naz + lb(i)*nbz + lc(i)*ncz
1321
1322
1323 nni = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
1324 ni2 = n1(i)*n1(i) + n2(i)*n2(i) + n3(i)*n3(i)
1325 IF(two*nni*nni < ni2)THEN
1326
1327 aaa = sqrt(ni2-nni*nni) - nni
1328 n1(i) = n1(i) + aaa*nqx
1329 n2(i) = n2(i) + aaa*nqy
1330 n3(i) = n3(i) + aaa*nqz
1331 ni2 = n1(i)*n1(i) + n2(i)*n2(i) + n3(i)*n3(i)
1332 ENDIF
1333 s2 = one/
max(em30,sqrt(ni2))
1334 n1(i) = n1(i)*s2
1335 n2(i) = n2(i)*s2
1336 n3(i) = n3(i)*s2
1337
1338 nni = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
1339 IF (nni<em20) THEN
1340 n1(i) = nqx
1341 n2(i) = nqy
1342 n3(i) = nqz
1343 END IF
1344
1345
1346
1347
1348 ENDIF
1349
1350 IF(pene(i) == zero)THEN
1351 icontact(i) = 0
1352 istep(i) = 0
1353 cycle
1354 ENDIF
1355
1356 ix1(i) = ixx(i,ic( 7,itq))
1357 ix2(i) = ixx(i,ic( 8,itq))
1358 ix3(i) = ixx(i,ic( 9,itq))
1359 ix4(i) = ixx(i,ic(10,itq))
1360 x1(i) = xx0(i,ic( 7,itq))
1361 x2(i) = xx0(i,ic( 8,itq))
1362 x3(i) = xx0(i,ic( 9,itq))
1363 x4(i) = xx0(i,ic(10,itq))
1364 y1(i) = yy0(i,ic( 7,itq))
1365 y2(i) = yy0(i,ic( 8,itq))
1366 y3(i) = yy0(i,ic( 9,itq))
1367 y4(i) = yy0(i,ic(10,itq))
1368 z1(i) = zz0(i,ic( 7,itq))
1369 z2(i) = zz0(i,ic( 8,itq))
1370 z3(i) = zz0(i,ic( 9,itq))
1371 z4(i) = zz0(i,ic(10,itq))
1372
1373 vx1(i) = vx(i,ic( 7,itq))
1374 vx2(i) = vx(i,ic( 8,itq))
1375 vx3(i) = vx(i,ic( 9,itq))
1376 vx4(i) = vx(i,ic(10,itq))
1377 vy1(i) = vy(i,ic( 7,itq))
1378 vy2(i) = vy(i,ic( 8,itq))
1379 vy3(i) = vy(i,ic( 9,itq))
1380 vy4(i) = vy(i,ic(10,itq))
1381 vz1(i) = vz(i,ic( 7,itq))
1382 vz2(i) = vz(i,ic( 8,itq))
1383 vz3(i) = vz(i,ic( 9,itq))
1384 vz4(i) = vz(i,ic(10,itq))
1385
1386 IF (ix1(i) == ix2(i))THEN
1387 h1(i) = la(i)
1388 h2(i) = zero
1389 h3(i) = lb(i)
1390 h4(i) = lc(i)
1391 ELSEIF(ix2(i) == ix3(i))THEN
1392 h1(i) = la(i)
1393 h2(i) = lb(i)
1394 h3(i) = zero
1395 h4(i) = lc(i)
1396
1397 ELSEIF(ix4(i) == ix1(i))THEN
1398 h1(i) = lc(i)
1399 h2(i) = la(i)
1400 h3(i) = lb(i)
1401 h4(i) = zero
1402 ELSE
1403 h0 = fourth * la(i)
1404 h1(i) = h0
1405 h2(i) = h0
1406 h3(i) = h0 + lb(i)
1407 h4(i) = h0 + lc(i)
1408 ENDIF
1409 ENDDO
1410
1411
1412
1413
1414 DO i=1,jlt
1415 IF(istep(i) /= 2)cycle
1416
1417
1418 xi5 = xx(i,5) - xi(i)
1419 yi5 = yy(i,5) - yi(i)
1420 zi5 = zz(i,5) - zi(i)
1421
1422 v12 = sx125(i)*xi5 + sy125(i)*yi5 + sz125(i)*zi5
1423 v23 = sx235(i)*xi5 + sy235(i)*yi5 + sz235(i)*zi5
1424 v34 = sx345(i)*xi5 + sy345(i)*yi5 + sz345(i)*zi5
1425 v41 = sx415(i)*xi5 + sy415(i)*yi5 + sz415(i)*zi5
1426
1427 IF(v12 < zero .and. v23 < zero .and.
1428 . v34 < zero .and. v41 < zero ) THEN
1429
1430 cycle
1431 ENDIF
1432
1433 n=cand_n(i)
1434
1435
1436
1437
1438
1439
1440 IF(intply > 0) THEN
1441
1442 dgap = dgap_nm(1,cand_e(i))
1443 nx1 = nxx(i,1)
1444 ny1 = nyy(i,1)
1445 nz1 = nzz(i,1)
1446 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1447 nx1 = nx1 * aaa
1448 ny1 = ny1 * aaa
1449 nz1 = nz1 * aaa
1450
1451 xo1 = xx(i,1) - dt1*vx(i,1) - dgap*nx1
1452 yo1 = yy(i,1) - dt1*vy(i,1) - dgap*ny1
1453 zo1 = zz(i,1) - dt1*vz(i,1) - dgap*nz1
1454
1455 nx1 = nxx(i,2)
1456 ny1 = nyy(i,2)
1457 nz1 = nzz(i,2)
1458 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1459 nx1 = nx1 * aaa
1460 ny1 = ny1 * aaa
1461 nz1 = nz1 * aaa
1462
1463 xo2 = xx(i,2) - dt1*vx(i,2)- dgap*nx1
1464 yo2 = yy(i,2) - dt1*vy(i,2)- dgap*ny1
1465 zo2 = zz(i,2) - dt1*vz(i,2)- dgap*nz1
1466 dgap = dgap_nm(3,cand_e(i))
1467 nx1 = nxx(i,3)
1468 ny1 = nyy(i,3)
1469 nz1 = nzz(i,3)
1470 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1471 nx1 = nx1 * aaa
1472 ny1 = ny1 * aaa
1473 nz1 = nz1 * aaa
1474
1475 xo3 = xx(i,3) - dt1*vx(i,3)- dgap*nx1
1476 yo3 = yy(i,3) - dt1*vy(i,3)- dgap*ny1
1477 zo3 = zz(i,3) - dt1*vz(i,3)- dgap*nz1
1478 dgap = dgap_nm(4,cand_e(i))
1479 nx1 = nxx(i,4)
1480 ny1 = nyy(i,4)
1481 nz1 = nzz(i,4)
1482 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1483 nx1 = nx1 * aaa
1484 ny1 = ny1 * aaa
1485 nz1 = nz1 * aaa
1486
1487 xo4 = xx(i,4) - dt1*vx(i,4)- dgap*nx1
1488 yo4 = yy(i,4) - dt1*vy(i,4)- dgap*ny1
1489 zo4 = zz(i,4) - dt1*vz(i,4)- dgap*nz1
1490 ELSE
1491 xo1 = xx(i,1) - dt1*vx(i,1)
1492 yo1 = yy(i,1) - dt1*vy(i,1)
1493 zo1 = zz(i,1) - dt1*vz(i,1)
1494
1495 xo2 = xx(i,2) - dt1*vx(i,2)
1496 yo2 = yy(i,2) - dt1*vy(i,2)
1497 zo2 = zz(i,2) - dt1*vz(i,2)
1498
1499 xo3 = xx(i,3) - dt1*vx(i,3)
1500 yo3 = yy(i,3) - dt1*vy(i,3)
1501 zo3 = zz(i,3) - dt1*vz(i,3)
1502
1503 xo4 = xx(i,4) - dt1*vx(i,4)
1504 yo4 = yy(i,4) - dt1*vy(i,4)
1505 zo4 = zz(i,4) - dt1*vz(i,4)
1506 ENDIF
1507
1508 xoi = xi(i) - dt1*vxi(i)
1509 yoi = yi(i) - dt1*vyi(i)
1510 zoi = zi(i) - dt1*vzi(i)
1511
1512 IF(ixx(i,3) /= ixx(i,4))THEN
1513 xo5 = fourth*(xo1+xo2+xo3+xo4)
1514 yo5 = fourth*(yo1+yo2+yo3+yo4)
1515 zo5 = fourth*(zo1+zo2+zo3+zo4)
1516 ELSE
1517 xo5 = xo3
1518 yo5 = yo3
1519 zo5 = zo3
1520 ENDIF
1521
1522
1523
1524 x51 = xo1 - xo5
1525 y51 = yo1 - yo5
1526 z51 = zo1 - zo5
1527
1528 x52 = xo2 - xo5
1529 y52 = yo2 - yo5
1530 z52 = zo2 - zo5
1531
1532 x53 = xo3 - xo5
1533 y53 = yo3 - yo5
1534 z53 = zo3 - zo5
1535
1536 x54 = xo4 - xo5
1537 y54 = yo4 - yo5
1538 z54 = zo4 - zo5
1539
1540 xi5 = xo5 - xoi
1541 yi5 = yo5 - yoi
1542 zi5 = zo5 - zoi
1543
1544 sox125 = y51*z52 - z51*y52
1545 soy125 = z51*x52 - x51*z52
1546 soz125 = x51*y52 - y51*x52
1547 vo12 = sox125*xi5 + soy125*yi5 + soz125*zi5
1548
1549 sox235 = y52*z53 - z52*y53
1550 soy235 = z52*x53 - x52*z53
1551 soz235 = x52*y53 - y52*x53
1552 vo23 = sox235*xi5 + soy235*yi5 + soz235*zi5
1553
1554 sox345 = y53*z54 - z53*y54
1555 soy345 = z53*x54 - x53*z54
1556 soz345 = x53*y54 - y53*x54
1557 vo34 = sox345*xi5 + soy345*yi5 + soz345*zi5
1558
1559 sox415 = y54*z51 - z54*y51
1560 soy415 = z54*x51 - x54*z51
1561 soz415 = x54*y51 - y54*x51
1562 vo41 = sox415*xi5 + soy415*yi5 + soz415*zi5
1563
1564 adt0(i) = -one
1565
1566
1567
1568
1569
1570 subtria(i)= 0
1571
1574 1 istep(i) ,subtria(i),ixx(i,3) ,ixx(i,4),
1576 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
1577 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2) ,
1578 4 xx(i,3 ) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4) ,
1579 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),adt0(i) ,
1580 6 vxi(i) ,vyi(i) ,vzi(i) ,xo1 ,xo2 ,
1581 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1582 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1583 9 zo3 ,zo4 ,zo5 ,nx(i) ,ny(i) ,
1584 a nz(i) ,sx125(i),sx235(i),sx345(i),sx415(i),
1585 b sy125(i) ,sy235(i),sy345(i),sy415(i),sz125(i),
1586 c sz235(i) ,sz345(i),sz415(i),xi(i) ,yi(i) ,
1587 d zi(i) ,zerom ,unp ,zeromt,unpt )
1588 IF(
intersect == 0.AND.ishel(i)>1.AND.ispt2(i)>0)
THEN
1589 xi5 = xx0(i,5) - xi(i)
1590 yi5 = yy0(i,5) - yi(i)
1591 zi5 = zz0(i,5) - zi(i)
1592 aaa = abs(sx1250(i)*xi5 + sy1250(i)*yi5 + sz1250(i)*zi5)
1593 aaa =
max(aaa,(sx2350(i)*xi5 + sy2350(i)*yi5 + sz2350(i)*zi5))
1594 aaa =
max(aaa,(sx3450(i)*xi5 + sy3450(i)*yi5 + sz3450(i)*zi5
1595 aaa =
max(aaa,(sx4150(i)*xi5 + sy4150(i)*yi5 + sz4150(i)*zi5)
1596 IF (aaa<five*em04)
1598 1 istep(i) ,subtria(i),ixx(i,3) ,ixx(i,4),
1600 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2) ,
1601 4 xx(i,3 ) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4) ,
1602 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),adt0(i) ,
1603 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
1604 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1605 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1606 9 zo3 ,zo4 ,zo5 ,nx(i) ,ny(i) ,
1607 a nz(i) ,sx125(i),sx235(i),sx345(i),sx415(i),
1608 b sy125(i) ,sy235(i),sy345(i),sy415(i),sz125(i),
1609 c sz235(i) ,sz345(i),sz415(i),xi(i) ,yi(i) ,
1610 d zi(i) ,sox125 ,sox235 ,sox345 ,sox415 ,
1611 b soy125 ,soy235 ,soy345 ,soy415 ,soz125 ,
1612 c soz235 ,soz345 ,soz415 ,mvoisin(1,cand_e(i)) )
1613 END IF
1614 ikeep = 1
1616 IF (icont_i(n) < 0 ) ikeep = 0
1619 ENDIF
1620
1621 IF(ikeep == 1.AND.ishel(i)>0.AND.ispt2(i)==0.AND.inacti/=5)THEN
1622 IF ((gaps(i)+gap_m(cand_e(i)))>tol_ints) ikeep = 0
1623 END IF
1626 1 istep(i) ,subtria(i),ixx(i,3) ,ixx(i,4),
1628 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
1629 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2) ,
1630 4 xx(i,3 ) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4) ,
1631 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),adt0(i) ,
1632 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
1633 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1634 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1635 9 zo3 ,zo4 ,zo5 ,nx(i) ,ny(i) ,
1636 a nz(i) ,sx125(i),sx235(i),sx345(i),sx415(i),
1637 b sy125(i) ,sy235(i),sy345(i),sy415(i),sz125(i),
1638 c sz235(i) ,sz345(i),sz415(i),xi(i) ,yi(i) ,
1639 d zi(i) ,zerom ,unp )
1640
1641 IF (tt==zero.AND.
intersect>0.AND.inacti==0)
THEN
1642 IF(n <= nsn)THEN
1643 icont_i(n) = -cand_e(i)
1644 ELSE
1646 ENDIF
1648 istep(i) = 2
1649 END IF
1652
1653 ENDDO
1654
1655
1656
1657
1658
1659
1660
1661
1662
1663
1664
1665
1666
1667
1668
1669
1670
1671
1672
1673
1674
1675
1676
1677
1678
1679
1680
1681
1682
1683
1684
1685
1686
1687
1688
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712
1713
1714
1715
1716
1717
1718
1719
1720
1721
1722
1723
1724
1725
1726
1727
1728
1729
1730
1731
1732
1733
1734
1735
1736
1737
1738
1739
1740
1741
1742
1743
1744
1745
1746
1747
1748
1749
1750
1751
1752
1753
1754
1755 DO i=1,jlt
1756
1757 IF(istep(i) /= 3)cycle
1758
1759
1760 nxx(i,6) =sx21140(i)
1761 nyy(i,6) =sy21140(i)
1762 nzz(i,6) =sz21140(i)
1763
1764 nxx(i,7) =sx21140(i)
1765 nyy(i,7) =sy21140(i)
1766 nzz(i,7) =sz21140(i)
1767
1768 nxx(i,8) =sx32150(i)
1769 nyy(i,8) =sy32150(i)
1770 nzz(i,8) =sz32150(i)
1771
1772 nxx(i,9) =sx32150(i)
1773 nyy(i,9) =sy32150(i)
1774 nzz(i,9) =sz32150(i)
1775
1776 nxx(i,15)=sx32150(i)
1777 nyy(i,15)=sy32150(i)
1778 nzz(i,15)=sz32150(i)
1779
1780 nxx(i,10)=sx43160(i)
1781 nyy(i,10)=sy43160(i)
1782 nzz(i,10)=sz43160(i)
1783
1784 nxx(i,11)=sx43160(i)
1785 nyy(i,11)=sy43160(i)
1786 nzz(i,11)=sz43160(i)
1787
1788 nxx(i,14)=sx21140(i)
1789 nyy(i,14)=sy21140(i)
1790 nzz(i,14)=sz21140(i)
1791
1792 nxx(i,16)=sx43160(i)
1793 nyy(i,16)=sy43160(i)
1794 nzz(i,16)=sz43160(i)
1795
1796 nxx(i,12)=sx14170(i)
1797 nyy(i,12)=sy14170(i)
1798 nzz(i,12)=sz14170(i)
1799
1800 nxx(i,13)=sx14170(i)
1801 nyy(i,13)=sy14170(i)
1802 nzz(i,13)=sz14170(i)
1803
1804 nxx(i,17)=sx14170(i)
1805 nyy(i,17)=sy14170(i)
1806 nzz(i,17)=sz14170(i)
1807
1808
1809
1810 n=cand_n(i)
1811
1812
1813
1814
1815
1816
1817 IF(ishel(i)>0) THEN
1818
1819
1820
1821
1822 bbb = sqrt(sx21140(i)*sx21140(i)
1823 + +sy21140(i)*sy21140(i)
1824 + +sz21140(i)*sz21140(i))
1825 aaa = gaps(i)+gap_nm(5,cand_e(i))
1826 aaa = aaa/bbb
1827 xx(i,6) = xx0(i,6) + sx21140(i) * aaa
1828 yy(i,6) = yy0(i,6) + sy21140(i) * aaa
1829 zz(i,6) = zz0(i,6) + sz21140(i) * aaa
1830 aaa = gaps(i)+gap_nm(6,cand_e(i))
1831 aaa = aaa/bbb
1832 xx(i,7) = xx0(i,7) + sx21140(i) * aaa
1833 yy(i,7) = yy0(i,7) + sy21140(i) * aaa
1834 zz(i,7) = zz0(i,7) + sz21140(i) * aaa
1835
1836 bbb = sqrt(sx32150(i)*sx32150(i)
1837 + +sy32150(i)*sy32150(i)
1838 + +sz32150(i)*sz32150(i))
1839 aaa = gaps(i)+gap_nm(7,cand_e(i))
1840 aaa = aaa/bbb
1841 xx(i,8) = xx0(i,8) + sx32150(i) * aaa
1842 yy(i,8) = yy0(i,8) + sy32150(i) * aaa
1843 zz(i,8) = zz0(i,8) + sz32150(i) * aaa
1844 aaa = gaps(i)+gap_nm(8,cand_e(i))
1845 aaa = aaa/bbb
1846 xx(i,9) = xx0(i,9) + sx32150(i) * aaa
1847 yy(i,9) = yy0(i,9) + sy32150(i) * aaa
1848 zz(i,9) = zz0(i,9) + sz32150(i) * aaa
1849
1850 IF(ixx(i,3) == ixx(i,4))THEN
1851
1852 xx(i,10) = xx0(i,10)
1853 yy(i,10) = yy0(i,10)
1854 zz(i,10) = zz0(i,10)
1855
1856 xx(i,11) = xx0(i,11)
1857 yy(i,11) = yy0(i,11)
1858 zz(i,11) = zz0(i,11)
1859
1860 ELSE
1861
1862 bbb = sqrt(sx43160(i)*sx43160(i)
1863 + +sy43160(i)*sy43160(i)
1864 + +sz43160(i)*sz43160(i))
1865 aaa = gaps(i)+gap_nm(9,cand_e(i))
1866 aaa = aaa/bbb
1867 xx(i,10) = xx0(i,10) + sx43160(i) * aaa
1868 yy(i,10) = yy0(i,10) + sy43160(i) * aaa
1869 zz(i,10) = zz0(i,10) + sz43160(i) * aaa
1870 aaa = gaps(i)+gap_nm(10,cand_e(i))
1871 aaa = aaa/bbb
1872 xx(i,11) = xx0(i,11) + sx43160(i) * aaa
1873 yy(i,11) = yy0(i,11) + sy43160(i) * aaa
1874 zz(i,11) = zz0(i,11) + sz43160(i) * aaa
1875 ENDIF
1876
1877 bbb = sqrt(sx14170(i)*sx14170(i)
1878 + +sy14170(i)*sy14170(i)
1879 + +sz14170(i)*sz14170(i))
1880 aaa = gaps(i)+gap_nm(11,cand_e(i))
1881 aaa = aaa/bbb
1882 xx(i,12) = xx0(i,12) + sx14170(i) * aaa
1883 yy(i,12) = yy0(i,12) + sy14170(i) * aaa
1884 zz(i,12) = zz0(i,12) + sz14170(i) * aaa
1885 aaa = gaps(i)+gap_nm(12,cand_e(i))
1886 aaa = aaa/bbb
1887 xx(i,13) = xx0(i,13) + sx14170(i) * aaa
1888 yy(i,13) = yy0(i,13) + sy14170(i) * aaa
1889 zz(i,13) = zz0(i,13) + sz14170(i) * aaa
1890
1891
1892 IF(ixx(i,6) == ixx(i,7))THEN
1893 xx(i,14) = xx(i,6)
1894 yy(i,14) = yy(i,6)
1895 zz(i,14) = zz(i,6)
1896 ELSE
1897 xx(i,14) = fourth*(xx(i,2)+xx(i,1)+xx(i,6)+xx(i,7))
1898 yy(i,14) = fourth*(yy(i,2)+yy(i,1)+yy(i,6)+yy(i,7))
1899 zz(i,14) = fourth*(zz(i,2)+zz(i,1)+zz(i,6)+zz(i,7))
1900 ENDIF
1901 IF(ixx(i,8)==ixx(i,9))THEN
1902 xx(i,15) = xx(i,8)
1903 yy(i,15) = yy(i,8)
1904 zz(i,15) = zz(i,8)
1905 ELSE
1906 xx(i,15) = fourth*(xx(i,3)+xx(i,2)+xx(i,8)+xx(i,9))
1907 yy(i,15) = fourth*(yy(i,3)+yy(i,2)+yy(i,8)+yy(i,9))
1908 zz(i,15) = fourth*(zz(i,3)+zz(i,2)+zz(i,8)+zz(i,9))
1909 ENDIF
1910 IF(ixx(i,10)==ixx(i,11))THEN
1911 xx(i,16) = xx(i,10)
1912 yy(i,16) = yy(i,10)
1913 zz(i,16) = zz(i,10)
1914 ELSE
1915 xx(i,16) = fourth*(xx(i,4)+xx(i,3)+xx(i,10)+xx(i,11))
1916 yy(i,16) = fourth*(yy(i,4)+yy(i,3)+yy(i,10)+yy(i,11))
1917 zz(i,16) = fourth*(zz(i,4)+zz(i,3)+zz(i,10)+zz(i,11))
1918 ENDIF
1919 IF(ixx(i,12)==ixx(i,13))THEN
1920 xx(i,17) = xx(i,12)
1921 yy(i,17) = yy(i,12)
1922 zz(i,17) = zz(i,12)
1923 ELSE
1924 xx(i,17) = fourth*(xx(i,1)+xx(i,4)+xx(i,12)+xx(i,13))
1925 yy(i,17) = fourth*(yy(i,1)+yy(i,4)+yy(i,12)+yy(i,13))
1926 zz(i,17) = fourth*(zz(i,1)+zz(i,4)+zz(i,12)+zz(i,13))
1927 ENDIF
1928
1929 x141 = xx(i,1) - xx(i,14)
1930 y141 = yy(i,1) - yy(i,14)
1931 z141 = zz(i,1) - zz(i,14)
1932
1933 x142 = xx(i,2) - xx(i,14)
1934 y142 = yy(i,2) - yy(i,14)
1935 z142 = zz(i,2) - zz(i,14)
1936
1937 x152 = xx(i,2) - xx(i,15)
1938 y152 = yy(i,2) - yy(i,15)
1939 z152 = zz(i,2) - zz(i,15)
1940
1941 x153 = xx(i,3) - xx(i,15)
1942 y153 = yy(i,3) - yy(i,15)
1943 z153 = zz(i,3) - zz(i,15)
1944
1945 x163 = xx(i,3) - xx(i,16)
1946 y163 = yy(i,3) - yy(i,16)
1947 z163 = zz(i,3) - zz(i,16)
1948
1949 x164 = xx(i,4) - xx(i,16)
1950 y164 = yy(i,4) - yy(i,16)
1951 z164 = zz(i,4) - zz(i,16)
1952
1953 x174 = xx(i,4) - xx(i,17)
1954 y174 = yy(i,4) - yy(i,17)
1955 z174 = zz(i,4) - zz(i,17)
1956
1957 x171 = xx(i,1) - xx(i,17)
1958 y171 = yy(i,1) - yy(i,17)
1959 z171 = zz(i,1) - zz(i,17)
1960
1961 IF(mvoisin(1,cand_e(i))/=0)THEN
1962 sx2114(i) = y142*z141 - z142*y141
1963 sy2114(i) = z142*x141 - x142*z141
1964 sz2114(i) = x142*y141 - y142*x141
1965 ELSE
1966 sx2114(i) = sx125(i)
1967 sy2114(i) = sy125(i)
1968 sz2114(i) = sz125(i)
1969 ENDIF
1970
1971 IF(mvoisin(2,cand_e(i))/=0)THEN
1972 sx3215(i) = y153*z152 - z153*y152
1973 sy3215(i) = z153*x152 - x153*z152
1974 sz3215(i) = x153*y152 - y153*x152
1975 ELSE
1976 sx3215(i) = sx235(i)
1977 sy3215(i) = sy235(i)
1978 sz3215(i) = sz235(i)
1979 ENDIF
1980
1981 IF(mvoisin(3,cand_e(i))/=0)THEN
1982 sx4316(i) = y164*z163 - z164*y163
1983 sy4316(i) = z164*x163 - x164*z163
1984 sz4316(i) = x164*y163 - y164*x163
1985 ELSE
1986 sx4316(i) = sx345(i)
1987 sy4316(i) = sy345(i)
1988 sz4316(i) = sz345(i)
1989 ENDIF
1990
1991 IF(mvoisin(4,cand_e(i))/=0)THEN
1992 sx1417(i) = y171*z174 - z171*y174
1993 sy1417(i) = z171*x174 - x171*z174
1994 sz1417(i) = x171*y174 - y171*x174
1995 ELSE
1996 sx1417(i) = sx415(i)
1997 sy1417(i) = sy415(i)
1998 sz1417(i) = sz415(i)
1999 ENDIF
2000
2001 ELSE
2002
2003 xx(i,6:17) = xx0(i,6:17)
2004 yy(i,6:17) = yy0(i,6:17)
2005 zz(i,6:17) = zz0(i,6:17)
2006
2007 sx2114(i) = sx21140(i)
2008 sy2114(i) = sy21140(i)
2009 sz2114(i) = sz21140(i)
2010
2011 sx3215(i) = sx32150(i)
2012 sy3215(i) = sy32150(i)
2013 sz3215(i) = sz32150(i)
2014
2015 sx4316(i) = sx43160(i)
2016 sy4316(i) = sy43160(i)
2017 sz4316(i) = sz43160(i)
2018
2019 sx1417(i) = sx14170(i)
2020 sy1417(i) = sy14170(i)
2021 sz1417(i) = sz14170(i)
2022
2023 ENDIF
2024
2025 ENDDO
2026
2027
2028
2029
2030
2031
2032
2033 nslid = 0
2034 DO i=1,jlt
2035 IF(istep(i) /= 3)cycle
2036 istep(i) = 5
2037
2038
2039
2040
2041
2042
2043
2044
2045! | 2 | 6 3 1 | 6 3 1 | - - - | - - - |
2046
2047
2048! +------+----------+----------+----------+----------+
2049
2050
2051
2052
2053
2054
2055
2056
2057
2058
2059
2060
2061
2062
2063
2064
2065
2066
2067
2068
2069 it(1:3,1:20,i)=it0(1:3,1:20)
2070
2071 IF(ixx(i,3) /= ixx(i,4))THEN
2072 IF(mvoisin(1,cand_e(i))==0)it( 1, 1,i) = 0
2073 IF(mvoisin(2,cand_e(i))==0)it( 1, 2,i) = 0
2074 IF(mvoisin(3,cand_e(i))==0)it( 1, 3,i) = 0
2075 IF(mvoisin(4,cand_e(i))==0)it( 1, 4,i) = 0
2076 ELSE
2077 it(2,1,i)=6
2078 it(3,1,i)=8
2079 it(1,5,i)=1
2080 it(1,6,i)=1
2081 it(1,8,i)=1
2082 IF(mvoisin(1,cand_e(i))==0)it( 1, 1,i) = 0
2083 IF(mvoisin(2,cand_e(i))==0)it( 2, 1,i) = 0
2084 IF(mvoisin(3,cand_e(i))==0)it( 3, 1,i) = 0
2085
2086 ENDIF
2087 IF(ixx(i,6)==ixx(i,7))THEN
2088 it(2,5,i)=-1
2089 it(3,5,i)=-1
2090 ENDIF
2091 IF(ixx(i,8)==ixx(i,9))THEN
2092 it(2,6,i)=-1
2093 it(3,6,i)=-1
2094 ENDIF
2095 IF(ixx(i,10)==ixxTHEN
2096
2097 it(3,7,i)=-1
2098 ENDIF
2099 IF(ixx(i,12)==ixx(i,13))THEN
2100 it(2,8,i)=-1
2101 it(3,8,i)=-1
2102 ENDIF
2103
2104 IF(ixx(i,3) /= ixx(i,4))THEN
2105
2106 IF(subtria(i)==1)THEN
2107 xa(i) = xx(i,5)
2108 ya(i) = yy(i,5)
2109 za(i) = zz(i,5)
2110 xb(i) = xx(i,1)
2111 yb(i) = yy(i,1)
2112 zb(i) = zz(i,1)
2113 xc(i) = xx(i,2)
2114 yc(i) = yy(i,2)
2115 zc(i) = zz(i,2)
2116
2117 nqx = sx125(i)
2118 nqy = sy125(i) ! \ /
2119 nqz = sz125(i)
2120 nsx = sx235(i)
2121 nsy = sy235(i)
2122 nsz = sz235(i)
2123 ntx = sx415(i)
2124 nty = sy415(i)
2125 ntz = sz415(i)
2126 nrx = sx2114(i)
2127 nry = sy2114(i)
2128 nrz = sz2114(i)
2129 itr = 5
2130 IF(mvoisin(1,cand_e(i))==0)THEN
2131 itr = -itr
2132 nrx = nqx
2133 nry = nqy
2134 nrz = nqz
2135 ENDIF
2136 its = 2
2137 itt = 4
2138 ELSEIF(subtria(i)==2)THEN
2139 xa(i) = xx(i,5)
2140 ya(i) = yy(i,5)
2141 za(i) = zz
2142 xb(i) = xx(i,2)
2143 yb(i) = yy(i,2)
2144 zb(i) = zz(i,2)
2145 xc(i) = xx(i,3)
2146 yc(i) = yy(i,3)
2147 zc(i) = zz(i,3)
2148 nqx = sx235(i)
2149 nqy = sy235(i)
2150 nqz = sz235(i)
2151 nsx = sx345(i)
2152 nsy = sy345(i)
2153 nsz = sz345(i)
2154 ntx = sx125(i)
2155 nty = sy125(i)
2156 ntz = sz125(i)
2157 nrx = sx3215(i)
2158 nry = sy3215(i)
2159 nrz = sz3215(i)
2160 itr = 6
2161 IF(mvoisin(2,cand_e(i))==0)THEN
2162 itr = -itr
2163 nrx = nqx
2164 nry = nqy
2165 nrz = nqz
2166 ENDIF
2167 its = 3
2168 itt = 1
2169 ELSEIF(subtria(i)==3)THEN
2170 xa(i) = xx(i,5)
2171 ya(i) = yy(i,5)
2172 za(i) = zz(i,5)
2173 xb(i) = xx(i,3)
2174 yb(i) = yy(i,3)
2175 zb(i) = zz(i,3)
2176 xc(i) = xx(i,4)
2177 yc(i) = yy(i,4)
2178 zc(i) = zz(i,4)
2179 nqx = sx345(i)
2180 nqy = sy345(i)
2181 nqz = sz345(i)
2182 nsx = sx415(i)
2183 nsy = sy415(i)
2184 nsz = sz415(i)
2185 ntx = sx235(i)
2186 nty = sy235(i)
2187 ntz = sz235(i)
2188 nrx = sx4316(i)
2189 nry = sy4316(i)
2190 nrz = sz4316(i)
2191 itr = 7
2192 IF(mvoisin(3,cand_e(i))==0)THEN
2193 itr = -itr
2194 nrx = nqx
2195 nry = nqy
2196 nrz = nqz
2197 ENDIF
2198 its = 4
2199 itt = 2
2200 ELSEIF(subtria(i)==4)THEN
2201 xa(i) = xx(i,5)
2202 ya(i) = yy(i,5)
2203 za(i) = zz(i,5)
2204 xb(i) = xx(i,4)
2205 yb(i) = yy(i,4)
2206 zb(i) = zz(i,4)
2207 xc(i) = xx(i,1)
2208 yc(i) = yy(i,1)
2209 zc(i) = zz(i,1)
2210 nqx = sx415(i)
2211 nqy = sy415(i)
2212 nqz = sz415(i)
2213 nsx = sx125(i)
2214 nsy = sy125(i)
2215 nsz = sz125(i)
2216 ntx = sx345(i)
2217 nty = sy345(i)
2218 ntz = sz345(i)
2219 nrx = sx1417(i)
2220 nry = sy1417(i)
2221 nrz = sz1417(i)
2222 itr = 8
2223 IF(mvoisin(4,cand_e(i))==0)THEN
2224 itr = -itr
2225 nrx = nqx
2226 nry = nqy
2227 nrz = nqz
2228 ENDIF
2229 its = 1
2230 itt = 3
2231 ENDIF
2232 else
2233 xa(i) = xx(i,3)
2234 ya(i) = yy(i,3)
2235 za(i) = zz(i,3)
2236 xb(i) = xx(i,1)
2237 yb(i) = yy(i,1)
2238 zb(i) = zz(i,1)
2239 xc(i) = xx(i,2)
2240 yc(i) = yy(i,2)
2241 zc(i) = zz(i,2)
2242 nqx = sx125(i)
2243 nqy = sy125(i)
2244 nqz = sz125(i)
2245 nsx = sx3215(i)
2246 nsy = sy3215(i)
2247 nsz = sz3215(i)
2248 ntx = sx1417(i)
2249 nty = sy1417(i)
2250 ntz = sz1417(i)
2251 nrx = sx2114(i)
2252 nry = sy2114(i)
2253 nrz = sz2114(i)
2254 itr = 5
2255 IF(mvoisin(1,cand_e(i))==0)THEN
2256 itr = -itr
2257 nrx = nqx
2258 nry = nqy
2259 nrz = nqz
2260 ENDIF
2261 its = 6
2262 IF(mvoisin(2,cand_e(i))==0)THEN
2263 its = -its
2264 nsx = nqx
2265 nsy = nqy
2266 nsz = nqz
2267 ENDIF
2268 itt = 8
2269 IF(mvoisin(4,cand_e(i))==0)THEN
2270 itt = -itt
2271 ntx = nqx
2272 nty = nqy
2273 ntz = nqz
2274 ENDIF
2275 ENDIF
2276
2277 nx(i) = nqx
2278 ny(i) = nqy
2279 nz(i) = nqz
2280
2281
2282 nqrx = nrx + nqx
2283 nqry = nry + nqy
2284 nqrz = nrz + nqz
2285 nqsx = nsx + nqx
2286 nqsy = nsy + nqy
2287 nqsz = nsz + nqz
2288 nqtx = ntx + nqx
2289 nqty = nty + nqy
2290 nqtz = ntz + nqz
2291 IF(abs(itr)==subtria_n(i))THEN
2292 prx = nqry*(zc(i)-zb(i)) - nqrz*(yc(i)-yb(i))
2293 pry = nqrz*(xc(i)-xb(i)) - nqrx*(zc(i)-zb(i))
2294 prz = nqrx*(yc(i)-yb(i)) - nqry*(xc(i)-xb(i))
2295 vr = prx*(xi(i)-xb(i))+ pry*(yi(i)-yb(i))+ prz*(zi(i)-zb(i))
2296 IF(vr < zero)THEN
2297 IF(itr>0)THEN
2298 istep(i) = 5
2299 nslid = nslid+1
2300 subtria(i) = itr
2301 ELSE
2302
2303 icontact(i) = 0
2304 istep(i) = 0
2305 subtria(i) = 0
2306 ENDIF
2307 nx(i) = nrx
2308 ny(i) = nry
2309 nz(i) = nrz
2310 ENDIF
2311 ELSEIF(abs(its)==subtria_n(i))THEN
2312 psx = nqsy*(za(i)-zc(i)) - nqsz*(ya(i)-yc(i))
2313 psy = nqsz*(xa(i)-xc(i)) - nqsx*(za(i)-zc(i))
2314 psz = nqsx*(ya(i)-yc(i)) - nqsy*(xa(i)-xc(i))
2315 vs = psx*(xi(i)-xc(i))+ psy*(yi(i)-yc(i))+ psz*(zi(i)-zc(i))
2316 IF(vs < zero)THEN
2317 IF(its>0)THEN
2318 istep(i) = 5
2319 nslid = nslid+1
2320 subtria(i) = its
2321 ELSE
2322
2323 icontact(i) = 0
2324 istep(i) = 0
2325 subtria(i) = 0
2326 ENDIF
2327 nx(i) = nsx
2328 ny(i) = nsy
2329 nz(i) = nsz
2330 ENDIF
2331 ELSEIF(abs(itt)==subtria_n(i))THEN
2332 ptx = nqty*(zb(i)-za(i)) - nqtz*(yb(i)-ya(i))
2333 pty = nqtz*(xb(i)-xa(i)) - nqtx*(zb(i)-za(i))
2334 ptz = nqtx*(yb(i)-ya(i)) - nqty*(xb(i)-xa(i))
2335 vt = ptx*(xi(i)-xa(i))+ pty*(yi(i)-ya(i))+ ptz*(zi(i)-za(i))
2336 IF(vt < zero)THEN
2337 IF(itt>0)THEN
2338 istep(i) = 5
2339 nslid = nslid+1
2340 subtria(i) = itt
2341 ELSE
2342
2343 icontact(i) = 0
2344 istep(i) = 0
2345 subtria(i) = 0
2346 ENDIF
2347 nx(i) = ntx
2348 ny(i) = nty
2349 nz(i) = ntz
2350 ENDIF
2351 ENDIF
2352 ENDDO
2353
2354
2355
2356
2357 nss = 1
2358
2359 nslid=0
2360
2361
2362
2363
2364#include "lockon.inc"
2365
2366 DO i=1,jlt
2367 IF(istep(i) /= 1)pene(i) = zero
2368 IF(istep(i) < 4 .OR.istep(i) == 6)cycle
2369 itq = subtria(i)
2370 k = ic(1,itq)
2371 xa(i) = xx(i,k)
2372 ya(i) = yy(i,k)
2373 za(i) = zz(i,k)
2374 nax = nxx(i,k)
2375 nay = nyy(i,k)
2376 naz = nzz(i,k)
2377
2378 k = ic(2,itq)
2379 xb(i) = xx(i,k)
2380 yb(i) = yy(i,k)
2381 zb(i) = zz(i,k)
2382 nbx = nxx(i,k)
2383 nby = nyy(i,k)
2384 nbz = nzz(i,k)
2385
2386 k = ic(3,itq)
2387 xc(i) = xx(i,k)
2388 yc(i) = yy(i,k)
2389 zc(i) = zz(i,k)
2390 ncx = nxx(i,k)
2391 ncy = nyy(i,k)
2392 ncz = nzz(i,k)
2393
2394 nqx = nx(i)
2395 nqy = ny(i)
2396 nqz = nz(i)
2397
2398 aaa =
max(em30,sqrt(nqx*nqx + nqy*nqy + nqz*nqz))
2400 s2 = one/aaa
2401 nqx = nqx*s2
2402 nqy = nqy*s2
2403 nqz = nqz*s2
2404
2405 bbb = (xi(i)-xa(i))*nqx+(yi(i)-ya(i))*nqy+(zi(i)-za(i))*nqz
2406 ikeep=0
2407
2408
2409 IF (bbb>=zero .AND.bbb<marge025.AND.eps>em02 .AND. impl_s==0
2410 + .AND. nsne==0) THEN
2411 pene(i) = em15
2412 ikeep=1
2413 ELSE
2414 pene(i) = -
min(zero,bbb)
2415 END IF
2416
2417
2418
2419 IF(icont_r(i) >0)THEN
2420 tole =eps02*aaa
2421 gap2 = gaps(i)+gap_m(cand_e(i))
2422 IF (gap2 >zero) tole =
min(tole,gap2*gap2)
2423 IF(pene(i)*pene(i) > tole) THEN
2424 pene(i) = zero
2425 END IF
2427
2428 IF(pene(i) == zero)THEN
2429 icontact(i) = 0
2430 cycle
2431 ENDIF
2432
2433
2434
2435
2436
2437
2438
2439
2440
2441
2442
2443
2444
2445
2446 ! / / \ \ normal at point a:
2447
2448
2449
2450 ! / / \ \
2451
2452
2453
2454
2455
2456 nni = (nqx*nax + nqy*nay + nqz*naz)
2457 ni2 = nax*nax + nay*nay + naz*naz
2458 IF(two*nni*nni < ni2)THEN
2459
2460 aaa = sqrt(ni2-nni*nni) - nni
2461 nax = nax + aaa*nqx
2462 nay = nay + aaa*nqy
2463 naz = naz + aaa*nqz
2464 nni = (nqx*nax + nqy*nay + nqz*naz)
2465 ENDIF
2466
2467 aaa = pene(i)/nni
2468 nax = nax*aaa
2469 nay = nay*aaa
2470 naz = naz*aaa
2471
2472 xpa = xa(i) - nax
2473 ypa = ya(i) - nay
2474 zpa = za(i) - naz
2475
2476 nni = (nqx*nbx + nqy*nby + nqz*nbz)
2477 ni2 = nbx*nbx + nby*nby + nbz*nbz
2478 IF(two*nni*nni < ni2)THEN
2479
2480 aaa = sqrt(ni2-nni*nni) - nni
2481 nbx = nbx + aaa*nqx
2482 nby = nby + aaa*nqy
2483 nbz = nbz + aaa*nqz
2484 nni = (nqx*nbx + nqy*nby + nqz*nbz)
2485 ENDIF
2486
2487 aaa = pene(i)/nni
2488 nbx = nbx*aaa
2489 nby = nby*aaa
2490 nbz = nbz*aaa
2491
2492 xpb = xb(i) - nbx
2493 ypb = yb(i) - nby
2494 zpb = zb(i) - nbz
2495
2496 nni = (nqx*ncx + nqy*ncy + nqz*ncz)
2497 ni2 = ncx*ncx + ncy*ncy + ncz*ncz
2498 IF(two*nni*nni < ni2)THEN
2499
2500 aaa = sqrt(ni2-nni*nni) - nni
2501 ncx = ncx + aaa*nqx
2502 ncy = ncy + aaa*nqy
2503 ncz = ncz + aaa*nqz
2504 nni = (nqx*ncx + nqy*ncy + nqz*ncz)
2505 ENDIF
2506
2507 aaa = pene(i)/nni
2508 ncx = ncx*aaa
2509 ncy = ncy*aaa
2510 ncz = ncz*aaa
2511
2512 xpc = xc(i) - ncx
2513 ypc = yc(i) - ncy
2514 zpc = zc(i) - ncz
2515
2516 nnx=(ypb-ypa)*(zpc-zpa) - (zpb-zpa)*(ypc-ypa)
2517 nny=(zpb-zpa)*(xpc-xpa) - (xpb-xpa)*(zpc-zpa)
2518 nnz=(xpb-xpa)*(ypc-ypa) - (ypb-ypa)*(xpc-xpa)
2519 aaa = nnx*nqx + nny*nqy + nnz*nqz
2520 IF(aaa <= zero)THEN
2521 nax = nqx*pene(i)
2522 nay = nqy*pene(i)
2523 naz = nqz*pene(i)
2524 nbx = nax
2525 nby = nay
2526 nbz = naz
2527 ncx = nax
2528 ncy = nay
2529 ncz = naz
2530 xpa = xa(i) - nax
2531 ypa = ya(i) - nay
2532 zpa = za(i) - naz
2533 xpb = xb(i) - nbx
2534 ypb = yb(i) - nby
2535 zpb = zb(i) - nbz
2536 xpc = xc(i) - ncx
2537 ypc = yc(i) - ncy
2538 zpc = zc(i) - ncz
2539 ENDIF
2540
2541 xab = xpb-xpa
2542 yab = ypb-ypa
2543 zab = zpb-zpa
2544 xbc = xpc-xpb
2545 ybc = ypc-ypb
2546 zbc = zpc-zpb
2547 xca = xpa-xpc
2548 yca = ypa-ypc
2549 zca = zpa-zpc
2550
2551 xia = xpa-xi(i)
2552 yia = ypa-yi(i)
2553 zia = zpa-zi(i)
2554 xib = xpb-xi(i)
2555 yib = ypb-yi(i)
2556 zib = zpb-zi(i)
2557 xic = xpc-xi(i)
2558 yic = ypc-yi(i)
2559 zic = zpc-zi(i)
2560 sx = - yab*zca + zab*yca
2561 sy = - zab*xca + xab*zca
2562 sz = - xab*yca + yab*xca
2563 s2 = sx*sx+sy*sy+sz*sz
2564 sax = yib*zic - zib*yic
2565 say = zib*xic - xib*zic
2566 saz = xib*yic - yib*xic
2567 la(i) = (sx*sax+sy*say+sz*saz)/s2
2568 sbx = yic*zia - zic*yia
2569 sby = zic*xia - xic*zia
2570 sbz = xic*yia - yic*xia
2571 lb(i) = (sx*sbx+sy*sby+sz*sbz)/s2
2572 lc(i) = one - la(i) - lb(i)
2573
2574 IF(la(i)>=zero.and.lb(i)>=zero.and.lc(i)>=zero)THEN
2575
2576
2577
2578
2579
2580 IF (ikeep==1) THEN
2581 pene(i) = zero
2582 icontact(i) = 0
2583 cycle
2584 END IF
2585 n1(i) = nqx
2586 n2(i) = nqy
2587 n3(i) = nqz
2588 ELSE
2589
2590 bbb = pene_tm(i)-pene(i)
2591 IF (bbb>tol_b.AND.la(i)<zero) THEN
2592 pene(i) = pene_tm(i)
2593 itq = subtria_tm(i)
2594
2595 subtria(i) = itq
2596 ix1(i) = ixx(i,ic( 7,itq))
2597 ix2(i) = ixx(i,ic( 8,itq))
2598 ix3(i) = ixx(i,ic( 9,itq))
2599 ix4(i) = ixx(i,ic(10,itq))
2600 x1(i) = xx0(i,ic( 7,itq))
2601 x2(i) = xx0(i,ic( 8,itq))
2602 x3(i) = xx0(i,ic( 9,itq))
2603 x4(i) = xx0(i,ic(10,itq))
2604 y1(i) = yy0(i,ic( 7,itq))
2605 y2(i) = yy0(i,ic( 8,itq))
2606 y3(i) = yy0(i,ic( 9,itq))
2607 y4(i) = yy0(i,ic(10,itq))
2608 z1(i) = zz0(i,ic( 7,itq))
2609 z2(i) = zz0(i,ic( 8,itq))
2610 z3(i) = zz0(i,ic( 9,itq))
2611 z4(i) = zz0(i,ic(10,itq))
2612
2613 vx1(i) = vx(i,ic( 7,itq))
2614 vx2(i) = vx(i,ic( 8,itq))
2615 vx3(i) = vx(i,ic( 9,itq))
2616 vx4(i) = vx(i,ic(10,itq))
2617 vy1(i) = vy(i,ic( 7,itq))
2618 vy2(i) = vy(i,ic( 8,itq))
2619 vy3(i) = vy(i,ic( 9,itq))
2620 vy4(i) = vy(i,ic(10,itq))
2621 vz1(i) = vz(i,ic( 7,itq))
2622 vz2(i) = vz(i,ic( 8,itq))
2623 vz3(i) = vz(i,ic( 9,itq))
2624 vz4(i) = vz(i,ic(10,itq))
2625
2626 la(i) = la_tm(i)
2627 lb(i) = lb_tm(i)
2628 lc(i) = one - la(i) - lb(i)
2629 IF(la(i)<zero)THEN
2630 IF(lb(i)<zero)THEN
2631 la(i) = zero
2632 lb(i) = zero
2633 lc(i) = one
2634 ELSEIF(lc(i)<zero)THEN
2635 lc(i) = zero
2636 la(i) = zero
2637 lb(i) = one
2638 ELSE
2639 la(i) = zero
2640 aaa = lb(i) + lc(i)
2641 lb(i) = lb(i)/aaa
2642 lc(i) = lc(i)/aaa
2643 ENDIF
2644 ELSEIF(lb(i)<zero)THEN
2645 IF(lc(i)<zero)THEN
2646 lb(i) = zero
2647 lc(i) = zero
2648 la(i) = one
2649 ELSE
2650 lb(i) = zero
2651 aaa = lc(i) + la(i)
2652 lc(i) = lc(i)/aaa
2653 la(i) = la(i)/aaa
2654 ENDIF
2655 ELSEIF(lc(i)<zero)THEN
2656 lc(i) = zero
2657 aaa = la(i) + lb(i)
2658 la(i) = la(i)/aaa
2659 lb(i) = lb(i)/aaa
2660 ENDIF
2661
2662 IF (ix1(i) == ix2(i))THEN
2663 h1(i) = la(i)
2664 h2(i) = zero
2665 h3(i) = lb(i)
2666 h4(i) = lc(i)
2667 ELSEIF(ix2(i) == ix3(i))THEN
2668 h1(i) = la(i)
2669 h2(i) = lb(i)
2670 h3(i) = zero
2671 h4(i) = lc(i)
2672
2673 ELSEIF(ix4(i) == ix1(i))THEN
2674 h1(i) = lc(i)
2675 h2(i) = la(i)
2676 h3(i) = lb(i)
2677 h4(i) = zero
2678 ELSE
2679 h0 = fourth * la(i)
2680 h1(i) = h0
2681 h2(i) = h0
2682 h3(i) = h0 + lb(i)
2683 h4(i) = h0 + lc(i)
2684 ENDIF
2685 n1(i) = n_tm(1,i)
2686 n2(i) = n_tm(2,i)
2687 n3(i) = n_tm(3,i)
2688 cycle
2689 END IF
2690 n1(i) = nqx
2691 n2(i) = nqy
2692 n3(i) = nqz
2693 IF(la(i)<zero.and.lb(i)<zero)THEN
2694 la(i) = zero
2695 lb(i) = zero
2696 lc(i) = one
2697 ELSEIF(lb(i)<zero.and.lc(i)<zero)THEN
2698 lb(i) = zero
2699 lc(i) = zero
2700 la(i) = one
2701 ELSEIF(lc(i)<zero.and.la(i)<zero)THEN
2702 lc(i) = zero
2703 la(i) = zero
2704 lb(i) = one
2705 ELSEIF(la(i)<zero)THEN
2706 la(i) = zero
2707 aaa = lb(i) + lc(i)
2708 lb(i) = lb(i)/aaa
2709 lc(i) = lc(i)/aaa
2710 ELSEIF(lb(i)<zero)THEN
2711 lb(i) = zero
2712 aaa = lc(i) + la(i)
2713 lc(i) = lc(i)/aaa
2714 la(i) = la(i)/aaa
2715 ELSEIF(lc(i)<zero)THEN
2716 lc(i) = zero
2717 aaa = la(i) + lb(i)
2718 la(i) = la(i)/aaa
2719 lb(i) = lb(i)/aaa
2720 ENDIF
2721 goto 3456
2722 xpi = la(i)*xa(i) + lb(i)*xb(i) + lc(i)*xc(i)
2723 ypi = la(i)*ya(i) + lb(i)*yb(i) + lc(i)*yc(i)
2724 zpi = la(i)*za(i) + lb(i)*zb(i) + lc(i)*zc(i)
2725 n1(i) = xpi - xi(i)
2726 n2(i) = ypi - yi(i)
2727 n3(i) = zpi - zi(i)
2728 pene(i) = sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2)
2729 s2 = one/
max(em30,pene(i))
2730 n1(i) = n1(i)*s2
2731 n2(i) = n2(i)*s2
2732 n3(i) = n3(i)*s2
2733 3456 continue
2734 ENDIF
2735
2736
2737 IF ((impl_s>0 .AND. ittoff>0)) THEN
2738
2739
2740
2741 overw = overw0
2742 n1(i) = la(i)*nax + lb(i)*nbx + lc(i)*ncx
2743 n2(i) = la(i)*nay + lb(i)*nby + lc(i)*ncy
2744 n3(i) = la(i)*naz + lb(i)*nbz + lc(i)*ncz
2745
2746
2747
2748
2749
2750
2751
2752 s2 = one/
max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
2753 n1(i) = n1(i)*s2
2754 n2(i) = n2(i)*s2
2755 n3(i) = n3(i)*s2
2756
2757 aaa = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
2758 aaa =
max(one,one/
max(aaa,em20))
2759 pene(i) = pene(i)*aaa
2760
2761
2762
2763
2764
2765
2766 END IF
2767
2768 ix1(i) = ixx(i,ic( 7,itq))
2769 ix2(i) = ixx(i,ic( 8,itq))
2770 ix3(i) = ixx(i,ic( 9,itq))
2771 ix4(i) = ixx(i,ic(10,itq))
2772
2773 x1(i) = xx0(i,ic( 7,itq))
2774 x2(i) = xx0(i,ic( 8,itq))
2775 x3(i) = xx0(i,ic( 9,itq))
2776 x4(i) = xx0(i,ic(10,itq))
2777 y1(i) = yy0(i,ic( 7,itq))
2778 y2(i) = yy0(i,ic( 8,itq))
2779 y3(i) = yy0(i,ic( 9,itq))
2780 y4(i) = yy0(i,ic(10,itq))
2781 z1(i) = zz0(i,ic( 7,itq))
2782 z2(i) = zz0(i,ic( 8,itq))
2783 z3(i) = zz0(i,ic( 9,itq))
2784 z4(i) = zz0(i,ic(10,itq))
2785
2786 vx1(i) = vx(i,ic( 7,itq))
2787 vx2(i) = vx(i,ic( 8,itq))
2788 vx3(i) = vx(i,ic( 9,itq))
2789 vx4(i) = vx(i,ic(10,itq))
2790 vy1(i) = vy(i,ic( 7,itq))
2791 vy2(i) = vy(i,ic( 8,itq))
2792 vy3(i) = vy(i,ic( 9,itq))
2793 vy4(i) = vy(i,ic(10,itq))
2794 vz1(i) = vz(i,ic( 7,itq))
2795 vz2(i) = vz(i,ic( 8,itq))
2796 vz3(i) = vz(i,ic( 9,itq))
2797 vz4(i) = vz(i,ic(10,itq))
2798
2799 IF (ix1(i) == ix2(i))THEN
2800 h1(i) = la(i)
2801 h2(i) = zero
2802 h3(i) = lb(i)
2803 h4(i) = lc(i)
2804 ELSEIF(ix2(i) == ix3(i))THEN
2805 h1(i) = la(i)
2806 h2(i) = lb(i)
2807 h3(i) = zero
2808 h4(i) = lc(i)
2809
2810 ELSEIF(ix4(i) == ix1(i))THEN
2811 h1(i) = lc(i)
2812 h2(i) = la(i)
2813 h3(i) = lb(i)
2814 h4(i) = zero
2815 ELSE
2816 h0 = fourth * la(i)
2817 h1(i) = h0
2818 h2(i) = h0
2819 h3(i) = h0 + lb(i)
2820 h4(i) = h0 + lc(i)
2821 ENDIF
2822
2823
2824 n=cand_n(i)
2825
2826
2827
2828
2829
2830! write(iout,*)'ix3=',ix3(i),'ix4=',ix4(i)
2831
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841
2842
2843
2844
2845
2846
2847
2848
2849 mglob = mseglo(cand_e(i))
2850 IF (itq==5 .or. itq==9 .or. itq==13 .or. itq==17)THEN
2851 mglob = mvoisin(1,cand_e(i))
2852 itq = itriv(1,i)
2853 ELSEIF(itq==6 .or. itq==10 .or. itq==14 .or. itq==18)THEN
2854 mglob = mvoisin(2,cand_e(i))
2855 itq = itriv(2,i)
2856 ELSEIF(itq==7 .or. itq==11 .or. itq==15 .or. itq==19)THEN
2857 mglob = mvoisin(3,cand_e(i))
2858 itq = itriv(3,i)
2859 ELSEIF(itq==8 .or. itq==12 .or. itq==16 .or. itq==20)THEN
2860 mglob = mvoisin(4,cand_e(i))
2861 itq = itriv(4,i)
2862 ENDIF
2863 IF(n <= nsn)THEN
2864 IF(irtlm(1,n)>0) THEN
2865 irtlm(1,n) = mglob
2866 time_s(n) = -ep20
2867 irtlm(2,n)= itq
2868
2869
2870
2871
2872
2873
2874 ELSEIF((time_s(n) < pene(i)) .or.
2875 . (time_s(n) == pene(i) .and.
2876 . -irtlm(1,n) < mglob ) )THEN
2877 irtlm(2,n)= itq
2878 irtlm(1,n) = -mglob
2879 time_s(n) = pene(i)
2880 ENDIF
2881 ELSE
2882 ii=n-nsn
2887 ELSEIF((
time_sfi(nin)%P(ii) < pene(i)).or.
2888 . (
time_sfi(nin)%P(ii)== pene(i).and.
2889 . -
irtlm_fi(nin)%P(1,ii) < mglob) )
THEN
2893 ENDIF
2894 ENDIF
2895 ENDDO
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925! | 15 | 16 11 4 |(12) 3 10 | 3 10 11 4 | | \ / |
2926
2927
2928
2929
2930
2931
2932
2933
2934 DO i=1,jlt
2935 IF(istep(i) /= 6 .OR. subtria(i)>20)cycle
2936
2937 itq = subtria(i)
2938 IF (ixx(i,4)==ixx(i,3)) itq =1
2939 l = cand_e(i)
2940 k = ic(1,itq)
2941 xa(i) = xx0(i,k)
2942 ya(i) = yy0(i,k)
2943 za(i) = zz0(i,k)
2944
2945 k = ic(2,itq)
2946 xb(i) = xx0(i,k)
2947 yb(i) = yy0(i,k)
2948 zb(i) = zz0(i,k)
2949
2950 k = ic(3,itq)
2951 xc(i) = xx0(i,k)
2952 yc(i) = yy0(i,k)
2953 zc(i) = zz0(i,k)
2954
2955 k = ic(5,itq)
2956 xd(i) = xx0(i,k)
2957 yd(i) = yy0(i,k)
2958 zd(i) = zz0(i,k)
2959
2960 k = ic(6,itq)
2961 xe(i) = xx0(i,k)
2962 ye(i) = yy0(i,k)
2963 ze(i) = zz0(i,k)
2964 xib = xb(i)-xi(i)
2965 yib = yb(i)-yi(i)
2966 zib = zb(i)-zi(i)
2967 xbc = xc(i)-xb(i)
2968 ybc = yc(i)-yb(i)
2969 zbc = zc(i)-zb(i)
2970 aaa = xib*xbc+yib*ybc+zib*zbc
2971 xbd = xd(i)-xb(i)
2972 ybd = yd(i)-yb(i)
2973 zbd = zd(i)-zb(i)
2974 xce = xe(i)-xc(i)
2975 yce = ye(i)-yc(i)
2976 zce = ze(i)-zc(i)
2977 bbb = xib*xbd+yib*ybd+zib*zbd
2978 subtria_n(i) = 0
2979
2980 IF (aaa>zero) THEN
2981 IF(itq ==1.AND.mvoisin(4,l)>0)THEN
2982 subtria_n(i) = 16
2983 ELSEIF(itq ==2.AND.mvoisin(1,l)>0)THEN
2984 subtria_n(i) = 13
2985 ELSEIF(itq ==3.AND.mvoisin(2,l)>0)THEN
2986 subtria_n(i) = 14
2987 ELSEIF(itq ==4.AND.mvoisin(3,l)>0)THEN
2988 subtria_n(i) = 15
2989
2990 ELSE
2991 icontact(i) = 0
2992 istep(i) = 0
2993 cycle
2994 ENDIF
2995
2996 ELSEIF (bbb<zero) THEN
2997 nni =zero
2998 IF(itq ==1.AND.mvoisin(4,l)==0)THEN
2999 nqx = sx4150(i)
3000 nqy = sy4150(i)
3001 nqz = sz4150(i)
3002 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3003 . (nqx*ybd - nqy*xbd)*zib
3004 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 4
3005 ELSEIF(itq ==2.AND.mvoisin(1,l)==0)THEN
3006 nqx = sx1250(i)
3007 nqy = sy1250(i)
3008 nqz = sz1250(i)
3009 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3010 . (nqx*ybd - nqy*xbd)*zib
3011 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 1
3012 ELSEIF(itq ==3.AND.mvoisin(2,l)==0)THEN
3013 nqx = sx3450(i)
3014 nqy = sy3450(i)
3015 nqz = sz3450(i)
3016 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3017 . (nqx*ybd - nqy*xbd)*zib
3018 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 2
3019 ELSEIF(itq ==4.AND.mvoisin(3,l)==0)THEN
3020 nqx = sx2350(i)
3021 nqy = sy2350(i)
3022 nqz = sz2350(i)
3023 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3024 . (nqx*ybd - nqy*xbd)*zib
3025 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 3
3026 ENDIF
3027 IF (nni>zero) THEN
3028 icontact(i) = 0
3029 istep(i) = 0
3030 cycle
3031 END IF
3032 END IF
3033
3034 xic = xc(i)-xi(i)
3035 yic = yc(i)-yi(i)
3036 zic = zc(i)-zi(i)
3037 aaa = xic*xbc+yic*ybc+zic*zbc
3038 bbb = xic*xce+yic*yce+zic*zce
3039
3040 IF (aaa<zero.AND.subtria_n(i)==0) THEN
3041 IF(itq ==1.AND.mvoisin(2,l)>0)THEN
3042 subtria_n(i) = 10
3043 ELSEIF(itq ==2.AND.mvoisin(3,l)>0)THEN
3044 subtria_n(i) = 11
3045 ELSEIF(itq ==3.AND.mvoisin(4,l)>0)THEN
3046 subtria_n(i) = 12
3047 ELSEIF(itq ==4.AND.mvoisin(1,l)>0)THEN
3048 subtria_n(i) = 9
3049
3050 ELSE
3051 icontact(i) = 0
3052 istep(i) = 0
3053 cycle
3054 ENDIF
3055
3056 ELSEIF (bbb>zero) THEN
3057 nni =zero
3058 IF(itq ==1.AND.mvoisin(2,l)==0)THEN
3059 nqx = sx2350(i)
3060 nqy = sy2350(i)
3061 nqz = sz2350(i)
3062 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3063 . (nqx*yce - nqy*xce)*zic
3064 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 2
3065 ELSEIF(itq ==2.AND.mvoisin(1,l)==0)THEN
3066 nqx = sx3450(i)
3067 nqy = sy3450(i)
3068 nqz = sz3450(i)
3069 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3070 . (nqx*yce - nqy*xce)*zic
3071 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 3
3072 ELSEIF(itq ==3.AND.mvoisin(2,l)==0)THEN
3073 nqx = sx4150(i)
3074 nqy = sy4150(i)
3075 nqz = sz4150(i)
3076 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3077 . (nqx*yce - nqy*xce)*zic
3078 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 4
3079 ELSEIF(itq ==4.AND.mvoisin(3,l)==0)THEN
3080 nqx = sx1250(i)
3081 nqy = sy1250(i)
3082 nqz = sz1250(i)
3083 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3084 . (nqx*yce - nqy*xce)*zic
3085 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 1
3086 ENDIF
3087 IF (nni>zero) THEN
3088 icontact(i) = 0
3089 istep(i) = 0
3090 cycle
3091 END IF
3092 END IF
3093 IF (subtria_n(i)>0) subtria(i) = subtria_n(i)
3094 subtria_n(i) = itq
3095 END DO
3096
3097 DO i=1,jlt
3098 IF(istep(i) /= 6 )cycle
3099
3100 itq = subtria(i)
3101 IF (itq>20) itq = itq -20
3102
3103 k = ic(1,itq)
3104 xa(i) = xx0(i,k)
3105 ya(i) = yy0(i,k)
3106 za(i) = zz0(i,k)
3107
3108 k = ic(2,itq)
3109 xb(i) = xx0(i,k)
3110 yb(i) = yy0(i,k)
3111 zb(i) = zz0(i,k)
3112
3113 k = ic(3,itq)
3114 xc(i) = xx0(i,k)
3115 yc(i) = yy0(i,k)
3116 zc(i) = zz0(i,k)
3117
3118 xbc = xc(i)-xb(i)
3119 ybc = yc(i)-yb(i)
3120 zbc = zc(i)-zb(i)
3121
3122 IF(itq ==1)THEN
3123 nqx = sx1250(i)
3124 nqy = sy1250(i)
3125 nqz = sz1250(i)
3126 ELSEIF(itq ==2)THEN
3127 nqx = sx2350(i)
3128 nqy = sy2350(i)
3129 nqz = sz2350(i)
3130 ELSEIF(itq ==3)THEN
3131 nqx = sx3450(i)
3132 nqy = sy3450(i)
3133 nqz = sz3450(i)
3134 ELSEIF(itq ==4)THEN
3135 nqx = sx4150(i)
3136 nqy = sy4150(i)
3137 nqz = sz4150(i)
3138 ELSE
3139
3140 xab = xb(i)-xa(i)
3141 yab = yb(i)-ya(i)
3142 zab = zb(i)-za(i)
3143 nqx = yab*zbc - zab*ybc
3144 nqy = zab*xbc - xab*zbc
3145 nqz = xab*ybc - yab*xbc
3146 ENDIF
3147 nx(i) = nqy*zbc - nqz*ybc
3148 ny(i) = nqz*xbc - nqx*zbc
3149 nz(i) = nqx*ybc - nqy*xbc
3150
3151 nqx = -nx(i)
3152 nqy = -ny(i)
3153 nqz = -nz(i)
3154
3155 aaa =
max(em30,sqrt(nqx*nqx + nqy*nqy + nqz*nqz))
3156 s2 = one/aaa
3157 nqx = nqx*s2
3158 nqy = nqy*s2
3159 nqz = nqz*s2
3160
3161 bbb = (xi(i)-xb(i))*nqx+(yi(i)-yb(i))*nqy+(zi(i)-zb(i))*nqz
3162 pene(i) = -
min(zero,bbb)
3163 IF(pene(i) == zero)THEN
3164 icontact(i) = 0
3165 cycle
3166 ENDIF
3167 n1(i) = nqx
3168 n2(i) = nqy
3169 n3(i) = nqz
3170 ix1(i) = ixx(i,ic( 7,itq))
3171 ix2(i) = ixx(i,ic( 8,itq))
3172 ix3(i) = ixx(i,ic( 9,itq))
3173 ix4(i) = ixx(i,ic(10,itq))
3174
3175 x1(i) = xx0(i,ic( 7,itq))
3176 x2(i) = xx0(i,ic( 8,itq))
3177 x3(i) = xx0(i,ic( 9,itq))
3178 x4(i) = xx0(i,ic(10,itq))
3179 y1(i) = yy0(i,ic( 7,itq))
3180 y2(i) = yy0(i,ic( 8,itq))
3181 y3(i) = yy0(i,ic( 9,itq))
3182 y4(i) = yy0(i,ic(10,itq))
3183 z1(i) = zz0(i,ic( 7,itq))
3184 z2(i) = zz0(i,ic( 8,itq))
3185 z3(i) = zz0(i,ic( 9,itq))
3186 z4(i) = zz0(i,ic(10,itq))
3187
3188 vx1(i) = vx(i,ic( 7,itq))
3189 vx2(i) = vx(i,ic( 8,itq))
3190 vx3(i) = vx(i,ic( 9,itq))
3191 vx4(i) = vx(i,ic(10,itq))
3192 vy1(i) = vy(i,ic( 7,itq))
3193 vy2(i) = vy(i,ic( 8,itq))
3194 vy3(i) = vy(i,ic( 9,itq))
3195 vy4(i) = vy(i,ic(10,itq))
3196 vz1(i) = vz(i,ic( 7,itq))
3197 vz2(i) = vz(i,ic( 8,itq))
3198 vz3(i) = vz(i,ic( 9,itq))
3199 vz4(i) = vz(i,ic(10,itq))
3200 h1(i) = zero
3201 h2(i) = zero
3202 h3(i) = half
3203 h4(i) = half
3204
3205 n=cand_n(i)
3206 mglob = mseglo(cand_e(i))
3207 IF (itq==5 .or. itq==9 .or. itq==13 .or. itq==17)THEN
3208 mglob = mvoisin(1,cand_e(i))
3209
3210 ELSEIF(itq==6 .or. itq==10 .or. itq==14 .or. itq==18)THEN
3211 mglob = mvoisin(2,cand_e(i))
3212
3213 ELSEIF(itq==7 .or. itq==11 .or. itq==15 .or. itq==19)THEN
3214 mglob = mvoisin(3,cand_e(i))
3215
3216 ELSEIF(itq==8 .or. itq==12 .or. itq==16 .or. itq==20)THEN
3217 mglob = mvoisin(4,cand_e(i))
3218
3219 ENDIF
3220 IF (itq>8) itq = subtria_n(i)
3221
3222 IF(n <= nsn)THEN
3223 IF(irtlm(1,n)>0) THEN
3224 irtlm(1,n) = mglob
3225 time_s(n) = -ep20
3226 irtlm(2,n)= itq+20
3227 ELSEIF((time_s(n) < pene(i)) .or.
3228 . (time_s(n) == pene(i) .and.
3229 . -irtlm(1,n) < mglob ) )THEN
3230 irtlm(2,n)= itq+20
3231 irtlm(1,n) = -mglob
3232 time_s(n) = pene(i)
3233 ENDIF
3234 ELSE
3235 ii=n-nsn
3240 ELSEIF((
time_sfi(nin)%P(ii) < pene(i)).or.
3241 . (
time_sfi(nin)%P(ii)== pene(i).and.
3242 . -
irtlm_fi(nin)%P(1,ii) < mglob) )
THEN
3246 ENDIF
3247 ENDIF
3248 END DO
3249
3250
3251
3252
3253
3254
3255 DO i=1,jlt
3256 ce_loc(i) = cand_e(i)
3257 cn_loc(i) = cand_n(i)
3258 IF(pene(i) == 0 )THEN
3259
3260 n=cand_n(i)
3261 IF(n <= nsn)THEN
3262
3263
3264 IF(irtlm(1,n) > 0)THEN
3265 irtlm(1,n)=0
3266 time_s(n) = -ep20
3267 secnd_fr(4,n)= zero
3268 secnd_fr(5,n)= zero
3269 secnd_fr(6,n)= zero
3270 IF (imconv==1) pene_old(2,n) = zero
3271 pene_old(5,n) = zero
3272 IF (imconv==1) stif_old(2,n) = zero
3273 ENDIF
3274 ELSE
3275 IF(
irtlm_fi(nin)%P(1,n-nsn) > 0)
THEN
3281 IF (imconv==1)
pene_oldfi(nin)%P(2,n-nsn) = zero
3283 IF (imconv==1)
stif_oldfi(nin)%P(2,n-nsn) = zero
3284 ENDIF
3285 ENDIF
3286 ELSE
3287 ns=nsvg(i)
3288 IF (nm1(i)==zero.AND.nm2(i)==zero.AND.nm3(i)==zero) THEN
3289 nm1(i) = n1(i)
3290 nm2(i) = n2(i)
3291 nm3(i) = n3(i)
3292 END IF
3293
3294 IF (ix1(i) == ix2(i))THEN
3295 ELSEIF(ix2(i) == ix3(i))THEN
3296
3297 ELSEIF(ix4(i) == ix1(i))THEN
3298 ELSE
3299
3300 ENDIF
3301 xs = h1(i)*x1(i) + h2(i)*x2(i) + h3(i)*x3(i) + h4(i)*x4(i)
3302 ys = h1(i)*y1(i) + h2(i)*y2(i) + h3(i)*y3(i) + h4(i)*y4(i)
3303 zs = h1(i)*z1(i) + h2(i)*z2(i) + h3(i)*z3(i) + h4(i)*z4(i)
3304 xs = xs - xi(i)
3305 ys = ys - yi(i)
3306 zs = zs - zi(i)
3307
3308 aaa = xs*n1(i) + ys*n2(i) + zs*n3(i)
3309
3310 IF(aaa > zero)THEN
3311
3312 aaa = xs**2 + ys**2 + zs**2
3313 aaa = onep01*sqrt(aaa)
3314 pmax_gap =
max(pmax_gap,aaa)
3315
3316 n=cand_n(i)
3317 IF(n <= nsn)THEN
3318 pene_old(3,n) =
max(pene_old(3,n),aaa)
3319 ELSE
3322 ENDIF
3323 ENDIF
3324 ENDIF
3325 ENDDO
3326#include "lockoff.inc"
3327
3328
3329 IF (inacti==5) THEN
3330 IF (tt==zero) THEN
3331 DO i=1,jlt
3332 IF(pene(i) == zero)cycle
3333 jg = nsvg(i)
3334 n = cand_n(i)
3335 IF(jg > 0)THEN
3336#include "lockon.inc"
3337 pene_old(5,n) =
max( pene(i) ,pene_old(5,n) )
3338 stif_old(1,n) =
max( stif(i) ,stif_old(1,n) )
3339#include "lockoff.inc"
3340 ELSE
3341 jg = -jg
3342#include "lockon.inc"
3345#include "lockoff.inc"
3346 ENDIF
3347 ENDDO
3348
3349 ELSE
3350 DO i=1,jlt
3351 IF(pene(i) == zero)cycle
3352 jg = nsvg(i)
3353 n = cand_n(i)
3354 IF(jg > 0)THEN
3355 IF (ipen0(i)==0)THEN
3356#include "lockon.inc"
3357 pene_old(5,n) =
max( pene(i) ,pene_old(5,n) )
3358 stif_old(1,n) =
max( stif(i) ,stif_old(1,n) )
3359#include "lockoff.inc"
3360 END IF
3361 ELSE
3362 jg = -jg
3363 IF (ipen0(i)==0)THEN
3364#include "lockon.inc"
3367#include "lockoff.inc"
3368 END IF
3369 ENDIF
3370 ENDDO
3371 END IF
3372
3373 DO i=1,jlt
3374 IF(pene(i) == zero)cycle
3375 jg = nsvg(i)
3376 n = cand_n(i)
3377 IF(jg > 0)THEN
3378 pene_sh =
max(zero,pene(i)-pene_old(5,n))
3379 ELSE
3380 jg = -jg
3382 ENDIF
3383
3384 pene(i) = pene_sh
3385 ENDDO
3386 END IF
3387
3388 penref(1:jlt) = zero
3389 IF (inacti==-1) THEN
3390 DO i=1,jlt
3391 IF(pene(i) == zero.OR.iknon(i)==0) cycle
3392 penref(i) = sqrt(eps2*
area(i))
3393 jg = nsvg(i)
3394 n = cand_n(i)
3395 IF(jg > 0)THEN
3396 pene_sh = em01*pene_old(5,n)
3397 ELSE
3398 jg = -jg
3400 ENDIF
3401 penref(i) =
max(penref(i),pene_sh)
3402 IF (iknon(i)<0) THEN
3403 f_pene = em01*pene(i)/pene_sh
3404 penref(i) = zero
3405 IF (f_pene>two_third) THEN
3406 fac_k =
min(twenty,six*f_pene)
3407 penref(i) = pene_sh/sqrt(fac_k)
3408 END IF
3409 END IF
3410 ENDDO
3411 ELSE
3412 DO i=1,jlt
3413 IF(pene(i) == zero.OR.iknon(i)==0) cycle
3414 n=cand_n(i)
3415 l=cand_e(i)
3416 penref(i) = one_fifth*sqrt(
area(i))
3417 gap2 = one_fifth*(gaps(i)+gap_m(l))
3418 IF (gap2 > em04) penref(i)=
min(penref(i),gap2)
3419 jg = nsvg(i)
3420 IF(jg > 0)THEN
3421 pene_sh = pene_old(5,n)
3422 ELSE
3423 jg = -jg
3425 ENDIF
3426 penref(i) =
max(penref(i),pene_sh)
3427 IF(iknon(i)==1) penref(i) = ten*penref(i)
3428 IF(iknon(i)>2) penref(i) = one_fifth*penref(i)
3429 ENDDO
3430 END IF
3431
3432 RETURN
if(complex_arithmetic) id
subroutine area(d1, x, x2, y, y2, eint, stif0)
subroutine intersecv0(istep, subtria, ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, adt0, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, nx, ny, nz, sx125, sx235, sx345, sx415, sy125, sy235, sy345, sy415, sz125, sz235, sz345, sz415, xi, yi, zi, zeromt, unpt)
subroutine s_in_m4(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xi, yi, zi, ier)
subroutine intersecv(istep, subtria, ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, adt0, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, nx, ny, nz, sx125, sx235, sx345, sx415, sy125, sy235, sy345, sy415, sz125, sz235, sz345, sz415, xi, yi, zi, zerom, unp, zeromt, unpt)
subroutine intersecsh(istep, subtria0, ixx3, ixx4, intersect, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, adt0, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, nx, ny, nz, sx125, sx235, sx345, sx415, sy125, sy235, sy345, sy415, sz125, sz235, sz345, sz415, xi, yi, zi, sox125, sox235, sox345, sox415, soy125, soy235, soy345, soy415, soz125, soz235, soz345, soz415, mvoisin)
subroutine i24nexttria(izlim, istep, subtria_n, subtria, largep, pene, xxl, yyl, zzl, sx125, sy125, sz125, sx235, sy235, sz235, sx345, sy345, sz345, sx415, sy415, sz415, xi, yi, zi, n1, n2, n3, la, lb, lc, gap2, eps0)
type(real_pointer2), dimension(:), allocatable stif_oldfi
type(real_pointer2), dimension(:), allocatable secnd_frfi
type(real_pointer), dimension(:), allocatable time_sfi
type(int_pointer2), dimension(:), allocatable irtlm_fi
type(real_pointer2), dimension(:), allocatable pene_oldfi
type(int_pointer), dimension(:), allocatable icont_i_fi
static2 void intersect(char *uplo, char *diag, Int j, Int start, Int end, Int action, Int *ptrsizebuff, complex **pptrbuff, complex *ptrblock, Int m, Int n, MDESC *ma, Int ia, Int ja, Int templateheight0, Int templatewidth0, MDESC *mb, Int ib, Int jb, Int templateheight1, Int templatewidth1)
1.15.0