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
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
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,al1den,al2den,al3den,al4den,
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(),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! | 2 | 5 2 3 | 15 4 1 | 4 1 2 3 | | \ / |
219
220
221
222
223
224
225
226
227
228
229
230
231! +----+----------+----------+-------------+ |/ 16 \|/ \|/ 10 \|
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(i)=1
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(mvoisin(3,cand_e(i))/=0)THEN
625 sx43160(i) = y164*z163 - z164*y163
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.
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
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(i,1:17)
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
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
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
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
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
1424 v34 = sx345(i)*xi5 + sy345(i)*yi5 + sz345
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 dgap = dgap_nm(2,cand_e(i))
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! 8o------1------2------o3 | /
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! | 3 | 5 3 4 | 16 1 2 | 1 2 3 4 | |\ 19 /|
1723
1724
1725! | 5 | 14 2 1 | 5 6 7 | 6 7 2 1 | | 16 |
1726
1727
1728
1729
1730
1731
1732
1733
1734! +----+----------+----------+-------------+ |20/ \ 8| 4/ \ | / \ |
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(iTHEN
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
2046
2047! | 4 | 8 1 3 | 8 1 3 | - - - | - - - |
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)==ixx(i,11))THEN
2096 it(2,7,i)=-1
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(i,5)
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
2333 pty = nqtz*(xb(i)-xa(i)) - nqtx*(zb(i)-za(i))
2334 ptz = nqtx*(yb(i)-ya(i)) - nqty*(xb(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
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 2345 continue
2581 IF (ikeep==1) THEN
2582 pene(i) = zero
2583 icontact(i) = 0
2584 cycle
2585 END IF
2586 n1(i) = nqx
2587 n2(i) = nqy
2588 n3(i) = nqz
2589 ELSE
2590
2591 bbb = pene_tm(i)-pene(i)
2592 IF (bbb>tol_b.AND.la(i)<zero) THEN
2593 pene(i) = pene_tm(i)
2594 itq = subtria_tm(i)
2595
2596 subtria(i) = itq
2597 ix1(i) = ixx(i,ic( 7,itq))
2598 ix2(i) = ixx(i,ic( 8,itq))
2599 ix3(i) = ixx(i,ic( 9,itq))
2600 ix4(i) = ixx(i,ic(10,itq))
2601 x1(i) = xx0(i,ic( 7,itq))
2602 x2(i) = xx0(i,ic( 8,itq))
2603 x3(i) = xx0(i,ic( 9,itq))
2604 x4(i) = xx0(i,ic(10,itq))
2605 y1(i) = yy0(i,ic( 7,itq))
2606 y2(i) = yy0(i,ic( 8,itq))
2607 y3(i) = yy0(i,ic( 9,itq))
2608 y4(i) = yy0(i,ic(10,itq))
2609 z1(i) = zz0(i,ic( 7,itq))
2610 z2(i) = zz0(i,ic( 8,itq))
2611 z3(i) = zz0(i,ic( 9,itq))
2612 z4(i) = zz0(i,ic(10,itq))
2613
2614 vx1(i) = vx(i,ic( 7,itq))
2615 vx2(i) = vx(i,ic( 8,itq))
2616 vx3(i) = vx(i,ic( 9,itq))
2617 vx4(i) = vx(i,ic(10,itq))
2618 vy1(i) = vy(i,ic( 7,itq))
2619 vy2(i) = vy(i,ic( 8,itq))
2620 vy3(i) = vy(i,ic( 9,itq))
2621 vy4(i) = vy(i,ic(10,itq))
2622 vz1(i) = vz(i,ic( 7,itq))
2623 vz2(i) = vz(i,ic( 8,itq))
2624 vz3(i) = vz(i,ic( 9,itq))
2625 vz4(i) = vz(i,ic(10,itq))
2626
2627 la(i) = la_tm(i)
2628 lb(i) = lb_tm(i)
2629 lc(i) = one - la(i) - lb(i)
2630 IF(la(i)<zero)THEN
2631 IF(lb(i)<zero)THEN
2632 la(i) = zero
2633 lb(i) = zero
2634 lc(i) = one
2635 ELSEIF(lc(i)<zero)THEN
2636 lc(i) = zero
2637 la(i) = zero
2638 lb(i) = one
2639 ELSE
2640 la(i) = zero
2641 aaa = lb(i) + lc(i)
2642 lb(i) = lb(i)/aaa
2643 lc(i) = lc(i)/aaa
2644 ENDIF
2645 ELSEIF(lb(i)<zero)THEN
2646 IF(lc(i)<zero)THEN
2647 lb(i) = zero
2648 lc(i) = zero
2649 la(i) = one
2650 ELSE
2651 lb(i) = zero
2652 aaa = lc(i) + la(i)
2653 lc(i) = lc(i)/aaa
2654 la(i) = la(i)/aaa
2655 ENDIF
2656 ELSEIF(lc(i)<zero)THEN
2657 lc(i) = zero
2658 aaa = la(i) + lb(i)
2659 la(i) = la(i)/aaa
2660 lb(i) = lb(i)/aaa
2661 ENDIF
2662
2663 IF (ix1(i) == ix2(i))THEN
2664 h1(i) = la(i)
2665 h2(i) = zero
2666 h3(i) = lb(i)
2667 h4(i) = lc(i)
2668 ELSEIF(ix2(i) == ix3(i))THEN
2669 h1(i) = la(i)
2670 h2(i) = lb(i)
2671 h3(i) = zero
2672 h4(i) = lc(i)
2673
2674 ELSEIF(ix4(i) == ix1(i))THEN
2675 h1(i) = lc(i)
2676 h2(i) = la(i)
2677 h3(i) = lb(i)
2678 h4(i) = zero
2679 ELSE
2680 h0 = fourth * la(i)
2681 h1(i) = h0
2682 h2(i) = h0
2683 h3(i) = h0 + lb(i)
2684 h4(i) = h0 + lc(i)
2685 ENDIF
2686 n1(i) = n_tm(1,i)
2687 n2(i) = n_tm(2,i)
2688 n3(i) = n_tm(3,i)
2689 cycle
2690 END IF
2691 n1(i) = nqx
2692 n2(i) = nqy
2693 n3(i) = nqz
2694 IF(la(i)<zero.and.lb(i)<zero)THEN
2695 la(i) = zero
2696 lb(i) = zero
2697 lc(i) = one
2698 ELSEIF(lb(i)<zero.and.lc(i)<zero)THEN
2699 lb(i) = zero
2700 lc(i) = zero
2701 la(i) = one
2702 ELSEIF(lc(i)<zero.and.la(i)<zero)THEN
2703 lc(i) = zero
2704 la(i) = zero
2705 lb(i) = one
2706 ELSEIF(la(i)<zero)THEN
2707 la(i) = zero
2708 aaa = lb(i) + lc(i)
2709 lb(i) = lb(i)/aaa
2710 lc(i) = lc(i)/aaa
2711 ELSEIF(lb(i)<zero)THEN
2712 lb(i) = zero
2713 aaa = lc(i) + la(i)
2714 lc(i) = lc(i)/aaa
2715 la(i) = la(i)/aaa
2716 ELSEIF(lc(i)<zero)THEN
2717 lc(i) = zero
2718 aaa = la(i) + lb(i)
2719 la(i) = la(i)/aaa
2720 lb(i) = lb(i)/aaa
2721 ENDIF
2722 goto 3456
2723 xpi = la(i)*xa(i) + lb(i)*xb(i) + lc(i)*xc(i)
2724 ypi = la(i)*ya(i) + lb(i)*yb(i) + lc(i)*yc(i)
2725 zpi = la(i)*za(i) + lb(i)*zb(i) + lc(i)*zc(i)
2726 n1(i) = xpi - xi(i)
2727 n2(i) = ypi - yi(i)
2728 n3(i) = zpi - zi(i)
2729 pene(i) = sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2)
2730 s2 = one/
max(em30,pene(i))
2731 n1(i) = n1(i)*s2
2732 n2(i) = n2(i)*s2
2733 n3(i) = n3(i)*s2
2734 3456 continue
2735 ENDIF
2736
2737
2738 IF ((impl_s>0 .AND. ittoff>0)) THEN
2739
2740
2741
2742 overw = overw0
2743 n1(i) = la(i)*nax + lb(i)*nbx + lc(i)*ncx
2744 n2(i) = la(i)*nay + lb(i)*nby + lc(i)*ncy
2745 n3(i) = la(i)*naz + lb(i)*nbz + lc(i)*ncz
2746
2747
2748
2749
2750
2751
2752
2753 s2 = one/
max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
2754 n1(i) = n1(i)*s2
2755 n2(i) = n2(i)*s2
2756 n3(i) = n3(i)*s2
2757
2758 aaa = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
2759 aaa =
max(one,one/
max(aaa,em20))
2760 pene(i) = pene(i)*aaa
2761
2762
2763
2764
2765
2766
2767 END IF
2768
2769 ix1(i) = ixx(i,ic( 7,itq))
2770 ix2(i) = ixx(i,ic( 8,itq))
2771 ix3(i) = ixx(i,ic( 9,itq))
2772 ix4(i) = ixx(i,ic(10,itq))
2773
2774 x1(i) = xx0(i,ic( 7,itq))
2775 x2(i) = xx0(i,ic( 8,itq))
2776 x3(i) = xx0(i,ic( 9,itq))
2777 x4(i) = xx0(i,ic(10,itq))
2778 y1(i) = yy0(i,ic( 7,itq))
2779 y2(i) = yy0(i,ic( 8,itq))
2780 y3(i) = yy0(i,ic( 9,itq))
2781 y4(i) = yy0(i,ic(10,itq))
2782 z1(i) = zz0(i,ic( 7,itq))
2783 z2(i) = zz0(i,ic( 8,itq))
2784 z3(i) = zz0(i,ic( 9,itq))
2785 z4(i) = zz0(i,ic(10,itq))
2786
2787 vx1(i) = vx(i,ic( 7,itq))
2788 vx2(i) = vx(i,ic( 8,itq))
2789 vx3(i) = vx(i,ic( 9,itq))
2790 vx4(i) = vx(i,ic(10,itq))
2791 vy1(i) = vy(i,ic( 7,itq))
2792 vy2(i) = vy(i,ic( 8,itq))
2793 vy3(i) = vy(i,ic( 9,itq))
2794 vy4(i) = vy(i,ic(10,itq))
2795 vz1(i) = vz(i,ic( 7,itq))
2796 vz2(i) = vz(i,ic( 8,itq))
2797 vz3(i) = vz(i,ic( 9,itq))
2798 vz4(i) = vz(i,ic(10,itq))
2799
2800 IF (ix1(i) == ix2(i))THEN
2801 h1(i) = la(i)
2802 h2(i) = zero
2803 h3(i) = lb(i)
2804 h4(i) = lc(i)
2805 ELSEIF(ix2(i) == ix3(i))THEN
2806 h1(i) = la(i)
2807 h2(i) = lb(i)
2808 h3(i) = zero
2809 h4(i) = lc(i)
2810
2811 ELSEIF(ix4(i) == ix1(i))THEN
2812 h1(i) = lc(i)
2813 h2(i) = la(i)
2814 h3(i) = lb(i)
2815 h4(i) = zero
2816 ELSE
2817 h0 = fourth * la(i)
2818 h1(i) = h0
2819 h2(i) = h0
2820 h3(i) = h0 + lb(i)
2821 h4(i) = h0 + lc(i)
2822 ENDIF
2823
2824
2825 n=cand_n(i)
2826
2827
2828
2829
2830
2831
2832
2833
2834
2835
2836
2837
2838
2839
2840
2841
2842
2843
2844
2845
2846
2847
2848
2849
2850 mglob = mseglo(cand_e(i))
2851 IF (itq==5 .or. itq==9 .or. itq==13 .or. itq==17)THEN
2852 mglob = mvoisin(1,cand_e(i))
2853 itq = itriv(1,i)
2854 ELSEIF(itq==6 .or. itq==10 .or. itq==14 .or. itq==18)THEN
2855 mglob = mvoisin(2,cand_e(i))
2856 itq = itriv(2,i)
2857 ELSEIF(itq==7 .or. itq==11 .or. itq==15 .or. itq==19)THEN
2858 mglob = mvoisin(3,cand_e(i))
2859 itq = itriv(3,i)
2860 ELSEIF(itq==8 .or. itq==12 .or. itq==16 .or. itq==20)THEN
2861 mglob = mvoisin(4,cand_e(i))
2862 itq = itriv(4,i)
2863 ENDIF
2864 IF(n <= nsn)THEN
2865 IF(irtlm(1,n)>0) THEN
2866 irtlm(1,n) = mglob
2867 time_s(n) = -ep20
2868 irtlm(2,n)= itq
2869
2870
2871
2872
2873
2874
2875 ELSEIF((time_s(n) < pene(i)) .or.
2876 . (time_s(n) == pene(i) .and.
2877 . -irtlm(1,n) < mglob ) )THEN
2878 irtlm(2,n)= itq
2879 irtlm(1,n) = -mglob
2880 time_s(n) = pene(i)
2881 ENDIF
2882 ELSE
2883 ii=n-nsn
2888 ELSEIF((
time_sfi(nin)%P(ii) < pene(i)).or.
2889 . (
time_sfi(nin)%P(ii)== pene(i).and.
2890 . -
irtlm_fi(nin)%P(1,ii) < mglob) )
THEN
2894 ENDIF
2895 ENDIF
2896 ENDDO
2897
2898
2899
2900 123 continue
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
2926
2927
2928
2929
2930
2931
2932
2933
2934
2935
2936 DO i=1,jlt
2937 IF(istep(i) /= 6 .OR. subtria(i)>20)cycle
2938
2939 itq = subtria(i)
2940 IF (ixx(i,4)==ixx(i,3)) itq =1
2941 l = cand_e(i)
2942 k = ic(1,itq)
2943 xa(i) = xx0(i,k)
2944 ya(i) = yy0(i,k)
2945 za(i) = zz0(i,k)
2946
2947 k = ic(2,itq)
2948 xb(i) = xx0(i,k)
2949 yb(i) = yy0(i,k)
2950 zb(i) = zz0(i,k)
2951
2952 k = ic(3,itq)
2953 xc(i) = xx0(i,k)
2954 yc(i) = yy0(i,k)
2955 zc(i) = zz0(i,k)
2956
2957 k = ic(5,itq)
2958 xd(i) = xx0(i,k)
2959 yd(i) = yy0(i
2960 zd(i) = zz0(i,k)
2961
2962 k = ic(6,itq)
2963 xe(i) = xx0(i,k)
2964 ye(i) = yy0(i,k)
2965 ze(i) = zz0(i,k)
2966 xib = xb(i)-xi(i)
2967 yib = yb(i)-yi(i)
2968 zib = zb(i)-zi(i)
2969 xbc = xc(i)-xb(i)
2970 ybc = yc(i)-yb(i)
2971 zbc = zc(i)-zb(i)
2972 aaa = xib*xbc+yib*ybc+zib*zbc
2973 xbd = xd(i)-xb(i)
2974 ybd = yd(i)-yb(i)
2975 zbd = zd(i)-zb(i)
2976 xce = xe(i)-xc(i)
2977 yce = ye(i)-yc(i)
2978 zce = ze(i)-zc(i)
2979 bbb = xib*xbd+yib*ybd+zib*zbd
2980 subtria_n(i) = 0
2981
2982 IF (aaa>zero) THEN
2983 IF(itq ==1.AND.mvoisin(4,l)>0)THEN
2984 subtria_n(i) = 16
2985 ELSEIF(itq ==2.AND.mvoisin(1,l)>0)THEN
2986 subtria_n(i) = 13
2987 ELSEIF(itq ==3.AND.mvoisin(2,l)>0)THEN
2988 subtria_n(i) = 14
2989 ELSEIF(itq ==4.AND.mvoisin(3,l)>0)THEN
2990 subtria_n(i) = 15
2991
2992 ELSE
2993 icontact(i) = 0
2994 istep(i) = 0
2995 cycle
2996 ENDIF
2997
2998 ELSEIF (bbb<zero) THEN
2999 nni =zero
3000 IF(itq ==1.AND.mvoisin(4,l)==0)THEN
3001 nqx = sx4150(i)
3002 nqy = sy4150(i)
3003 nqz = sz4150(i)
3004 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3005 . (nqx*ybd - nqy*xbd)*zib
3006 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 4
3007 ELSEIF(itq ==2.AND.mvoisin(1,l)==0)THEN
3008 nqx = sx1250(i)
3009 nqy = sy1250(i)
3010 nqz = sz1250(i)
3011 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3012 . (nqx*ybd - nqy*xbd)*zib
3013 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 1
3014 ELSEIF(itq ==3.AND.mvoisin(2,l)==0)THEN
3015 nqx = sx3450(i)
3016 nqy = sy3450(i)
3017 nqz = sz3450(i)
3018 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3019 . (nqx*ybd - nqy*xbd)*zib
3020 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 2
3021 ELSEIF(itq ==4.AND.mvoisin(3,l)==0)THEN
3022 nqx = sx2350(i)
3023 nqy = sy2350(i)
3024 nqz = sz2350(i)
3025 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib
3026 . (nqx*ybd - nqy*xbd)*zib
3027 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 3
3028 ENDIF
3029 IF (nni>zero) THEN
3030 icontact(i) = 0
3031 istep(i) = 0
3032 cycle
3033 END IF
3034 END IF
3035
3036 xic = xc(i)-xi(i)
3037 yic = yc(i)-yi(i)
3038 zic = zc(i)-zi(i)
3039 aaa = xic*xbc+yic*ybc+zic*zbc
3040 bbb = xic*xce+yic*yce+zic*zce
3041
3042 IF (aaa<zero.AND.subtria_n(i)==0) THEN
3043 IF(itq ==1.AND.mvoisin(2,l)>0)THEN
3044 subtria_n(i) = 10
3045 ELSEIF(itq ==2.AND.mvoisin(3,l)>0)THEN
3046 subtria_n(i) = 11
3047 ELSEIF(itq ==3.AND.mvoisin(4,l)>0)THEN
3048 subtria_n(i) = 12
3049 ELSEIF(itq ==4.AND.mvoisin(1,l)>0)THEN
3050 subtria_n(i) = 9
3051
3052 ELSE
3053 icontact(i) = 0
3054 istep(i) = 0
3055 cycle
3056 ENDIF
3057
3058 ELSEIF (bbb>zero) THEN
3059 nni =zero
3060 IF(itq ==1.AND.mvoisin(2,l)==0)THEN
3061 nqx = sx2350(i)
3062 nqy = sy2350(i)
3063 nqz = sz2350(i)
3064 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3065 . (nqx*yce - nqy*xce)*zic
3066 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 2
3067 ELSEIF(itq ==2.AND.mvoisin(1,l)==0)THEN
3068 nqx = sx3450(i)
3069 nqy = sy3450(i)
3070 nqz = sz3450(i)
3071 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3072 . (nqx*yce - nqy*xce)*zic
3073 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 3
3074 ELSEIF(itq ==3.AND.mvoisin(2,l)==0)THEN
3075 nqx = sx4150(i)
3076 nqy = sy4150(i)
3077 nqz = sz4150(i)
3078 nni = (nqy*zce - nqz*yce)*xic+(nqz
3079 . (nqx*yce - nqy*xce)*zic
3080 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 4
3081 ELSEIF(itq ==4.AND.mvoisin(3,l)==0)THEN
3082 nqx = sx1250(i)
3083 nqy = sy1250(i)
3084 nqz = sz1250(i)
3085 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3086 . (nqx*yce - nqy*xce)*zic
3087 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 1
3088 ENDIF
3089 IF (nni>zero) THEN
3090 icontact(i) = 0
3091 istep(i) = 0
3092 cycle
3093 END IF
3094 END IF
3095 IF (subtria_n(i)>0) subtria(i) = subtria_n(i)
3096 subtria_n(i) = itq
3097 END DO
3098
3099 DO i=1,jlt
3100 IF(istep(i) /= 6 )cycle
3101
3102 itq = subtria(i)
3103 IF (itq>20) itq = itq -20
3104
3105 k = ic(1,itq)
3106 xa(i) = xx0(i,k)
3107 ya(i) = yy0(i,k)
3108 za(i) = zz0(i,k)
3109
3110 k = ic(2,itq)
3111 xb(i) = xx0(i,k)
3112 yb(i) = yy0(i,k)
3113 zb(i) = zz0(i,k)
3114
3115 k = ic(3,itq)
3116 xc(i) = xx0(i,k)
3117 yc(i) = yy0(i,k)
3118 zc(i) = zz0(i,k)
3119
3120 xbc = xc(i)-xb(i)
3121 ybc = yc(i)-yb(i)
3122 zbc = zc(i)-zb(i)
3123
3124 IF(itq ==1)THEN
3125 nqx = sx1250(i)
3126 nqy = sy1250(i)
3127 nqz = sz1250(i)
3128 ELSEIF(itq ==2)THEN
3129 nqx = sx2350(i)
3130 nqy = sy2350(i)
3131 nqz = sz2350(i)
3132 ELSEIF(itq ==3)THEN
3133 nqx = sx3450(i)
3134 nqy = sy3450(i)
3135 nqz = sz3450(i)
3136 ELSEIF(itq ==4)THEN
3137 nqx = sx4150(i)
3138 nqy = sy4150(i)
3139 nqz = sz4150(i)
3140 ELSE
3141
3142 xab = xb(i)-xa(i)
3143 yab = yb(i)-ya(i)
3144 zab = zb(i)-za(i)
3145 nqx = yab*zbc - zab*ybc
3146 nqy = zab*xbc - xab*zbc
3147 nqz = xab*ybc - yab*xbc
3148 ENDIF
3149 nx(i) = nqy*zbc - nqz*ybc
3150 ny(i) = nqz*xbc - nqx*zbc
3151 nz(i) = nqx*ybc - nqy
3152
3153 nqx = -nx(i)
3154 nqy = -ny(i)
3155 nqz = -nz(i)
3156
3158 s2 = one/aaa
3159 nqx = nqx*s2
3160 nqy = nqy*s2
3161 nqz = nqz*s2
3162
3163 bbb = (xi(i)-xb(i))*nqx+(yi(i)-yb(i))*nqy+(zi(i)-zb(i))*nqz
3164 pene(i) = -
min(zero,bbb)
3165 IF(pene(i) == zero)THEN
3166 icontact(i) = 0
3167 cycle
3168 ENDIF
3169 n1(i) = nqx
3170 n2(i) = nqy
3171 n3(i) = nqz
3172 ix1(i) = ixx(i,ic( 7,itq))
3173 ix2(i) = ixx(i,ic( 8,itq))
3174 ix3(i) = ixx(i,ic( 9,itq))
3175 ix4(i) = ixx(i,ic(10,itq))
3176
3177 x1(i) = xx0(i,ic( 7,itq))
3178 x2(i) = xx0(i,ic( 8,itq))
3179 x3(i) = xx0(i,ic( 9,itq))
3180 x4(i) = xx0(i,ic(10,itq))
3181 y1(i) = yy0(i,ic( 7,itq))
3182 y2(i) = yy0(i,ic( 8,itq))
3183 y3(i) = yy0(i,ic( 9,itq))
3184 y4(i) = yy0(i,ic(10,itq))
3185 z1(i) = zz0(i,ic( 7,itq))
3186 z2(i) = zz0(i,ic( 8,itq))
3187 z3(i) = zz0(i,ic( 9,itq))
3188 z4(i) = zz0(i,ic(10,itq))
3189
3190 vx1(i) = vx(i,ic( 7,itq))
3191 vx2(i) = vx(i,ic( 8,itq))
3192 vx3(i) = vx(i,ic( 9,itq))
3193 vx4(i) = vx(i,ic(10,itq))
3194 vy1(i) = vy(i,ic( 7,itq))
3195 vy2(i) = vy(i,ic( 8,itq))
3196 vy3(i) = vy(i,ic( 9,itq))
3197 vy4(i) = vy(i,ic(10,itq))
3198 vz1(i) = vz(i,ic( 7,itq))
3199 vz2(i) = vz(i,ic( 8,itq))
3200 vz3(i) = vz(i,ic( 9,itq))
3201 vz4(i) = vz(i,ic(10,itq))
3202 h1(i) = zero
3203 h2(i) = zero
3204 h3(i) = half
3205 h4(i) = half
3206
3207 n=cand_n(i)
3208 mglob = mseglo(cand_e(i))
3209 IF (itq==5 .or. itq==9 .or. itq==13 .or. itq==17)THEN
3210 mglob = mvoisin(1,cand_e(i))
3211
3212 ELSEIF(itq==6 .or. itq==10 .or. itq==14 .or. itq==18)THEN
3213 mglob = mvoisin(2,cand_e(i))
3214
3215 ELSEIF(itq==7 .or. itq==11 .or. itq==15 .or. itq==19)THEN
3216 mglob = mvoisin(3,cand_e(i))
3217
3218 ELSEIF(itq==8 .or. itq==12 .or. itq==16 .or. itq==20)THEN
3219 mglob = mvoisin(4,cand_e(i))
3220
3221 ENDIF
3222 IF (itq>8) itq = subtria_n(i)
3223
3224 IF(n <= nsn)THEN
3225 IF(irtlm(1,n)>0) THEN
3226 irtlm(1,n) = mglob
3227 time_s(n) = -ep20
3228 irtlm(2,n)= itq+20
3229 ELSEIF((time_s(n) < pene(i)) .or.
3230 . (time_s(n) == pene(i) .and.
3231 . -irtlm(1,n) < mglob ) )THEN
3232 irtlm(2,n)= itq+20
3233 irtlm(1,n) = -mglob
3234 time_s(n) = pene(i)
3235 ENDIF
3236 ELSE
3237 ii=n-nsn
3242 ELSEIF((
time_sfi(nin)%P(ii) < pene(i)).or.
3243 . (
time_sfi(nin)%P(ii)== pene(i).and.
3244 . -
irtlm_fi(nin)%P(1,ii) < mglob) )
THEN
3248 ENDIF
3249 ENDIF
3250 END DO
3251
3252
3253
3254
3255
3256
3257 DO i=1,jlt
3258 ce_loc(i) = cand_e(i)
3259 cn_loc(i) = cand_n(i)
3260 IF(pene(i) == 0 )THEN
3261
3262 n=cand_n(i)
3263 IF(n <= nsn)THEN
3264
3265
3266 IF(irtlm(1,n) > 0)THEN
3267 irtlm(1,n)=0
3268 time_s(n) = -ep20
3269 secnd_fr(4,n)= zero
3270 secnd_fr(5,n)= zero
3271 secnd_fr(6,n)= zero
3272 IF (imconv==1) pene_old(2,n) = zero
3273 pene_old(5,n) = zero
3274 IF (imconv==1) stif_old(2,n) = zero
3275 ENDIF
3276 ELSE
3277 IF(
irtlm_fi(nin)%P(1,n-nsn) > 0)
THEN
3283 IF (imconv==1)
pene_oldfi(nin)%P(2,n-nsn) = zero
3285 IF (imconv==1)
stif_oldfi(nin)%P(2,n-nsn) = zero
3286 ENDIF
3287 ENDIF
3288 ELSE
3289 ns=nsvg(i)
3290 IF (nm1(i)==zero.AND.nm2(i)==zero.AND.nm3(i)==zero) THEN
3291 nm1(i) = n1(i)
3292 nm2(i) = n2(i)
3293 nm3(i) = n3(i)
3294 END IF
3295
3296 IF (ix1(i) == ix2(i))THEN
3297 ELSEIF(ix2(i) == ix3(i))THEN
3298
3299 ELSEIF(ix4(i) == ix1(i))THEN
3300 ELSE
3301
3302 ENDIF
3303 xs = h1(i)*x1(i) + h2(i)*x2(i) + h3(i)*x3(i) + h4(i)*x4(i)
3304 ys = h1(i)*y1(i) + h2(i)*y2(i) + h3(i)*y3(i) + h4(i)*y4(i)
3305 zs = h1(i)*z1(i) + h2(i)*z2(i) + h3(i)*z3(i) + h4(i)*z4(i)
3306 xs = xs - xi(i)
3307 ys = ys - yi(i)
3308 zs = zs - zi(i)
3309
3310 aaa = xs*n1(i) + ys*n2(i) + zs*n3(i)
3311
3312 IF(aaa > zero)THEN
3313
3314 aaa = xs**2 + ys**2 + zs**2
3315 aaa = onep01*sqrt(aaa)
3316 pmax_gap =
max(pmax_gap,aaa)
3317
3318 n=cand_n(i)
3319 IF(n <= nsn)THEN
3320 pene_old(3,n) =
max(pene_old(3,n),aaa)
3321 ELSE
3324 ENDIF
3325 ENDIF
3326 ENDIF
3327 ENDDO
3328#include "lockoff.inc"
3329
3330
3331 IF (inacti==5) THEN
3332 IF (tt==zero) THEN
3333 DO i=1,jlt
3334 IF(pene(i) == zero)cycle
3335 jg = nsvg(i)
3336 n = cand_n(i)
3337 IF(jg > 0)THEN
3338#include "lockon.inc"
3339 pene_old(5,n) =
max( pene(i) ,pene_old(5,n) )
3340 stif_old(1,n) =
max( stif(i) ,stif_old(1,n) )
3341#include "lockoff.inc"
3342 ELSE
3343 jg = -jg
3344#include "lockon.inc"
3347#include "lockoff.inc"
3348 ENDIF
3349 ENDDO
3350
3351 ELSE
3352 DO i=1,jlt
3353 IF(pene(i) == zero)cycle
3354 jg = nsvg(i)
3355 n = cand_n(i)
3356 IF(jg > 0)THEN
3357 IF (ipen0(i)==0)THEN
3358#include "lockon.inc"
3359 pene_old(5,n) =
max( pene(i) ,pene_old(5,n) )
3360 stif_old(1,n) =
max( stif(i) ,stif_old(1,n) )
3361#include "lockoff.inc"
3362 END IF
3363 ELSE
3364 jg = -jg
3365 IF (ipen0(i)==0)THEN
3366#include "lockon.inc"
3369#include "lockoff.inc"
3370 END IF
3371 ENDIF
3372 ENDDO
3373 END IF
3374
3375 DO i=1,jlt
3376 IF(pene(i) == zero)cycle
3377
3378 n = cand_n(i)
3379 IF(jg > 0)THEN
3380 pene_sh =
max(zero,pene(i)-pene_old(5,n))
3381 ELSE
3382 jg = -jg
3384 ENDIF
3385
3386 pene(i) = pene_sh
3387 ENDDO
3388 END IF
3389
3390 penref(1:jlt) = zero
3391 IF (inacti==-1) THEN
3392 DO i=1,jlt
3393 IF(pene(i) == zero.OR.iknon(i)==0) cycle
3394 penref(i) = sqrt(eps2*
area(i))
3395 jg = nsvg(i)
3396 n = cand_n(i)
3397 IF(jg > 0)THEN
3398 pene_sh = em01*pene_old(5,n)
3399 ELSE
3400 jg = -jg
3402 ENDIF
3403 penref(i) =
max(penref(i),pene_sh)
3404 IF (iknon(i)<0) THEN
3405 f_pene = em01*pene(i)/pene_sh
3406 penref(i) = zero
3407 IF (f_pene>two_third) THEN
3408 fac_k =
min(twenty,six*f_pene)
3409 penref(i) = pene_sh/sqrt(fac_k)
3410 END IF
3411 END IF
3412 ENDDO
3413 ELSE
3414 DO i=1,jlt
3415 IF(pene(i) == zero.OR.iknon(i)==0) cycle
3416 n=cand_n(i)
3417 l=cand_e(i)
3418 penref(i) = one_fifth*sqrt
3419 gap2 = one_fifth*(gaps(i)+gap_m(l))
3420 IF (gap2 > em04) penref(i)=
min(penref(i),gap2)
3421 jg = nsvg(i)
3422 IF(jg > 0)THEN
3423 pene_sh = pene_old(5,n)
3424 ELSE
3425 jg = -jg
3427 ENDIF
3428 penref(i) =
max(penref(i),pene_sh)
3429 IF(iknon(i)==1) penref(i) = ten*penref(i)
3430 IF(iknon(i)>2) penref(i) = one_fifth*penref(i)
3431 ENDDO
3432 END IF
3433
3434 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