OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i15tot1.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com04_c.inc"
#include "units_c.inc"
#include "vectorize.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i15tot1 (noint, ndeb, nsc, x, ksurf, igrsurf, bufsf, ksc, ksi, nold, xp1, xp2, xp3, xp4, gx, xtk, ytk, ztk, ntx, nty, ntz, penet, depth, xi, yi, zi, nxi, nyi, nzi, ansmx, hold, iactiv, itab)

Function/Subroutine Documentation

◆ i15tot1()

subroutine i15tot1 ( integer noint,
integer ndeb,
integer nsc,
x,
integer ksurf,
type (surf_), dimension(nsurf) igrsurf,
bufsf,
integer, dimension(*) ksc,
integer, dimension(4,*) ksi,
nold,
xp1,
xp2,
xp3,
xp4,
gx,
xtk,
ytk,
ztk,
ntx,
nty,
ntz,
penet,
depth,
xi,
yi,
zi,
nxi,
nyi,
nzi,
ansmx,
hold,
integer, dimension(4,*) iactiv,
integer, dimension(*) itab )

Definition at line 30 of file i15tot1.F.

37C-----------------------------------------------
38C M o d u l e s
39C-----------------------------------------------
40 USE groupdef_mod
41C-----------------------------------------------
42C I m p l i c i t T y p e s
43C-----------------------------------------------
44#include "implicit_f.inc"
45#include "comlock.inc"
46C-----------------------------------------------
47C G l o b a l P a r a m e t e r s
48C-----------------------------------------------
49#include "mvsiz_p.inc"
50C-----------------------------------------------
51C C o m m o n B l o c k s
52C-----------------------------------------------
53#include "com04_c.inc"
54#include "units_c.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER KSURF ,KSI(4,*) ,
59 . NOINT ,NDEB, NSC, IACTIV(4,*), KSC(*), ITAB(*)
60C REAL
61 my_real
62 . x(3,*) ,bufsf(*) ,
63 . nold(3,*) ,xtk(4,*) ,ytk(4,*) ,ztk(4,*) ,
64 . ntx(4,*) ,nty(4,*) ,ntz(4,*) ,
65 . penet(4,*),depth(4,*),
66 . xi(4,*) ,yi(4,*) ,zi(4,*) ,
67 . nxi(4,*) ,nyi(4,*) ,nzi(4,*) ,
68 . xp1(3,mvsiz), xp2(3,mvsiz), xp3(3,mvsiz), xp4(3,mvsiz),
69 . gx(3,mvsiz), ansmx, hold(3,*)
70 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
71C-----------------------------------------------
72C L o c a l V a r i a b l e s
73C-----------------------------------------------
74 INTEGER ADRBUF, I, IL, NLS, IDG,
75 . IN1, IN2, IN3, IN4,
76 . INSIDE1,INSIDE2
77 my_real
78 . a, b, c, an, bn, cn, rot(9), dgr, expn,
79 . xg, yg, zg,
80 . ntn,
81 . n1, n2, n3, n, n11, n12, n13, n21, n22, n23, nr1, nr2,
82 . xkn1, ykn1, zkn1, sgnxkn, sgnykn, sgnzkn, xkn, ykn, zkn, eh,
83 . lambda1, lambda2, alp, bet,
84 . xh , yh , zh , mu, xh1, yh1, zh1, mu1, xh2, yh2, zh2, mu2,
85 . dx1, dy1, dz1, dx2, dy2, dz2, dx3, dy3, dz3,
86 . side1, side2, oldnx, oldny, oldnz, out, ps, nr,
87 . lambda, em, ep
88 my_real
89 . x1(mvsiz), y1(mvsiz), z1(mvsiz),
90 . x2(mvsiz), y2(mvsiz), z2(mvsiz),
91 . x3(mvsiz), y3(mvsiz), z3(mvsiz),
92 . x12(mvsiz), y12(mvsiz), z12(mvsiz),
93 . x13(mvsiz), y13(mvsiz), z13(mvsiz),
94 . x23(mvsiz), y23(mvsiz), z23(mvsiz),
95 . xm(mvsiz) , ym(mvsiz) , zm(mvsiz) ,
96 . xtk2(mvsiz), ytk2(mvsiz) , ztk2(mvsiz) ,
97 . d(mvsiz)
98C-----------------------------------------------
99 adrbuf=igrsurf(ksurf)%IAD_BUFR
100 a =bufsf(adrbuf+1)
101 b =bufsf(adrbuf+2)
102 c =bufsf(adrbuf+3)
103 dgr =bufsf(adrbuf+36)
104 idg =nint(dgr)
105 dgr =idg
106 dgr =max(two,dgr)
107 expn=dgr/(dgr-1.)
108 an=a**dgr
109 bn=b**dgr
110 cn=c**dgr
111 an=one/max(em20,an)
112 bn=one/max(em20,bn)
113 cn=one/max(em20,cn)
114 DO i=1,9
115 rot(i)=bufsf(adrbuf+7+i-1)
116 END DO
117C-------------------------------
118C Passage au repere local :
119C-------------------------------
120 DO 75 nls=ndeb+1,min(ndeb+mvsiz,nsc)
121 il =ksc(nls)
122 i =nls-ndeb
123 in1=ksi(1,il)
124 in2=ksi(2,il)
125 in3=ksi(3,il)
126 in4=ksi(4,il)
127C Passage au repere local.
128 xg=x(1,in1)-bufsf(adrbuf+4)
129 yg=x(2,in1)-bufsf(adrbuf+5)
130 zg=x(3,in1)-bufsf(adrbuf+6)
131 xp1(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
132 xp1(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
133 xp1(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
134 xg=x(1,in2)-bufsf(adrbuf+4)
135 yg=x(2,in2)-bufsf(adrbuf+5)
136 zg=x(3,in2)-bufsf(adrbuf+6)
137 xp2(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
138 xp2(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
139 xp2(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
140 xg=x(1,in3)-bufsf(adrbuf+4)
141 yg=x(2,in3)-bufsf(adrbuf+5)
142 zg=x(3,in3)-bufsf(adrbuf+6)
143 xp3(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
144 xp3(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
145 xp3(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
146 xg=x(1,in4)-bufsf(adrbuf+4)
147 yg=x(2,in4)-bufsf(adrbuf+5)
148 zg=x(3,in4)-bufsf(adrbuf+6)
149 xp4(1,i)=rot(1)*xg+rot(2)*yg+rot(3)*zg
150 xp4(2,i)=rot(4)*xg+rot(5)*yg+rot(6)*zg
151 xp4(3,i)=rot(7)*xg+rot(8)*yg+rot(9)*zg
152 gx(1,i)=fourth*(xp1(1,i)+xp2(1,i)+xp3(1,i)+xp4(1,i))
153 gx(2,i)=fourth*(xp1(2,i)+xp2(2,i)+xp3(2,i)+xp4(2,i))
154 gx(3,i)=fourth*(xp1(3,i)+xp2(3,i)+xp3(3,i)+xp4(3,i))
155 75 CONTINUE
156C-------------------------------
157C 1ER TRIANGLE.
158C-------------------------------
159 DO 100 nls=ndeb+1,min(ndeb+mvsiz,nsc)
160 il =ksc(nls)
161 i =nls-ndeb
162C-----
163 x1(i)=gx(1,i)
164 y1(i)=gx(2,i)
165 z1(i)=gx(3,i)
166 x2(i)=xp1(1,i)
167 y2(i)=xp1(2,i)
168 z2(i)=xp1(3,i)
169 x3(i)=xp2(1,i)
170 y3(i)=xp2(2,i)
171 z3(i)=xp2(3,i)
172 x12(i)=x2(i)-x1(i)
173 y12(i)=y2(i)-y1(i)
174 z12(i)=z2(i)-z1(i)
175 x13(i)=x3(i)-x1(i)
176 y13(i)=y3(i)-y1(i)
177 z13(i)=z3(i)-z1(i)
178 n1=y12(i)*z13(i)-z12(i)*y13(i)
179 n2=z12(i)*x13(i)-x12(i)*z13(i)
180 n3=x12(i)*y13(i)-y12(i)*x13(i)
181 ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
182 ntx(1,i)=ntn*n1
183 nty(1,i)=ntn*n2
184 ntz(1,i)=ntn*n3
185 d(i) =-ntx(1,i)*x1(i)-nty(1,i)*y1(i)-ntz(1,i)*z1(i)
186C-----
187 x23(i)=x3(i)-x2(i)
188 y23(i)=y3(i)-y2(i)
189 z23(i)=z3(i)-z2(i)
190 100 CONTINUE
191C-------------------------------
192C POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
193C-------------------------------
194 DO 125 nls=ndeb+1,min(ndeb+mvsiz,nsc)
195 il =ksc(nls)
196 i =nls-ndeb
197 eh =(abs(ntx(1,i)/(dgr*an))**expn)*an
198 . +(abs(nty(1,i)/(dgr*bn))**expn)*bn
199 . +(abs(ntz(1,i)/(dgr*cn))**expn)*cn
200C X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
201C LAMBDA1=EXP(LOG(1./EH)/EXPN)
202 lambda1=(max(em20,eh))**(-one/expn)
203 xh1=abs(lambda1*ntx(1,i)/(dgr*an))**(one/(dgr-one))
204 IF (ntx(1,i)<zero) xh1=-xh1
205 yh1=abs(lambda1*nty(1,i)/(dgr*bn))**(one/(dgr-one))
206 IF (nty(1,i)<zero) yh1=-yh1
207 zh1=abs(lambda1*ntz(1,i)/(dgr*cn))**(one/(dgr-one))
208 IF (ntz(1,i)<zero) zh1=-zh1
209 mu1 =-ntx(1,i)*xh1-nty(1,i)*yh1-ntz(1,i)*zh1-d(i)
210C LAMBDA2=-LAMBDA1
211 xh2 =-xh1
212 yh2 =-yh1
213 zh2 =-zh1
214C MU2 =-NTX(1,I)*XH2-NTY(1,I)*YH2-NTZ(1,I)*ZH2-D(I)
215 mu2=-mu1-two*d(i)
216 xtk(1,i) =xh1+mu1*ntx(1,i)
217 ytk(1,i) =yh1+mu1*nty(1,i)
218 ztk(1,i) =zh1+mu1*ntz(1,i)
219 xtk2(i) =xh2+mu2*ntx(1,i)
220 ytk2(i) =yh2+mu2*nty(1,i)
221 ztk2(i) =zh2+mu2*ntz(1,i)
222 125 CONTINUE
223C-------------------------------
224C RAMENER LE PT PK SUR LE TRIANGLE.
225C-------------------------------
226 DO 150 nls=ndeb+1,min(ndeb+mvsiz,nsc)
227 il =ksc(nls)
228 i =nls-ndeb
229C-----
230 dx1=xtk(1,i)-x1(i)
231 dy1=ytk(1,i)-y1(i)
232 dz1=ztk(1,i)-z1(i)
233 dx2=xtk(1,i)-x2(i)
234 dy2=ytk(1,i)-y2(i)
235 dz2=ztk(1,i)-z2(i)
236 out = (dy1*dz2-dy2*dz1)*ntx(1,i)
237 . +(dz1*dx2-dz2*dx1)*nty(1,i)
238 . +(dx1*dy2-dy1*dx2)*ntz(1,i)
239 IF (out<0.) THEN
240C PROJECTION SUR 1,2
241 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
242 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
243 bet=ps/max(em20,nr)
244 bet=max(bet,zero)
245 bet=min(bet,one)
246 alp=one-bet
247 xtk(1,i)=alp*x1(i)+bet*x2(i)
248 ytk(1,i)=alp*y1(i)+bet*y2(i)
249 ztk(1,i)=alp*z1(i)+bet*z2(i)
250 ENDIF
251 dx3=xtk(1,i)-x3(i)
252 dy3=ytk(1,i)-y3(i)
253 dz3=ztk(1,i)-z3(i)
254 out = (dy2*dz3-dy3*dz2)*ntx(1,i)
255 . +(dz2*dx3-dz3*dx2)*nty(1,i)
256 . +(dx2*dy3-dy2*dx3)*ntz(1,i)
257 IF (out<zero) THEN
258C PROJECTION SUR 2,3
259 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
260 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
261 bet=ps/max(em20,nr)
262 bet=max(bet,zero)
263 bet=min(bet,one)
264 alp=one-bet
265 xtk(1,i)=alp*x2(i)+bet*x3(i)
266 ytk(1,i)=alp*y2(i)+bet*y3(i)
267 ztk(1,i)=alp*z2(i)+bet*z3(i)
268 ENDIF
269 out = (dy3*dz1-dy1*dz3)*ntx(1,i)
270 . +(dz3*dx1-dz1*dx3)*nty(1,i)
271 . +(dx3*dy1-dy3*dx1)*ntz(1,i)
272 IF (out<zero) THEN
273C PROJECTION SUR 3,1
274 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
275 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
276 bet=ps/max(em20,nr)
277 bet=max(bet,zero)
278 bet=min(bet,one)
279 alp=one-bet
280 xtk(1,i)=alp*x3(i)+bet*x1(i)
281 ytk(1,i)=alp*y3(i)+bet*y1(i)
282 ztk(1,i)=alp*z3(i)+bet*z1(i)
283 ENDIF
284C-----
285 dx1=xtk2(i)-x1(i)
286 dy1=ytk2(i)-y1(i)
287 dz1=ztk2(i)-z1(i)
288 dx2=xtk2(i)-x2(i)
289 dy2=ytk2(i)-y2(i)
290 dz2=ztk2(i)-z2(i)
291 out = (dy1*dz2-dy2*dz1)*ntx(1,i)
292 . +(dz1*dx2-dz2*dx1)*nty(1,i)
293 . +(dx1*dy2-dy1*dx2)*ntz(1,i)
294 IF (out<zero) THEN
295C PROJECTION SUR 1,2
296 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
297 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
298 bet=ps/max(em20,nr)
299 bet=max(bet,zero)
300 bet=min(bet,one)
301 alp=one-bet
302 xtk2(i)=alp*x1(i)+bet*x2(i)
303 ytk2(i)=alp*y1(i)+bet*y2(i)
304 ztk2(i)=alp*z1(i)+bet*z2(i)
305 ENDIF
306 dx3=xtk2(i)-x3(i)
307 dy3=ytk2(i)-y3(i)
308 dz3=ztk2(i)-z3(i)
309 out = (dy2*dz3-dy3*dz2)*ntx(1,i)
310 . +(dz2*dx3-dz3*dx2)*nty(1,i)
311 . +(dx2*dy3-dy2*dx3)*ntz(1,i)
312 IF (out<zero) THEN
313C PROJECTION SUR 2,3
314 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
315 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
316 bet=ps/max(em20,nr)
317 bet=max(bet,zero)
318 bet=min(bet,one)
319 alp=1.-bet
320 xtk2(i)=alp*x2(i)+bet*x3(i)
321 ytk2(i)=alp*y2(i)+bet*y3(i)
322 ztk2(i)=alp*z2(i)+bet*z3(i)
323 ENDIF
324 out = (dy3*dz1-dy1*dz3)*ntx(1,i)
325 . +(dz3*dx1-dz1*dx3)*nty(1,i)
326 . +(dx3*dy1-dy3*dx1)*ntz(1,i)
327 IF (out<zero) THEN
328C PROJECTION SUR 3,1
329 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
330 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
331 bet=ps/max(em20,nr)
332 bet=max(bet,zero)
333 bet=min(bet,one)
334 alp=1.-bet
335 xtk2(i)=alp*x3(i)+bet*x1(i)
336 ytk2(i)=alp*y3(i)+bet*y1(i)
337 ztk2(i)=alp*z3(i)+bet*z1(i)
338 ENDIF
339C-------------------------------
340 150 CONTINUE
341C-------------------------------
342C PROJECTION DE PK SUR L ET PENETRATION.
343C-------------------------------
344 DO 175 nls=ndeb+1,min(ndeb+mvsiz,nsc)
345 il =ksc(nls)
346 i =nls-ndeb
347 IF (iactiv(1,il)==-1) GOTO 175
348C-----
349 xh =xtk(1,i)
350 yh =ytk(1,i)
351 zh =ztk(1,i)
352 xkn1 =xh**(idg-1)
353 sgnxkn=-one
354 IF (xkn1*xh>=zero) sgnxkn=one
355 ykn1 =yh**(idg-1)
356 sgnykn=-one
357 IF (ykn1*yh>=zero) sgnykn=one
358 zkn1 =zh**(idg-1)
359 sgnzkn=-one
360 IF (zkn1*zh>=zero) sgnzkn=one
361 n11 =sgnxkn*xkn1*an
362 n12 =sgnykn*ykn1*bn
363 n13 =sgnzkn*zkn1*cn
364 nr1 =n11*n11+n12*n12+n13*n13
365 nr1 =sqrt(nr1)
366 em =n11*xtk(1,i)+n12*ytk(1,i)+n13*ztk(1,i)
367 IF (em<=one) THEN
368 lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
369 . / max(exp((dgr-one)*log(em20)/dgr),nr1)
370 inside1=1
371 ELSE
372 inside1=0
373 ENDIF
374C-----
375 xh =xtk2(i)
376 yh =ytk2(i)
377 zh =ztk2(i)
378 xkn1 =xh**(idg-1)
379 sgnxkn=-one
380 IF (xkn1*xh>=zero) sgnxkn=one
381 ykn1 =yh**(idg-1)
382 sgnykn=-one
383 IF (ykn1*yh>=zero) sgnykn=one
384 zkn1 =zh**(idg-1)
385 sgnzkn=-one
386 IF (zkn1*zh>=zero) sgnzkn=one
387 n21 =sgnxkn*xkn1*an
388 n22 =sgnykn*ykn1*bn
389 n23 =sgnzkn*zkn1*cn
390 nr2 =n21*n21+n22*n22+n23*n23
391 nr2 =sqrt(nr2)
392 em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
393 IF (em<=one) THEN
394 lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
395 . / max(exp((dgr-one)*log(em20)/dgr),nr2)
396 inside2=1
397 ELSE
398 inside2=0
399 ENDIF
400C-----
401 IF (inside1==0.AND.inside2==0) THEN
402 iactiv(1,il)=0
403 ELSE
404C-----
405 IF (iactiv(1,il)==0) THEN
406 IF (inside1/=0.AND.inside2/=0) THEN
407 IF (abs(lambda1)>=abs(lambda2)) THEN
408 xm(i)=xtk(1,i)-lambda1*n11/max(em20,nr1)
409 ym(i)=ytk(1,i)-lambda1*n12/max(em20,nr1)
410 zm(i)=ztk(1,i)-lambda1*n13/max(em20,nr1)
411 ELSE
412 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
413 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
414 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
415 xtk(1,i)=xtk2(i)
416 ytk(1,i)=ytk2(i)
417 ztk(1,i)=ztk2(i)
418 ENDIF
419 ELSEIF(inside1/=0) THEN
420 xm(i)=xtk(1,i)-lambda1*n11/max(em20,nr1)
421 ym(i)=ytk(1,i)-lambda1*n12/max(em20,nr1)
422 zm(i)=ztk(1,i)-lambda1*n13/max(em20,nr1)
423 ELSEIF(inside2/=0) THEN
424 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
425 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
426 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
427 xtk(1,i)=xtk2(i)
428 ytk(1,i)=ytk2(i)
429 ztk(1,i)=ztk2(i)
430 ENDIF
431C-----
432 ELSE
433 xh=hold(1,4*(il-1)+1)
434 yh=hold(2,4*(il-1)+1)
435 zh=hold(3,4*(il-1)+1)
436 n1=nold(1,4*(il-1)+1)
437 n2=nold(2,4*(il-1)+1)
438 n3=nold(3,4*(il-1)+1)
439 lambda1=(xh-xtk(1,i))*n1
440 . +(yh-ytk(1,i))*n2
441 . +(zh-ztk(1,i))*n3
442 lambda2=(xh-xtk2(i))*n1
443 . +(yh-ytk2(i))*n2
444 . +(zh-ztk2(i))*n3
445 IF (inside1/=0.AND.inside2/=0) THEN
446 IF (abs(lambda1)>=abs(lambda2)) THEN
447 xm(i)=xtk(1,i)+lambda1*n1
448 ym(i)=ytk(1,i)+lambda1*n2
449 zm(i)=ztk(1,i)+lambda1*n3
450 ELSE
451 xm(i)=xtk2(i)+lambda2*n1
452 ym(i)=ytk2(i)+lambda2*n2
453 zm(i)=ztk2(i)+lambda2*n3
454 xtk(1,i)=xtk2(i)
455 ytk(1,i)=ytk2(i)
456 ztk(1,i)=ztk2(i)
457 ENDIF
458 ELSEIF(inside1/=0) THEN
459 xm(i)=xtk(1,i)+lambda1*n1
460 ym(i)=ytk(1,i)+lambda1*n2
461 zm(i)=ztk(1,i)+lambda1*n3
462 ELSEIF(inside2/=0) THEN
463 xm(i)=xtk2(i)+lambda2*n1
464 ym(i)=ytk2(i)+lambda2*n2
465 zm(i)=ztk2(i)+lambda2*n3
466 xtk(1,i)=xtk2(i)
467 ytk(1,i)=ytk2(i)
468 ztk(1,i)=ztk2(i)
469 ENDIF
470C one more iteration.
471 xkn1 =xm(i)**(idg-1)
472 sgnxkn=-one
473 IF (xkn1*xm(i)>=zero) sgnxkn=one
474 ykn1 =ym(i)**(idg-1)
475 sgnykn=-one
476 IF (ykn1*ym(i)>=zero) sgnykn=one
477 zkn1 =zm(i)**(idg-1)
478 sgnzkn=-one
479 IF (zkn1*zm(i)>=zero) sgnzkn=one
480 n1 =sgnxkn*xkn1*an
481 n2 =sgnykn*ykn1*bn
482 n3 =sgnzkn*zkn1*cn
483 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
484 xm(i)=xm(i)/max(em20,em**(one/dgr))
485 ym(i)=ym(i)/max(em20,em**(one/dgr))
486 zm(i)=zm(i)/max(em20,em**(one/dgr))
487 xkn1 =xm(i)**(idg-1)
488 sgnxkn=-one
489 IF (xkn1*xm(i)>=zero) sgnxkn=one
490 ykn1 =ym(i)**(idg-1)
491 sgnykn=-one
492 IF (ykn1*ym(i)>=zero) sgnykn=one
493 zkn1 =zm(i)**(idg-1)
494 sgnzkn=-one
495 IF (zkn1*zm(i)>=zero) sgnzkn=one
496 n1 =sgnxkn*xkn1*an
497 n2 =sgnykn*ykn1*bn
498 n3 =sgnzkn*zkn1*cn
499 nr =n1*n1+n2*n2+n3*n3
500 nr =one/max(em20,sqrt(nr))
501 n1 =n1*nr
502 n2 =n2*nr
503 n3 =n3*nr
504 lambda1=(xm(i)-xtk(1,i))*n1
505 . +(ym(i)-ytk(1,i))*n2
506 . +(zm(i)-ztk(1,i))*n3
507 xm(i)=xtk(1,i)+lambda1*n1
508 ym(i)=ytk(1,i)+lambda1*n2
509 zm(i)=ztk(1,i)+lambda1*n3
510 ENDIF
511 iactiv(1,il)=iactiv(1,il)+1
512 ENDIF
513 175 CONTINUE
514C-----
515#include "vectorize.inc"
516 DO 200 nls=ndeb+1,min(ndeb+mvsiz,nsc)
517 il =ksc(nls)
518 i =nls-ndeb
519C-----
520 IF (iactiv(1,il)<=0) GOTO 200
521C-----
522C projection radiale de M.
523 xkn1 =xm(i)**(idg-1)
524 sgnxkn=-one
525 IF (xkn1*xm(i)>=zero) sgnxkn=one
526 ykn1 =ym(i)**(idg-1)
527 sgnykn=-one
528 IF (ykn1*ym(i)>=zero) sgnykn=one
529 zkn1 =zm(i)**(idg-1)
530 sgnzkn=-one
531 IF (zkn1*zm(i)>=zero) sgnzkn=one
532 n1 =sgnxkn*xkn1*an
533 n2 =sgnykn*ykn1*bn
534 n3 =sgnzkn*zkn1*cn
535 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
536C------
537 xm(i)=xm(i)/max(em20,em**(one/dgr))
538 ym(i)=ym(i)/max(em20,em**(one/dgr))
539 zm(i)=zm(i)/max(em20,em**(one/dgr))
540C normale a l'ellipsoide en M.
541 xkn1 =xm(i)**(idg-1)
542 sgnxkn=-one
543 IF (xkn1*xm(i)>=zero) sgnxkn=one
544 ykn1 =ym(i)**(idg-1)
545 sgnykn=-one
546 IF (ykn1*ym(i)>=zero) sgnykn=one
547 zkn1 =zm(i)**(idg-1)
548 sgnzkn=-one
549 IF (zkn1*zm(i)>=zero) sgnzkn=one
550 n1 =sgnxkn*xkn1*an
551 n2 =sgnykn*ykn1*bn
552 n3 =sgnzkn*zkn1*cn
553 nr =n1*n1+n2*n2+n3*n3
554 nr =one/max(em20,sqrt(nr))
555 nxi(1,i)=n1*nr
556 nyi(1,i)=n2*nr
557 nzi(1,i)=n3*nr
558 xi(1,i)=xm(i)
559 yi(1,i)=ym(i)
560 zi(1,i)=zm(i)
561 depth(1,i)=xm(i)*nxi(1,i)+ym(i)*nyi(1,i)+zm(i)*nzi(1,i)
562C------
563 penet(1,i)=(xtk(1,i)-xm(i))*nxi(1,i)
564 . +(ytk(1,i)-ym(i))*nyi(1,i)
565 . +(ztk(1,i)-zm(i))*nzi(1,i)
566 penet(1,i)=-penet(1,i)
567 penet(1,i)=max(zero,penet(1,i))
568 IF (depth(1,i)-penet(1,i)<em10*depth(1,i)) THEN
569 penet(1,i)=zero
570 iactiv(1,il)=-2
571 ENDIF
572 ansmx=max(ansmx,penet(1,i))
573 200 CONTINUE
574C-------------------------------
575C MESSAGES DESACTIVATION
576C-------------------------------
577 DO 210 nls=ndeb+1,min(ndeb+mvsiz,nsc)
578 il =ksc(nls)
579 i =nls-ndeb
580C-----
581 IF (iactiv(1,il)==-2) THEN
582 in1=itab(ksi(1,il))
583 in2=itab(ksi(2,il))
584 in3=itab(ksi(3,il))
585 in4=itab(ksi(4,il))
586#include "lockon.inc"
587 WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
588 WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
589 . ' FROM 4NODES ELEMENT/SEGMENT,',
590 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
591 . in1,in2,'ISOCENTER(',in1,in2,in3,in4,')'
592 WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
593 WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
594 . ' FROM 4NODES ELEMENT/SEGMENT,',
595 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
596 . in1,in2,'ISOCENTER(',in1,in2,in3,in4,')'
597#include "lockoff.inc"
598 iactiv(1,il)=-1
599 ENDIF
600 210 CONTINUE
601C-------------------------------
602C 2EME TRIANGLE.
603C-------------------------------
604 DO 225 nls=ndeb+1,min(ndeb+mvsiz,nsc)
605 il =ksc(nls)
606 i =nls-ndeb
607C-----
608 x1(i)=gx(1,i)
609 y1(i)=gx(2,i)
610 z1(i)=gx(3,i)
611 x2(i)=xp2(1,i)
612 y2(i)=xp2(2,i)
613 z2(i)=xp2(3,i)
614 x3(i)=xp3(1,i)
615 y3(i)=xp3(2,i)
616 z3(i)=xp3(3,i)
617 x12(i)=x2(i)-x1(i)
618 y12(i)=y2(i)-y1(i)
619 z12(i)=z2(i)-z1(i)
620 x13(i)=x3(i)-x1(i)
621 y13(i)=y3(i)-y1(i)
622 z13(i)=z3(i)-z1(i)
623 n1=y12(i)*z13(i)-z12(i)*y13(i)
624 n2=z12(i)*x13(i)-x12(i)*z13(i)
625 n3=x12(i)*y13(i)-y12(i)*x13(i)
626 ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
627 ntx(2,i)=ntn*n1
628 nty(2,i)=ntn*n2
629 ntz(2,i)=ntn*n3
630 d(i) =-ntx(2,i)*x1(i)-nty(2,i)*y1(i)-ntz(2,i)*z1(i)
631C-----
632 x23(i)=x3(i)-x2(i)
633 y23(i)=y3(i)-y2(i)
634 z23(i)=z3(i)-z2(i)
635 225 CONTINUE
636C-------------------------------
637C POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
638C-------------------------------
639 DO 250 nls=ndeb+1,min(ndeb+mvsiz,nsc)
640 il =ksc(nls)
641 i =nls-ndeb
642C-----
643 eh =(abs(ntx(2,i)/(dgr*an))**expn)*an
644 . +(abs(nty(2,i)/(dgr*bn))**expn)*bn
645 . +(abs(ntz(2,i)/(dgr*cn))**expn)*cn
646C X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
647C LAMBDA1=EXP(LOG(1./EH)/EXPN)
648 lambda1=(max(em20,eh))**(-one/expn)
649 xh1 =abs(lambda1*ntx(2,i)/(dgr*an))**(one/(dgr-one))
650 IF (ntx(2,i)<zero) xh1=-xh1
651 yh1 =abs(lambda1*nty(2,i)/(dgr*bn))**(one/(dgr-one))
652 IF (nty(2,i)<zero) yh1=-yh1
653 zh1 =abs(lambda1*ntz(2,i)/(dgr*cn))**(one/(dgr-one))
654 IF (ntz(2,i)<zero) zh1=-zh1
655 mu1 =-ntx(2,i)*xh1-nty(2,i)*yh1-ntz(2,i)*zh1-d(i)
656C LAMBDA2=-LAMBDA1
657 xh2 =-xh1
658 yh2 =-yh1
659 zh2 =-zh1
660C MU2 =-NTX(2,I)*XH2-NTY(2,I)*YH2-NTZ(2,I)*ZH2-D(I)
661 mu2=-mu1-two*d(i)
662 xtk(2,i) =xh1+mu1*ntx(2,i)
663 ytk(2,i) =yh1+mu1*nty(2,i)
664 ztk(2,i) =zh1+mu1*ntz(2,i)
665 xtk2(i) =xh2+mu2*ntx(2,i)
666 ytk2(i) =yh2+mu2*nty(2,i)
667 ztk2(i) =zh2+mu2*ntz(2,i)
668C-------------------------------
669 250 CONTINUE
670C-------------------------------
671C RAMENER LE PT PK SUR LE TRIANGLE.
672C-------------------------------
673 DO 275 nls=ndeb+1,min(ndeb+mvsiz,nsc)
674 il =ksc(nls)
675 i =nls-ndeb
676C-----
677 dx1=xtk(2,i)-x1(i)
678 dy1=ytk(2,i)-y1(i)
679 dz1=ztk(2,i)-z1(i)
680 dx2=xtk(2,i)-x2(i)
681 dy2=ytk(2,i)-y2(i)
682 dz2=ztk(2,i)-z2(i)
683 out = (dy1*dz2-dy2*dz1)*ntx(2,i)
684 . +(dz1*dx2-dz2*dx1)*nty(2,i)
685 . +(dx1*dy2-dy1*dx2)*ntz(2,i)
686 IF (out<zero) THEN
687C PROJECTION SUR 1,2
688 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
689 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
690 bet=ps/max(em20,nr)
691 bet=max(bet,zero)
692 bet=min(bet,one)
693 alp=1.-bet
694 xtk(2,i)=alp*x1(i)+bet*x2(i)
695 ytk(2,i)=alp*y1(i)+bet*y2(i)
696 ztk(2,i)=alp*z1(i)+bet*z2(i)
697 ENDIF
698 dx3=xtk(2,i)-x3(i)
699 dy3=ytk(2,i)-y3(i)
700 dz3=ztk(2,i)-z3(i)
701 out = (dy2*dz3-dy3*dz2)*ntx(2,i)
702 . +(dz2*dx3-dz3*dx2)*nty(2,i)
703 . +(dx2*dy3-dy2*dx3)*ntz(2,i)
704 IF (out<zero) THEN
705C PROJECTION SUR 2,3
706 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
707 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
708 bet=ps/max(em20,nr)
709 bet=max(bet,zero)
710 bet=min(bet,one)
711 alp=one-bet
712 xtk(2,i)=alp*x2(i)+bet*x3(i)
713 ytk(2,i)=alp*y2(i)+bet*y3(i)
714 ztk(2,i)=alp*z2(i)+bet*z3(i)
715 ENDIF
716 out = (dy3*dz1-dy1*dz3)*ntx(2,i)
717 . +(dz3*dx1-dz1*dx3)*nty(2,i)
718 . +(dx3*dy1-dy3*dx1)*ntz(2,i)
719 IF (out<zero) THEN
720C PROJECTION SUR 3,1
721 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
722 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
723 bet=ps/max(em20,nr)
724 bet=max(bet,zero)
725 bet=min(bet,one)
726 alp=1.-bet
727 xtk(2,i)=alp*x3(i)+bet*x1(i)
728 ytk(2,i)=alp*y3(i)+bet*y1(i)
729 ztk(2,i)=alp*z3(i)+bet*z1(i)
730 ENDIF
731C-----
732 dx1=xtk2(i)-x1(i)
733 dy1=ytk2(i)-y1(i)
734 dz1=ztk2(i)-z1(i)
735 dx2=xtk2(i)-x2(i)
736 dy2=ytk2(i)-y2(i)
737 dz2=ztk2(i)-z2(i)
738 out = (dy1*dz2-dy2*dz1)*ntx(2,i)
739 . +(dz1*dx2-dz2*dx1)*nty(2,i)
740 . +(dx1*dy2-dy1*dx2)*ntz(2,i)
741 IF (out<zero) THEN
742C PROJECTION SUR 1,2
743 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
744 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
745 bet=ps/max(em20,nr)
746 bet=max(bet,zero)
747 bet=min(bet,one)
748 alp=1.-bet
749 xtk2(i)=alp*x1(i)+bet*x2(i)
750 ytk2(i)=alp*y1(i)+bet*y2(i)
751 ztk2(i)=alp*z1(i)+bet*z2(i)
752 ENDIF
753 dx3=xtk2(i)-x3(i)
754 dy3=ytk2(i)-y3(i)
755 dz3=ztk2(i)-z3(i)
756 out = (dy2*dz3-dy3*dz2)*ntx(2,i)
757 . +(dz2*dx3-dz3*dx2)*nty(2,i)
758 . +(dx2*dy3-dy2*dx3)*ntz(2,i)
759 IF (out<zero) THEN
760C PROJECTION SUR 2,3
761 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
762 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
763 bet=ps/max(em20,nr)
764 bet=max(bet,zero)
765 bet=min(bet,one)
766 alp=one-bet
767 xtk2(i)=alp*x2(i)+bet*x3(i)
768 ytk2(i)=alp*y2(i)+bet*y3(i)
769 ztk2(i)=alp*z2(i)+bet*z3(i)
770 ENDIF
771 out = (dy3*dz1-dy1*dz3)*ntx(2,i)
772 . +(dz3*dx1-dz1*dx3)*nty(2,i)
773 . +(dx3*dy1-dy3*dx1)*ntz(2,i)
774 IF (out<zero) THEN
775C PROJECTION SUR 3,1
776 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
777 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
778 bet=ps/max(em20,nr)
779 bet=max(bet,zero)
780 bet=min(bet,one)
781 alp=one-bet
782 xtk2(i)=alp*x3(i)+bet*x1(i)
783 ytk2(i)=alp*y3(i)+bet*y1(i)
784 ztk2(i)=alp*z3(i)+bet*z1(i)
785 ENDIF
786C-------------------------------
787 275 CONTINUE
788C-------------------------------
789C PROJECTION DE PK SUR L ET PENETRATION.
790C-------------------------------
791 DO 300 nls=ndeb+1,min(ndeb+mvsiz,nsc)
792 il =ksc(nls)
793 i =nls-ndeb
794 IF (iactiv(2,il)==-1) GOTO 300
795C-----
796 xh =xtk(2,i)
797 yh =ytk(2,i)
798 zh =ztk(2,i)
799 xkn1 =xh**(idg-1)
800 sgnxkn=-one
801 IF (xkn1*xh>=zero) sgnxkn=one
802 ykn1 =yh**(idg-1)
803 sgnykn=-1.
804 IF (ykn1*yh>=0.) sgnykn=1.
805 zkn1 =zh**(idg-1)
806 sgnzkn=-one
807 IF (zkn1*zh>=zero) sgnzkn=one
808 n11 =sgnxkn*xkn1*an
809 n12 =sgnykn*ykn1*bn
810 n13 =sgnzkn*zkn1*cn
811 nr1 =n11*n11+n12*n12+n13*n13
812 nr1 =sqrt(nr1)
813 em =n11*xtk(2,i)+n12*ytk(2,i)+n13*ztk(2,i)
814 IF (em<=one) THEN
815 lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
816 . / max(exp((dgr-one)*log(em20)/dgr),nr1)
817 inside1=1
818 ELSE
819 inside1=0
820 ENDIF
821C-----
822 xh =xtk2(i)
823 yh =ytk2(i)
824 zh =ztk2(i)
825 xkn1 =xh**(idg-1)
826 sgnxkn=-one
827 IF (xkn1*xh>=zero) sgnxkn=one
828 ykn1 =yh**(idg-1)
829 sgnykn=-one
830 IF (ykn1*yh>=zero) sgnykn=one
831 zkn1 =zh**(idg-1)
832 sgnzkn=-one
833 IF (zkn1*zh>=zero) sgnzkn=one
834 n21 =sgnxkn*xkn1*an
835 n22 =sgnykn*ykn1*bn
836 n23 =sgnzkn*zkn1*cn
837 nr2 =n21*n21+n22*n22+n23*n23
838 nr2 =sqrt(nr2)
839 em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
840 IF (em<=one) THEN
841 lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
842 . / max(exp((dgr-one)*log(em20)/dgr),nr2)
843 inside2=1
844 ELSE
845 inside2=0
846 ENDIF
847C-----
848 IF (inside1==0.AND.inside2==0) THEN
849 iactiv(2,il)=0
850 ELSE
851C-----
852 IF (iactiv(2,il)==0) THEN
853 IF (inside1/=0.AND.inside2/=0) THEN
854 IF (abs(lambda1)>=abs(lambda2)) THEN
855 xm(i)=xtk(2,i)-lambda1*n11/max(em20,nr1)
856 ym(i)=ytk(2,i)-lambda1*n12/max(em20,nr1)
857 zm(i)=ztk(2,i)-lambda1*n13/max(em20,nr1)
858 ELSE
859 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
860 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
861 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
862 xtk(2,i)=xtk2(i)
863 ytk(2,i)=ytk2(i)
864 ztk(2,i)=ztk2(i)
865 ENDIF
866 ELSEIF(inside1/=0) THEN
867 xm(i)=xtk(2,i)-lambda1*n11/max(em20,nr1)
868 ym(i)=ytk(2,i)-lambda1*n12/max(em20,nr1)
869 zm(i)=ztk(2,i)-lambda1*n13/max(em20,nr1)
870 ELSEIF(inside2/=0) THEN
871 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
872 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
873 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
874 xtk(2,i)=xtk2(i)
875 ytk(2,i)=ytk2(i)
876 ztk(2,i)=ztk2(i)
877 ENDIF
878C-----
879 ELSE
880 xh=hold(1,4*(il-1)+2)
881 yh=hold(2,4*(il-1)+2)
882 zh=hold(3,4*(il-1)+2)
883 n1=nold(1,4*(il-1)+2)
884 n2=nold(2,4*(il-1)+2)
885 n3=nold(3,4*(il-1)+2)
886 lambda1=(xh-xtk(2,i))*n1
887 . +(yh-ytk(2,i))*n2
888 . +(zh-ztk(2,i))*n3
889 lambda2=(xh-xtk2(i))*n1
890 . +(yh-ytk2(i))*n2
891 . +(zh-ztk2(i))*n3
892 IF (inside1/=0.AND.inside2/=0) THEN
893 IF (abs(lambda1)>=abs(lambda2)) THEN
894 xm(i)=xtk(2,i)+lambda1*n1
895 ym(i)=ytk(2,i)+lambda1*n2
896 zm(i)=ztk(2,i)+lambda1*n3
897 ELSE
898 xm(i)=xtk2(i)+lambda2*n1
899 ym(i)=ytk2(i)+lambda2*n2
900 zm(i)=ztk2(i)+lambda2*n3
901 xtk(2,i)=xtk2(i)
902 ytk(2,i)=ytk2(i)
903 ztk(2,i)=ztk2(i)
904 ENDIF
905 ELSEIF(inside1/=0) THEN
906 xm(i)=xtk(2,i)+lambda1*n1
907 ym(i)=ytk(2,i)+lambda1*n2
908 zm(i)=ztk(2,i)+lambda1*n3
909 ELSEIF(inside2/=0) THEN
910 xm(i)=xtk2(i)+lambda2*n1
911 ym(i)=ytk2(i)+lambda2*n2
912 zm(i)=ztk2(i)+lambda2*n3
913 xtk(2,i)=xtk2(i)
914 ytk(2,i)=ytk2(i)
915 ztk(2,i)=ztk2(i)
916 ENDIF
917C one more iteration.
918 xkn1 =xm(i)**(idg-1)
919 sgnxkn=-one
920 IF (xkn1*xm(i)>=zero) sgnxkn=one
921 sgnykn=-one
922 IF (ykn1*ym(i)>=zero) sgnykn=one
923 zkn1 =zm(i)**(idg-1)
924 sgnzkn=-one
925 IF (zkn1*zm(i)>=zero) sgnzkn=one
926 n1 =sgnxkn*xkn1*an
927 n2 =sgnykn*ykn1*bn
928 n3 =sgnzkn*zkn1*cn
929 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
930 xm(i)=xm(i)/max(em20,em**(one/dgr))
931 ym(i)=ym(i)/max(em20,em**(one/dgr))
932 zm(i)=zm(i)/max(em20,em**(one/dgr))
933 xkn1 =xm(i)**(idg-1)
934 sgnxkn=-one
935 IF (xkn1*xm(i)>=zero) sgnxkn=one
936 ykn1 =ym(i)**(idg-1)
937 sgnykn=-one
938 IF (ykn1*ym(i)>=zero) sgnykn=one
939 zkn1 =zm(i)**(idg-1)
940 sgnzkn=-one
941 IF (zkn1*zm(i)>=zero) sgnzkn=one
942 n1 =sgnxkn*xkn1*an
943 n2 =sgnykn*ykn1*bn
944 n3 =sgnzkn*zkn1*cn
945 nr =n1*n1+n2*n2+n3*n3
946 nr =one/max(em20,sqrt(nr))
947 n1 =n1*nr
948 n2 =n2*nr
949 n3 =n3*nr
950 lambda1=(xm(i)-xtk(2,i))*n1
951 . +(ym(i)-ytk(2,i))*n2
952 . +(zm(i)-ztk(2,i))*n3
953 xm(i)=xtk(2,i)+lambda1*n1
954 ym(i)=ytk(2,i)+lambda1*n2
955 zm(i)=ztk(2,i)+lambda1*n3
956 ENDIF
957 iactiv(2,il)=iactiv(2,il)+1
958 ENDIF
959 300 CONTINUE
960C-----
961#include "vectorize.inc"
962 DO 325 nls=ndeb+1,min(ndeb+mvsiz,nsc)
963 il =ksc(nls)
964 i =nls-ndeb
965C-----
966 IF (iactiv(2,il)<=0) GOTO 325
967C-----
968C projection radiale de M.
969 xkn1 =xm(i)**(idg-1)
970 sgnxkn=-one
971 IF (xkn1*xm(i)>=zero) sgnxkn=one
972 ykn1 =ym(i)**(idg-1)
973 sgnykn=-one
974 IF (ykn1*ym(i)>=zero) sgnykn=one
975 zkn1 =zm(i)**(idg-1)
976 sgnzkn=-one
977 IF (zkn1*zm(i)>=zero) sgnzkn=one
978 n1 =sgnxkn*xkn1*an
979 n2 =sgnykn*ykn1*bn
980 n3 =sgnzkn*zkn1*cn
981 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
982C------
983 xm(i)=xm(i)/max(em20,em**(one/dgr))
984 ym(i)=ym(i)/max(em20,em**(one/dgr))
985 zm(i)=zm(i)/max(em20,em**(one/dgr))
986C------
987C normale a l'ellipsoide en M.
988 xkn1 =xm(i)**(idg-1)
989 sgnxkn=-one
990 IF (xkn1*xm(i)>=zero) sgnxkn=one
991 ykn1 =ym(i)**(idg-1)
992 sgnykn=-one
993 IF (ykn1*ym(i)>=zero) sgnykn=one
994 zkn1 =zm(i)**(idg-1)
995 sgnzkn=-one
996 IF (zkn1*zm(i)>=zero) sgnzkn=one
997 n1 =sgnxkn*xkn1*an
998 n2 =sgnykn*ykn1*bn
999 n3 =sgnzkn*zkn1*cn
1000 nr =n1*n1+n2*n2+n3*n3
1001 nr =one/max(em20,sqrt(nr))
1002 nxi(2,i)=n1*nr
1003 nyi(2,i)=n2*nr
1004 nzi(2,i)=n3*nr
1005 xi(2,i)=xm(i)
1006 yi(2,i)=ym(i)
1007 zi(2,i)=zm(i)
1008 depth(2,i)=xm(i)*nxi(2,i)+ym(i)*nyi(2,i)+zm(i)*nzi(2,i)
1009C------
1010 penet(2,i)=(xtk(2,i)-xm(i))*nxi(2,i)
1011 . +(ytk(2,i)-ym(i))*nyi(2,i)
1012 . +(ztk(2,i)-zm(i))*nzi(2,i)
1013 penet(2,i)=-penet(2,i)
1014 penet(2,i)=max(zero,penet(2,i))
1015 IF (depth(2,i)-penet(2,i)<em10*depth(2,i)) THEN
1016 penet(2,i)=zero
1017 iactiv(2,il)=-2
1018 ENDIF
1019 ansmx=max(ansmx,penet(2,i))
1020 325 CONTINUE
1021C-------------------------------
1022C MESSAGES DESACTIVATION
1023C-------------------------------
1024 DO 335 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1025 il =ksc(nls)
1026 i =nls-ndeb
1027C-----
1028 IF (iactiv(2,il)==-2) THEN
1029 in1=itab(ksi(1,il))
1030 in2=itab(ksi(2,il))
1031 in3=itab(ksi(3,il))
1032 in4=itab(ksi(4,il))
1033#include "lockon.inc"
1034 WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
1035 WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
1036 . ' FROM 4NODES ELEMENT/SEGMENT,',
1037 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
1038 . in2,in3,'ISOCENTER(',in1,in2,in3,in4,')'
1039 WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
1040 WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
1041 . ' FROM 4NODES ELEMENT/SEGMENT,',
1042 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
1043 . in2,in3,'ISOCENTER(',in1,in2,in3,in4,')'
1044#include "lockoff.inc"
1045 iactiv(2,il)=-1
1046 ENDIF
1047 335 CONTINUE
1048C-------------------------------
1049C 3EME TRIANGLE.
1050C-------------------------------
1051 DO 350 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1052 il =ksc(nls)
1053 i =nls-ndeb
1054C-----
1055 x1(i)=gx(1,i)
1056 y1(i)=gx(2,i)
1057 z1(i)=gx(3,i)
1058 x2(i)=xp3(1,i)
1059 y2(i)=xp3(2,i)
1060 z2(i)=xp3(3,i)
1061 x3(i)=xp4(1,i)
1062 y3(i)=xp4(2,i)
1063 z3(i)=xp4(3,i)
1064 x12(i)=x2(i)-x1(i)
1065 y12(i)=y2(i)-y1(i)
1066 z12(i)=z2(i)-z1(i)
1067 x13(i)=x3(i)-x1(i)
1068 y13(i)=y3(i)-y1(i)
1069 z13(i)=z3(i)-z1(i)
1070 n1=y12(i)*z13(i)-z12(i)*y13(i)
1071 n2=z12(i)*x13(i)-x12(i)*z13(i)
1072 n3=x12(i)*y13(i)-y12(i)*x13(i)
1073 ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
1074 ntx(3,i)=ntn*n1
1075 nty(3,i)=ntn*n2
1076 ntz(3,i)=ntn*n3
1077 d(i) =-ntx(3,i)*x1(i)-nty(3,i)*y1(i)-ntz(3,i)*z1(i)
1078C-----
1079 x23(i)=x3(i)-x2(i)
1080 y23(i)=y3(i)-y2(i)
1081 z23(i)=z3(i)-z2(i)
1082 350 CONTINUE
1083C-------------------------------
1084C POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
1085C-------------------------------
1086 DO 375 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1087 il =ksc(nls)
1088 i =nls-ndeb
1089C-----
1090 eh =(abs(ntx(3,i)/(dgr*an))**expn)*an
1091 . +(abs(nty(3,i)/(dgr*bn))**expn)*bn
1092 . +(abs(ntz(3,i)/(dgr*cn))**expn)*cn
1093C X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
1094C LAMBDA1=EXP(LOG(1./EH)/EXPN)
1095 lambda1=(max(em20,eh))**(-one/expn)
1096 xh1 =abs(lambda1*ntx(3,i)/(dgr*an))**(one/(dgr-one))
1097 IF (ntx(3,i)<zero) xh1=-xh1
1098 yh1 =abs(lambda1*nty(3,i)/(dgr*bn))**(one/(dgr-one))
1099 IF (nty(3,i)<zero) yh1=-yh1
1100 zh1 =abs(lambda1*ntz(3,i)/(dgr*cn))**(one/(dgr-one))
1101 IF (ntz(3,i)<zero) zh1=-zh1
1102 mu1 =-ntx(3,i)*xh1-nty(3,i)*yh1-ntz(3,i)*zh1-d(i)
1103C LAMBDA2=-LAMBDA1
1104 xh2 =-xh1
1105 yh2 =-yh1
1106 zh2 =-zh1
1107C MU2 =-NTX(3,I)*XH2-NTY(3,I)*YH2-NTZ(3,I)*ZH2-D(I)
1108 mu2=-mu1-two*d(i)
1109 xtk(3,i) =xh1+mu1*ntx(3,i)
1110 ytk(3,i) =yh1+mu1*nty(3,i)
1111 ztk(3,i) =zh1+mu1*ntz(3,i)
1112 xtk2(i)=xh2+mu2*ntx(3,i)
1113 ytk2(i)=yh2+mu2*nty(3,i)
1114 ztk2(i)=zh2+mu2*ntz(3,i)
1115C-------------------------------
1116 375 CONTINUE
1117C-------------------------------
1118C RAMENER LE PT PK SUR LE TRIANGLE.
1119C-------------------------------
1120 DO 400 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1121 il =ksc(nls)
1122 i =nls-ndeb
1123C-----
1124 dx1=xtk(3,i)-x1(i)
1125 dy1=ytk(3,i)-y1(i)
1126 dz1=ztk(3,i)-z1(i)
1127 dx2=xtk(3,i)-x2(i)
1128 dy2=ytk(3,i)-y2(i)
1129 dz2=ztk(3,i)-z2(i)
1130 out = (dy1*dz2-dy2*dz1)*ntx(3,i)
1131 . +(dz1*dx2-dz2*dx1)*nty(3,i)
1132 . +(dx1*dy2-dy1*dx2)*ntz(3,i)
1133 IF (out<zero) THEN
1134C PROJECTION SUR 1,2
1135 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
1136 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
1137 bet=ps/max(em20,nr)
1138 bet=max(bet,zero)
1139 bet=min(bet,one)
1140 alp=1.-bet
1141 xtk(3,i)=alp*x1(i)+bet*x2(i)
1142 ytk(3,i)=alp*y1(i)+bet*y2(i)
1143 ztk(3,i)=alp*z1(i)+bet*z2(i)
1144 ENDIF
1145 dx3=xtk(3,i)-x3(i)
1146 dy3=ytk(3,i)-y3(i)
1147 dz3=ztk(3,i)-z3(i)
1148 out = (dy2*dz3-dy3*dz2)*ntx(3,i)
1149 . +(dz2*dx3-dz3*dx2)*nty(3,i)
1150 . +(dx2*dy3-dy2*dx3)*ntz(3,i)
1151 IF (out<zero) THEN
1152C PROJECTION SUR 2,3
1153 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
1154 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
1155 bet=ps/max(em20,nr)
1156 bet=max(bet,zero)
1157 bet=min(bet,one)
1158 alp=1.-bet
1159 xtk(3,i)=alp*x2(i)+bet*x3(i)
1160 ytk(3,i)=alp*y2(i)+bet*y3(i)
1161 ztk(3,i)=alp*z2(i)+bet*z3(i)
1162 ENDIF
1163 out = (dy3*dz1-dy1*dz3)*ntx(3,i)
1164 . +(dz3*dx1-dz1*dx3)*nty(3,i)
1165 . +(dx3*dy1-dy3*dx1)*ntz(3,i)
1166 IF (out<0.) THEN
1167C PROJECTION SUR 3,1
1168 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
1169 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
1170 bet=ps/max(em20,nr)
1171 bet=max(bet,zero)
1172 bet=min(bet,one)
1173 alp=1.-bet
1174 xtk(3,i)=alp*x3(i)+bet*x1(i)
1175 ytk(3,i)=alp*y3(i)+bet*y1(i)
1176 ztk(3,i)=alp*z3(i)+bet*z1(i)
1177 ENDIF
1178C-----
1179 dx1=xtk2(i)-x1(i)
1180 dy1=ytk2(i)-y1(i)
1181 dz1=ztk2(i)-z1(i)
1182 dx2=xtk2(i)-x2(i)
1183 dy2=ytk2(i)-y2(i)
1184 dz2=ztk2(i)-z2(i)
1185 out = (dy1*dz2-dy2*dz1)*ntx(3,i)
1186 . +(dz1*dx2-dz2*dx1)*nty(3,i)
1187 . +(dx1*dy2-dy1*dx2)*ntz(3,i)
1188 IF (out<zero) THEN
1189C PROJECTION SUR 1,2
1190 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
1191 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
1192 bet=ps/max(em20,nr)
1193 bet=max(bet,zero)
1194 bet=min(bet,one)
1195 alp=1.-bet
1196 xtk2(i)=alp*x1(i)+bet*x2(i)
1197 ytk2(i)=alp*y1(i)+bet*y2(i)
1198 ztk2(i)=alp*z1(i)+bet*z2(i)
1199 ENDIF
1200 dx3=xtk2(i)-x3(i)
1201 dy3=ytk2(i)-y3(i)
1202 dz3=ztk2(i)-z3(i)
1203 out = (dy2*dz3-dy3*dz2)*ntx(3,i)
1204 . +(dz2*dx3-dz3*dx2)*nty(3,i)
1205 . +(dx2*dy3-dy2*dx3)*ntz(3,i)
1206 IF (out<zero) THEN
1207C PROJECTION SUR 2,3
1208 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
1209 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
1210 bet=ps/max(em20,nr)
1211 bet=max(bet,zero)
1212 bet=min(bet,one)
1213 alp=one-bet
1214 xtk2(i)=alp*x2(i)+bet*x3(i)
1215 ytk2(i)=alp*y2(i)+bet*y3(i)
1216 ztk2(i)=alp*z2(i)+bet*z3(i)
1217 ENDIF
1218 out = (dy3*dz1-dy1*dz3)*ntx(3,i)
1219 . +(dz3*dx1-dz1*dx3)*nty(3,i)
1220 . +(dx3*dy1-dy3*dx1)*ntz(3,i)
1221 IF (out<zero) THEN
1222C PROJECTION SUR 3,1
1223 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
1224 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
1225 bet=ps/max(em20,nr)
1226 bet=max(bet,zero)
1227 bet=min(bet,one)
1228 alp=one-bet
1229 xtk2(i)=alp*x3(i)+bet*x1(i)
1230 ytk2(i)=alp*y3(i)+bet*y1(i)
1231 ztk2(i)=alp*z3(i)+bet*z1(i)
1232 ENDIF
1233C-------------------------------
1234 400 CONTINUE
1235C-------------------------------
1236C PROJECTION DE PK SUR L ET PENETRATION.
1237C-------------------------------
1238 DO 425 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1239 il =ksc(nls)
1240 i =nls-ndeb
1241 IF (iactiv(3,il)==-1) GOTO 425
1242C-----
1243 xh =xtk(3,i)
1244 yh =ytk(3,i)
1245 zh =ztk(3,i)
1246 xkn1 =xh**(idg-1)
1247 sgnxkn=-one
1248 IF (xkn1*xh>=zero) sgnxkn=one
1249 ykn1 =yh**(idg-1)
1250 sgnykn=-one
1251 IF (ykn1*yh>=zero) sgnykn=one
1252 zkn1 =zh**(idg-1)
1253 sgnzkn=-one
1254 IF (zkn1*zh>=zero) sgnzkn=one
1255 n11 =sgnxkn*xkn1*an
1256 n12 =sgnykn*ykn1*bn
1257 n13 =sgnzkn*zkn1*cn
1258 nr1 =n11*n11+n12*n12+n13*n13
1259 nr1 =sqrt(nr1)
1260 em =n11*xtk(3,i)+n12*ytk(3,i)+n13*ztk(3,i)
1261 IF (em<=one) THEN
1262 lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
1263 . / max(exp((dgr-one)*log(em20)/dgr),nr1)
1264 inside1=1
1265 ELSE
1266 inside1=0
1267 ENDIF
1268C-----
1269 xh =xtk2(i)
1270 yh =ytk2(i)
1271 zh =ztk2(i)
1272 xkn1 =xh**(idg-1)
1273 sgnxkn=-one
1274 IF (xkn1*xh>=zero) sgnxkn=one
1275 ykn1 =yh**(idg-1)
1276 sgnykn=-one
1277 IF (ykn1*yh>=zero) sgnykn=one
1278 zkn1 =zh**(idg-1)
1279 sgnzkn=-one
1280 IF (zkn1*zh>=zero) sgnzkn=one
1281 n21 =sgnxkn*xkn1*an
1282 n22 =sgnykn*ykn1*bn
1283 n23 =sgnzkn*zkn1*cn
1284 nr2 =n21*n21+n22*n22+n23*n23
1285 nr2 =sqrt(nr2)
1286 em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
1287 IF (em<=one) THEN
1288 lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
1289 . / max(exp((dgr-one)*log(em20)/dgr),nr2)
1290 inside2=1
1291 ELSE
1292 inside2=0
1293 ENDIF
1294C-----
1295 IF (inside1==0.AND.inside2==0) THEN
1296 iactiv(3,il)=0
1297 ELSE
1298C-----
1299 IF (iactiv(3,il)==0) THEN
1300 IF (inside1/=0.AND.inside2/=0) THEN
1301 IF (abs(lambda1)>=abs(lambda2)) THEN
1302 xm(i)=xtk(3,i)-lambda1*n11/max(em20,nr1)
1303 ym(i)=ytk(3,i)-lambda1*n12/max(em20,nr1)
1304 zm(i)=ztk(3,i)-lambda1*n13/max(em20,nr1)
1305 ELSE
1306 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
1307 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
1308 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
1309 xtk(3,i)=xtk2(i)
1310 ytk(3,i)=ytk2(i)
1311 ztk(3,i)=ztk2(i)
1312 ENDIF
1313 ELSEIF(inside1/=0) THEN
1314 xm(i)=xtk(3,i)-lambda1*n11/max(em20,nr1)
1315 ym(i)=ytk(3,i)-lambda1*n12/max(em20,nr1)
1316 zm(i)=ztk(3,i)-lambda1*n13/max(em20,nr1)
1317 ELSEIF(inside2/=0) THEN
1318 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
1319 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
1320 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
1321 xtk(3,i)=xtk2(i)
1322 ytk(3,i)=ytk2(i)
1323 ztk(3,i)=ztk2(i)
1324 ENDIF
1325C-----
1326 ELSE
1327 xh=hold(1,4*(il-1)+3)
1328 yh=hold(2,4*(il-1)+3)
1329 zh=hold(3,4*(il-1)+3)
1330 n1=nold(1,4*(il-1)+3)
1331 n2=nold(2,4*(il-1)+3)
1332 n3=nold(3,4*(il-1)+3)
1333 lambda1=(xh-xtk(3,i))*n1
1334 . +(yh-ytk(3,i))*n2
1335 . +(zh-ztk(3,i))*n3
1336 lambda2=(xh-xtk2(i))*n1
1337 . +(yh-ytk2(i))*n2
1338 . +(zh-ztk2(i))*n3
1339 IF (inside1/=0.AND.inside2/=0) THEN
1340 IF (abs(lambda1)>=abs(lambda2)) THEN
1341 xm(i)=xtk(3,i)+lambda1*n1
1342 ym(i)=ytk(3,i)+lambda1*n2
1343 zm(i)=ztk(3,i)+lambda1*n3
1344 ELSE
1345 xm(i)=xtk2(i)+lambda2*n1
1346 ym(i)=ytk2(i)+lambda2*n2
1347 zm(i)=ztk2(i)+lambda2*n3
1348 xtk(3,i)=xtk2(i)
1349 ytk(3,i)=ytk2(i)
1350 ztk(3,i)=ztk2(i)
1351 ENDIF
1352 ELSEIF(inside1/=0) THEN
1353 xm(i)=xtk(3,i)+lambda1*n1
1354 ym(i)=ytk(3,i)+lambda1*n2
1355 zm(i)=ztk(3,i)+lambda1*n3
1356 ELSEIF(inside2/=0) THEN
1357 xm(i)=xtk2(i)+lambda2*n1
1358 ym(i)=ytk2(i)+lambda2*n2
1359 zm(i)=ztk2(i)+lambda2*n3
1360 xtk(3,i)=xtk2(i)
1361 ytk(3,i)=ytk2(i)
1362 ztk(3,i)=ztk2(i)
1363 ENDIF
1364C one more iteration.
1365 xkn1 =xm(i)**(idg-1)
1366 sgnxkn=-one
1367 IF (xkn1*xm(i)>=zero) sgnxkn=one
1368 ykn1 =ym(i)**(idg-1)
1369 sgnykn=-one
1370 IF (ykn1*ym(i)>=zero) sgnykn=one
1371 zkn1 =zm(i)**(idg-1)
1372 sgnzkn=-one
1373 IF (zkn1*zm(i)>=zero) sgnzkn=one
1374 n1 =sgnxkn*xkn1*an
1375 n2 =sgnykn*ykn1*bn
1376 n3 =sgnzkn*zkn1*cn
1377 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
1378 xm(i)=xm(i)/max(em20,em**(one/dgr))
1379 ym(i)=ym(i)/max(em20,em**(one/dgr))
1380 zm(i)=zm(i)/max(em20,em**(one/dgr))
1381 xkn1 =xm(i)**(idg-1)
1382 sgnxkn=-one
1383 IF (xkn1*xm(i)>=zero) sgnxkn=one
1384 ykn1 =ym(i)**(idg-1)
1385 sgnykn=-one
1386 IF (ykn1*ym(i)>=zero) sgnykn=one
1387 zkn1 =zm(i)**(idg-1)
1388 sgnzkn=-one
1389 IF (zkn1*zm(i)>=zero) sgnzkn=one
1390 n1 =sgnxkn*xkn1*an
1391 n2 =sgnykn*ykn1*bn
1392 n3 =sgnzkn*zkn1*cn
1393 nr =n1*n1+n2*n2+n3*n3
1394 nr =one/max(em20,sqrt(nr))
1395 n1 =n1*nr
1396 n2 =n2*nr
1397 n3 =n3*nr
1398 lambda1=(xm(i)-xtk(3,i))*n1
1399 . +(ym(i)-ytk(3,i))*n2
1400 . +(zm(i)-ztk(3,i))*n3
1401 xm(i)=xtk(3,i)+lambda1*n1
1402 ym(i)=ytk(3,i)+lambda1*n2
1403 zm(i)=ztk(3,i)+lambda1*n3
1404 ENDIF
1405 iactiv(3,il)=iactiv(3,il)+1
1406 ENDIF
1407 425 CONTINUE
1408C-----
1409#include "vectorize.inc"
1410 DO 450 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1411 il =ksc(nls)
1412 i =nls-ndeb
1413C-----
1414 IF (iactiv(3,il)<=0) GOTO 450
1415C-----
1416C projection radiale de M.
1417 xkn1 =xm(i)**(idg-1)
1418 sgnxkn=-one
1419 IF (xkn1*xm(i)>=zero) sgnxkn=one
1420 ykn1 =ym(i)**(idg-1)
1421 sgnykn=-one
1422 IF (ykn1*ym(i)>=zero) sgnykn=one
1423 zkn1 =zm(i)**(idg-1)
1424 sgnzkn=-one
1425 IF (zkn1*zm(i)>=zero) sgnzkn=one
1426 n1 =sgnxkn*xkn1*an
1427 n2 =sgnykn*ykn1*bn
1428 n3 =sgnzkn*zkn1*cn
1429 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
1430C------
1431 xm(i)=xm(i)/max(em20,em**(one/dgr))
1432 ym(i)=ym(i)/max(em20,em**(one/dgr))
1433 zm(i)=zm(i)/max(em20,em**(one/dgr))
1434C------
1435C normale a l'ellipsoide en M.
1436 xkn1 =xm(i)**(idg-1)
1437 sgnxkn=-one
1438 IF (xkn1*xm(i)>=zero) sgnxkn=one
1439 ykn1 =ym(i)**(idg-1)
1440 sgnykn=-one
1441 IF (ykn1*ym(i)>=zero) sgnykn=one
1442 zkn1 =zm(i)**(idg-1)
1443 sgnzkn=-one
1444 IF (zkn1*zm(i)>=zero) sgnzkn=one
1445 n1 =sgnxkn*xkn1*an
1446 n2 =sgnykn*ykn1*bn
1447 n3 =sgnzkn*zkn1*cn
1448 nr =n1*n1+n2*n2+n3*n3
1449 nr =one/max(em20,sqrt(nr))
1450 nxi(3,i)=n1*nr
1451 nyi(3,i)=n2*nr
1452 nzi(3,i)=n3*nr
1453 xi(3,i)=xm(i)
1454 yi(3,i)=ym(i)
1455 zi(3,i)=zm(i)
1456 depth(3,i)=xm(i)*nxi(3,i)+ym(i)*nyi(3,i)+zm(i)*nzi(3,i)
1457C------
1458 penet(3,i)=(xtk(3,i)-xm(i))*nxi(3,i)
1459 . +(ytk(3,i)-ym(i))*nyi(3,i)
1460 . +(ztk(3,i)-zm(i))*nzi(3,i)
1461 penet(3,i)=-penet(3,i)
1462 penet(3,i)=max(zero,penet(3,i))
1463 IF (depth(3,i)-penet(3,i)<em10*depth(3,i)) THEN
1464 penet(3,i)=zero
1465 iactiv(3,il)=-2
1466 ENDIF
1467 ansmx=max(ansmx,penet(3,i))
1468 450 CONTINUE
1469C-------------------------------
1470C MESSAGES DESACTIVATION
1471C-------------------------------
1472 DO 460 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1473 il =ksc(nls)
1474 i =nls-ndeb
1475C-----
1476 IF (iactiv(3,il)==-2) THEN
1477 in1=itab(ksi(1,il))
1478 in2=itab(ksi(2,il))
1479 in3=itab(ksi(3,il))
1480 in4=itab(ksi(4,il))
1481#include "lockon.inc"
1482 WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
1483 WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
1484 . ' FROM 4NODES ELEMENT/SEGMENT,',
1485 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
1486 . in3,in4,'ISOCENTER(',in1,in2,in3,in4,')'
1487 WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
1488 WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
1489 . ' FROM 4NODES ELEMENT/SEGMENT,',
1490 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
1491 . in3,in4,'ISOCENTER(',in1,in2,in3,in4,')'
1492#include "lockoff.inc"
1493 iactiv(3,il)=-1
1494 ENDIF
1495 460 CONTINUE
1496C-------------------------------
1497C 4EME TRIANGLE.
1498C-------------------------------
1499 DO 475 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1500 il =ksc(nls)
1501 i =nls-ndeb
1502C-----
1503 x1(i)=gx(1,i)
1504 y1(i)=gx(2,i)
1505 z1(i)=gx(3,i)
1506 x2(i)=xp4(1,i)
1507 y2(i)=xp4(2,i)
1508 z2(i)=xp4(3,i)
1509 x3(i)=xp1(1,i)
1510 y3(i)=xp1(2,i)
1511 z3(i)=xp1(3,i)
1512 x12(i)=x2(i)-x1(i)
1513 y12(i)=y2(i)-y1(i)
1514 z12(i)=z2(i)-z1(i)
1515 x13(i)=x3(i)-x1(i)
1516 y13(i)=y3(i)-y1(i)
1517 z13(i)=z3(i)-z1(i)
1518 n1=y12(i)*z13(i)-z12(i)*y13(i)
1519 n2=z12(i)*x13(i)-x12(i)*z13(i)
1520 n3=x12(i)*y13(i)-y12(i)*x13(i)
1521 ntn=one/max(em20,sqrt(n1*n1+n2*n2+n3*n3))
1522 ntx(4,i)=ntn*n1
1523 nty(4,i)=ntn*n2
1524 ntz(4,i)=ntn*n3
1525 d(i) =-ntx(4,i)*x1(i)-nty(4,i)*y1(i)-ntz(4,i)*z1(i)
1526C-----
1527 x23(i)=x3(i)-x2(i)
1528 y23(i)=y3(i)-y2(i)
1529 z23(i)=z3(i)-z2(i)
1530 475 CONTINUE
1531C-------------------------------
1532C POINTS K,H SUR E,L / LA DISTANCE D(K,H) EST LOCALEMENT MAXIMUM
1533C-------------------------------
1534 DO 500 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1535 il =ksc(nls)
1536 i =nls-ndeb
1537C-----
1538 eh =(abs(ntx(4,i)/(dgr*an))**expn)*an
1539 . +(abs(nty(4,i)/(dgr*bn))**expn)*bn
1540 . +(abs(ntz(4,i)/(dgr*cn))**expn)*cn
1541C X EST DU SIGNE DE LAMBDA*NTX, IDEM EN Y ET Z.
1542C LAMBDA1=EXP(LOG(1./EH)/EXPN)
1543 lambda1=(max(em20,eh))**(-one/expn)
1544 xh1 =abs(lambda1*ntx(4,i)/(dgr*an))**(one/(dgr-one))
1545 IF (ntx(4,i)<zero) xh1=-xh1
1546 yh1 =abs(lambda1*nty(4,i)/(dgr*bn))**(one/(dgr-one))
1547 IF (nty(4,i)<zero) yh1=-yh1
1548 zh1 =abs(lambda1*ntz(4,i)/(dgr*cn))**(one/(dgr-one))
1549 IF (ntz(4,i)<zero) zh1=-zh1
1550 mu1 =-ntx(4,i)*xh1-nty(4,i)*yh1-ntz(4,i)*zh1-d(i)
1551C LAMBDA2=-LAMBDA1
1552 xh2 =-xh1
1553 yh2 =-yh1
1554 zh2 =-zh1
1555C MU2 =-NTX(4,I)*XH2-NTY(4,I)*YH2-NTZ(4,I)*ZH2-D(I)
1556 mu2=-mu1-two*d(i)
1557 xtk(4,i) =xh1+mu1*ntx(4,i)
1558 ytk(4,i) =yh1+mu1*nty(4,i)
1559 ztk(4,i) =zh1+mu1*ntz(4,i)
1560 xtk2(i)=xh2+mu2*ntx(4,i)
1561 ytk2(i)=yh2+mu2*nty(4,i)
1562 ztk2(i)=zh2+mu2*ntz(4,i)
1563C-------------------------------
1564 500 CONTINUE
1565C-------------------------------
1566C RAMENER LE PT PK SUR LE TRIANGLE.
1567C-------------------------------
1568 DO 525 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1569 il =ksc(nls)
1570 i =nls-ndeb
1571C-----
1572 dx1=xtk(4,i)-x1(i)
1573 dy1=ytk(4,i)-y1(i)
1574 dz1=ztk(4,i)-z1(i)
1575 dx2=xtk(4,i)-x2(i)
1576 dy2=ytk(4,i)-y2(i)
1577 dz2=ztk(4,i)-z2(i)
1578 out = (dy1*dz2-dy2*dz1)*ntx(4,i)
1579 . +(dz1*dx2-dz2*dx1)*nty(4,i)
1580 . +(dx1*dy2-dy1*dx2)*ntz(4,i)
1581 IF (out<zero) THEN
1582C PROJECTION SUR 1,2
1583 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
1584 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
1585 bet=ps/max(em20,nr)
1586 bet=max(bet,zero)
1587 bet=min(bet,one)
1588 alp=1.-bet
1589 xtk(4,i)=alp*x1(i)+bet*x2(i)
1590 ytk(4,i)=alp*y1(i)+bet*y2(i)
1591 ztk(4,i)=alp*z1(i)+bet*z2(i)
1592 ENDIF
1593 dx3=xtk(4,i)-x3(i)
1594 dy3=ytk(4,i)-y3(i)
1595 dz3=ztk(4,i)-z3(i)
1596 out = (dy2*dz3-dy3*dz2)*ntx(4,i)
1597 . +(dz2*dx3-dz3*dx2)*nty(4,i)
1598 . +(dx2*dy3-dy2*dx3)*ntz(4,i)
1599 IF (out<0.) THEN
1600C PROJECTION SUR 2,3
1601 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
1602 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
1603 bet=ps/max(em20,nr)
1604 bet=max(bet,zero)
1605 bet=min(bet,one)
1606 alp=1.-bet
1607 xtk(4,i)=alp*x2(i)+bet*x3(i)
1608 ytk(4,i)=alp*y2(i)+bet*y3(i)
1609 ztk(4,i)=alp*z2(i)+bet*z3(i)
1610 ENDIF
1611 out = (dy3*dz1-dy1*dz3)*ntx(4,i)
1612 . +(dz3*dx1-dz1*dx3)*nty(4,i)
1613 . +(dx3*dy1-dy3*dx1)*ntz(4,i)
1614 IF (out<0.) THEN
1615C PROJECTION SUR 3,1
1616 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
1617 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
1618 bet=ps/max(em20,nr)
1619 bet=max(bet,zero)
1620 bet=min(bet,one)
1621 alp=1.-bet
1622 xtk(4,i)=alp*x3(i)+bet*x1(i)
1623 ytk(4,i)=alp*y3(i)+bet*y1(i)
1624 ztk(4,i)=alp*z3(i)+bet*z1(i)
1625 ENDIF
1626C-----
1627 dx1=xtk2(i)-x1(i)
1628 dy1=ytk2(i)-y1(i)
1629 dz1=ztk2(i)-z1(i)
1630 dx2=xtk2(i)-x2(i)
1631 dy2=ytk2(i)-y2(i)
1632 dz2=ztk2(i)-z2(i)
1633 out = (dy1*dz2-dy2*dz1)*ntx(4,i)
1634 . +(dz1*dx2-dz2*dx1)*nty(4,i)
1635 . +(dx1*dy2-dy1*dx2)*ntz(4,i)
1636 IF (out<zero) THEN
1637C PROJECTION SUR 1,2
1638 ps =dx1*x12(i)+dy1*y12(i)+dz1*z12(i)
1639 nr =x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i)
1640 bet=ps/max(em20,nr)
1641 bet=max(bet,zero)
1642 bet=min(bet,one)
1643 alp=1.-bet
1644 xtk2(i)=alp*x1(i)+bet*x2(i)
1645 ytk2(i)=alp*y1(i)+bet*y2(i)
1646 ztk2(i)=alp*z1(i)+bet*z2(i)
1647 ENDIF
1648 dx3=xtk2(i)-x3(i)
1649 dy3=ytk2(i)-y3(i)
1650 dz3=ztk2(i)-z3(i)
1651 out = (dy2*dz3-dy3*dz2)*ntx(4,i)
1652 . +(dz2*dx3-dz3*dx2)*nty(4,i)
1653 . +(dx2*dy3-dy2*dx3)*ntz(4,i)
1654 IF (out<zero) THEN
1655C PROJECTION SUR 2,3
1656 ps =dx2*x23(i)+dy2*y23(i)+dz2*z23(i)
1657 nr =x23(i)*x23(i)+y23(i)*y23(i)+z23(i)*z23(i)
1658 bet=ps/max(em20,nr)
1659 bet=max(bet,zero)
1660 bet=min(bet,one)
1661 alp=one-bet
1662 xtk2(i)=alp*x2(i)+bet*x3(i)
1663 ytk2(i)=alp*y2(i)+bet*y3(i)
1664 ztk2(i)=alp*z2(i)+bet*z3(i)
1665 ENDIF
1666 out = (dy3*dz1-dy1*dz3)*ntx(4,i)
1667 . +(dz3*dx1-dz1*dx3)*nty(4,i)
1668 . +(dx3*dy1-dy3*dx1)*ntz(4,i)
1669 IF (out<zero) THEN
1670C PROJECTION SUR 3,1
1671 ps =-dx3*x13(i)-dy3*y13(i)-dz3*z13(i)
1672 nr =x13(i)*x13(i)+y13(i)*y13(i)+z13(i)*z13(i)
1673 bet=ps/max(em20,nr)
1674 bet=max(bet,zero)
1675 bet=min(bet,one)
1676 alp=one-bet
1677 xtk2(i)=alp*x3(i)+bet*x1(i)
1678 ytk2(i)=alp*y3(i)+bet*y1(i)
1679 ztk2(i)=alp*z3(i)+bet*z1(i)
1680 ENDIF
1681C-------------------------------
1682 525 CONTINUE
1683C-------------------------------
1684C PROJECTION DE PK SUR L ET PENETRATION.
1685C-------------------------------
1686 DO 550 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1687 il =ksc(nls)
1688 i =nls-ndeb
1689 IF (iactiv(4,il)==-1) GOTO 550
1690C-----
1691 xh =xtk(4,i)
1692 yh =ytk(4,i)
1693 zh =ztk(4,i)
1694 xkn1 =xh**(idg-1)
1695 sgnxkn=-one
1696 IF (xkn1*xh>=zero) sgnxkn=one
1697 ykn1 =yh**(idg-1)
1698 sgnykn=-one
1699 IF (ykn1*yh>=zero) sgnykn=one
1700 zkn1 =zh**(idg-1)
1701 sgnzkn=-one
1702 IF (zkn1*zh>=zero) sgnzkn=one
1703 n11 =sgnxkn*xkn1*an
1704 n12 =sgnykn*ykn1*bn
1705 n13 =sgnzkn*zkn1*cn
1706 nr1 =n11*n11+n12*n12+n13*n13
1707 nr1 =sqrt(nr1)
1708 em =n11*xtk(4,i)+n12*ytk(4,i)+n13*ztk(4,i)
1709 IF (em<=one) THEN
1710 lambda1=(em-exp((dgr-one)*log(max(em,em20))/dgr))
1711 . / max(exp((dgr-one)*log(em20)/dgr),nr1)
1712 inside1=1
1713 ELSE
1714 inside1=0
1715 ENDIF
1716C-----
1717 xh =xtk2(i)
1718 yh =ytk2(i)
1719 zh =ztk2(i)
1720 xkn1 =xh**(idg-1)
1721 sgnxkn=-one
1722 IF (xkn1*xh>=zero) sgnxkn=one
1723 ykn1 =yh**(idg-1)
1724 sgnykn=-one
1725 IF (ykn1*yh>=one) sgnykn=one
1726 zkn1 =zh**(idg-1)
1727 sgnzkn=-one
1728 IF (zkn1*zh>=zero) sgnzkn=one
1729 n21 =sgnxkn*xkn1*an
1730 n22 =sgnykn*ykn1*bn
1731 n23 =sgnzkn*zkn1*cn
1732 nr2 =n21*n21+n22*n22+n23*n23
1733 nr2 =sqrt(nr2)
1734 em =n21*xtk2(i)+n22*ytk2(i)+n23*ztk2(i)
1735 IF (em<=one) THEN
1736 lambda2=(em-exp((dgr-one)*log(max(em,em20))/dgr))
1737 . / max(exp((dgr-one)*log(em20)/dgr),nr2)
1738 inside2=1
1739 ELSE
1740 inside2=0
1741 ENDIF
1742C-----
1743 IF (inside1==0.AND.inside2==0) THEN
1744 iactiv(4,il)=0
1745 ELSE
1746C-----
1747 IF (iactiv(4,il)==0) THEN
1748 IF (inside1/=0.AND.inside2/=0) THEN
1749 IF (abs(lambda1)>=abs(lambda2)) THEN
1750 xm(i)=xtk(4,i)-lambda1*n11/max(em20,nr1)
1751 ym(i)=ytk(4,i)-lambda1*n12/max(em20,nr1)
1752 zm(i)=ztk(4,i)-lambda1*n13/max(em20,nr1)
1753 ELSE
1754 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
1755 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
1756 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
1757 xtk(4,i)=xtk2(i)
1758 ytk(4,i)=ytk2(i)
1759 ztk(4,i)=ztk2(i)
1760 ENDIF
1761 ELSEIF(inside1/=0) THEN
1762 xm(i)=xtk(4,i)-lambda1*n11/max(em20,nr1)
1763 ym(i)=ytk(4,i)-lambda1*n12/max(em20,nr1)
1764 zm(i)=ztk(4,i)-lambda1*n13/max(em20,nr1)
1765 ELSEIF(inside2/=0) THEN
1766 xm(i)=xtk2(i)-lambda2*n21/max(em20,nr2)
1767 ym(i)=ytk2(i)-lambda2*n22/max(em20,nr2)
1768 zm(i)=ztk2(i)-lambda2*n23/max(em20,nr2)
1769 xtk(4,i)=xtk2(i)
1770 ytk(4,i)=ytk2(i)
1771 ztk(4,i)=ztk2(i)
1772 ENDIF
1773C-----
1774 ELSE
1775 xh=hold(1,4*(il-1)+4)
1776 yh=hold(2,4*(il-1)+4)
1777 zh=hold(3,4*(il-1)+4)
1778 n1=nold(1,4*(il-1)+4)
1779 n2=nold(2,4*(il-1)+4)
1780 n3=nold(3,4*(il-1)+4)
1781 lambda1=(xh-xtk(4,i))*n1
1782 . +(yh-ytk(4,i))*n2
1783 . +(zh-ztk(4,i))*n3
1784 lambda2=(xh-xtk2(i))*n1
1785 . +(yh-ytk2(i))*n2
1786 . +(zh-ztk2(i))*n3
1787 IF (inside1/=0.AND.inside2/=0) THEN
1788 IF (abs(lambda1)>=abs(lambda2)) THEN
1789 xm(i)=xtk(4,i)+lambda1*n1
1790 ym(i)=ytk(4,i)+lambda1*n2
1791 zm(i)=ztk(4,i)+lambda1*n3
1792 ELSE
1793 xm(i)=xtk2(i)+lambda2*n1
1794 ym(i)=ytk2(i)+lambda2*n2
1795 zm(i)=ztk2(i)+lambda2*n3
1796 xtk(4,i)=xtk2(i)
1797 ytk(4,i)=ytk2(i)
1798 ztk(4,i)=ztk2(i)
1799 ENDIF
1800 ELSEIF(inside1/=0) THEN
1801 xm(i)=xtk(4,i)+lambda1*n1
1802 ym(i)=ytk(4,i)+lambda1*n2
1803 zm(i)=ztk(4,i)+lambda1*n3
1804 ELSEIF(inside2/=0) THEN
1805 xm(i)=xtk2(i)+lambda2*n1
1806 ym(i)=ytk2(i)+lambda2*n2
1807 zm(i)=ztk2(i)+lambda2*n3
1808 xtk(4,i)=xtk2(i)
1809 ytk(4,i)=ytk2(i)
1810 ztk(4,i)=ztk2(i)
1811 ENDIF
1812C one more iteration.
1813 xkn1 =xm(i)**(idg-1)
1814 sgnxkn=-one
1815 IF (xkn1*xm(i)>=zero) sgnxkn=one
1816 ykn1 =ym(i)**(idg-1)
1817 sgnykn=-one
1818 IF (ykn1*ym(i)>=zero) sgnykn=one
1819 zkn1 =zm(i)**(idg-1)
1820 sgnzkn=-one
1821 IF (zkn1*zm(i)>=zero) sgnzkn=one
1822 n1 =sgnxkn*xkn1*an
1823 n2 =sgnykn*ykn1*bn
1824 n3 =sgnzkn*zkn1*cn
1825 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
1826 xm(i)=xm(i)/max(em20,em**(one/dgr))
1827 ym(i)=ym(i)/max(em20,em**(one/dgr))
1828 zm(i)=zm(i)/max(em20,em**(one/dgr))
1829 xkn1 =xm(i)**(idg-1)
1830 sgnxkn=-one
1831 IF (xkn1*xm(i)>=zero) sgnxkn=one
1832 ykn1 =ym(i)**(idg-1)
1833 sgnykn=-one
1834 IF (ykn1*ym(i)>=zero) sgnykn=one
1835 zkn1 =zm(i)**(idg-1)
1836 sgnzkn=-one
1837 IF (zkn1*zm(i)>=zero) sgnzkn=one
1838 n1 =sgnxkn*xkn1*an
1839 n2 =sgnykn*ykn1*bn
1840 n3 =sgnzkn*zkn1*cn
1841 nr =n1*n1+n2*n2+n3*n3
1842 nr =1./max(em20,sqrt(nr))
1843 n1 =n1*nr
1844 n2 =n2*nr
1845 n3 =n3*nr
1846 lambda1=(xm(i)-xtk(4,i))*n1
1847 . +(ym(i)-ytk(4,i))*n2
1848 . +(zm(i)-ztk(4,i))*n3
1849 xm(i)=xtk(4,i)+lambda1*n1
1850 ym(i)=ytk(4,i)+lambda1*n2
1851 zm(i)=ztk(4,i)+lambda1*n3
1852 ENDIF
1853 iactiv(4,il)=iactiv(4,il)+1
1854 ENDIF
1855 550 CONTINUE
1856C-----
1857#include "vectorize.inc"
1858 DO 575 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1859 il =ksc(nls)
1860 i =nls-ndeb
1861C-----
1862 IF (iactiv(4,il)<=0) GOTO 575
1863C-----
1864C projection radiale de M.
1865 xkn1 =xm(i)**(idg-1)
1866 sgnxkn=-one
1867 IF (xkn1*xm(i)>=zero) sgnxkn=one
1868 ykn1 =ym(i)**(idg-1)
1869 sgnykn=-one
1870 IF (ykn1*ym(i)>=zero) sgnykn=one
1871 zkn1 =zm(i)**(idg-1)
1872 sgnzkn=-one
1873 IF (zkn1*zm(i)>=zero) sgnzkn=one
1874 n1 =sgnxkn*xkn1*an
1875 n2 =sgnykn*ykn1*bn
1876 n3 =sgnzkn*zkn1*cn
1877 em=n1*xm(i)+n2*ym(i)+n3*zm(i)
1878C------
1879 xm(i)=xm(i)/max(em20,em**(one/dgr))
1880 ym(i)=ym(i)/max(em20,em**(one/dgr))
1881 zm(i)=zm(i)/max(em20,em**(one/dgr))
1882C------
1883C normale a l'ellipsoide en M.
1884 xkn1 =xm(i)**(idg-1)
1885 sgnxkn=-one
1886 IF (xkn1*xm(i)>=zero) sgnxkn=one
1887 ykn1 =ym(i)**(idg-1)
1888 sgnykn=-one
1889 IF (ykn1*ym(i)>=zero) sgnykn=one
1890 zkn1 =zm(i)**(idg-1)
1891 sgnzkn=-one
1892 IF (zkn1*zm(i)>=zero) sgnzkn=one
1893 n1 =sgnxkn*xkn1*an
1894 n2 =sgnykn*ykn1*bn
1895 n3 =sgnzkn*zkn1*cn
1896 nr =n1*n1+n2*n2+n3*n3
1897 nr =one/max(em20,sqrt(nr))
1898 nxi(4,i)=n1*nr
1899 nyi(4,i)=n2*nr
1900 nzi(4,i)=n3*nr
1901 xi(4,i)=xm(i)
1902 yi(4,i)=ym(i)
1903 zi(4,i)=zm(i)
1904 depth(4,i)=xm(i)*nxi(4,i)+ym(i)*nyi(4,i)+zm(i)*nzi(4,i)
1905C------
1906 penet(4,i)=(xtk(4,i)-xm(i))*nxi(4,i)
1907 . +(ytk(4,i)-ym(i))*nyi(4,i)
1908 . +(ztk(4,i)-zm(i))*nzi(4,i)
1909 penet(4,i)=-penet(4,i)
1910 penet(4,i)=max(zero,penet(4,i))
1911 IF (depth(4,i)-penet(4,i)<em10*depth(4,i)) THEN
1912 penet(4,i)=zero
1913 iactiv(4,il)=-2
1914 ENDIF
1915 ansmx=max(ansmx,penet(4,i))
1916 575 CONTINUE
1917C-------------------------------
1918C MESSAGES DESACTIVATION
1919C-------------------------------
1920 DO 585 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1921 il =ksc(nls)
1922 i =nls-ndeb
1923C-----
1924 IF (iactiv(4,il)==-2) THEN
1925 in1=itab(ksi(1,il))
1926 in2=itab(ksi(2,il))
1927 in3=itab(ksi(3,il))
1928 in4=itab(ksi(4,il))
1929#include "lockon.inc"
1930 WRITE(istdo,'(A,I8)')' WARNING INTERFACE ',noint
1931 WRITE(istdo,'(A,A,A,/,2I8,A,4I8,A)')' 3NODES FACET',
1932 . ' FROM 4NODES ELEMENT/SEGMENT,',
1933 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
1934 . in4,in1,'ISOCENTER(',in1,in2,in3,in4,')'
1935 WRITE(iout ,'(A,I8)')' WARNING INTERFACE ',noint
1936 WRITE(iout,'(A,A,A,/,2I8,A,4I8,A)') ' 3NODES FACET',
1937 . ' FROM 4NODES ELEMENT/SEGMENT,',
1938 . ' DE-ACTIVATED FROM INTERFACE; FACET VERTICES:',
1939 . in4,in1,'ISOCENTER(',in1,in2,in3,in4,')'
1940#include "lockoff.inc"
1941 iactiv(4,il)=-1
1942 ENDIF
1943 585 CONTINUE
1944C------------------------------------------------------------
1945 RETURN
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
groupdef_mod
Definition groupdef_mod.F:662