OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i15for1.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i15for1 ../engine/source/interfaces/int15/i15for1.F
25!||--- called by ------------------------------------------------------
26!|| i15cmp ../engine/source/interfaces/int15/i15cmp.F
27!||--- uses -----------------------------------------------------
28!|| groupdef_mod ../common_source/modules/groupdef_mod.F
29!||====================================================================
30 SUBROUTINE i15for1(NDEB, NSC, STFAC,X ,V ,
31 2 KSURF ,IGRSURF ,BUFSF ,KSC ,KSI ,
32 3 IACTIV ,IOLD ,HOLD ,NOLD ,DOLD ,
33 4 XP1 ,XP2 ,XP3 ,XP4 ,GX ,
34 5 XTK ,YTK ,ZTK ,NTX ,NTY ,
35 6 NTZ ,PENET ,DEPTH ,XI ,YI ,
36 7 ZI ,NXI ,NYI ,NZI ,MS ,
37 8 DE ,NPC ,PLD ,WNF ,WTF ,
38 9 WNS ,FNORMX,FNORMY,FNORMZ,FTANGX,
39 A FTANGY,FTANGZ,DT2T ,NOINT ,NELTST ,
40 B ITYPTST ,VFRIC )
41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE groupdef_mod
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.inc"
49#include "comlock.inc"
50C-----------------------------------------------
51C G l o b a l P a r a m e t e r s
52C-----------------------------------------------
53#include "mvsiz_p.inc"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "com04_c.inc"
58#include "com08_c.inc"
59#include "units_c.inc"
60C-----------------------------------------------
61C D u m m y A r g u m e n t s
62C-----------------------------------------------
63 INTEGER NDEB, NSC
64 INTEGER KSURF,
65 . KSI(4,*),IACTIV(4,*),NPC(*),KSC(*),
66 . NOINT,NELTST,ITYPTST
67C REAL
68 my_real STFAC, BUFSF(*),
69 . X(3,*) , IOLD(3,*), HOLD(3,*), NOLD(3,*) ,DOLD(3,*),
70 . MS(*) , V(3,*), DE, PLD(*), XTK(4,*),
71 . YTK(4,*) ,ZTK(4,*) ,NTX(4,*) ,NTY(4,*) ,NTZ(4,*) ,
72 . PENET(4,*),DEPTH(4,*),
73 . xi(4,*) ,yi(4,*) ,zi(4,*) ,nxi(4,*) ,
74 . nyi(4,*) ,nzi(4,*) ,wnf(3,*) ,wtf(3,*) , wns(*),
75 . xp1(3,*) ,xp2(3,*) ,xp3(3,*) ,xp4(3,*) ,gx(3,*),
76 . fnormx,fnormy,fnormz,ftangx,ftangy,ftangz,
77 . dt2t,vfric
78 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
79C-----------------------------------------------
80C L o c a l V a r i a b l e s
81C-----------------------------------------------
82 INTEGER ADRBUF, I, IL, IN1, IN2, IN3, IN4, NLS
83 my_real
84 . ROT(9), X1, X2, X3, XM, YM, ZM,
85 . V1, V2, V3, VXM, VYM, VZM, VRX, VRY, VRZ,
86 . V1X2, V2X1, V1X3, V3X1, V2X3, V3X2,
87 . DELTX, DELTY, DELTZ, ND, PN, SCALE, NV,
88 . vxk, vyk, vzk, dvx, dvy, dvz, xh, yh, zh,
89 . dx1, dy1, dz1, dx2, dy2, dz2, dx3, dy3, dz3,
90 . dx, dy, dz,
91 . s1, s2, s3, s,
92 . f1, f2, f3, f4, f1g, f11, f12, f2g, f22, f23, f3g, f33, f34,
93 . f4g, f44, f41, fg,
94 . st1, st2, st3, st4, st1g, st11, st12, st2g, st22, st23, st3g,
95 . st33, st34, st4g, st44, st41, stg,
96 . ftx, fty, ftz, fnx, fny, fnz, ftn, fmax,fn,
97 . h, dti
98 my_real
99 . fxn(4,mvsiz),fyn(4,mvsiz),fzn(4,mvsiz),stf(4,mvsiz),
100 . fxt(4,mvsiz),fyt(4,mvsiz),fzt(4,mvsiz),
101 . l1(mvsiz),l2(mvsiz),l3(mvsiz),
102 . vxg(mvsiz),vyg(mvsiz),vzg(mvsiz),
103 . vx1(mvsiz),vy1(mvsiz),vz1(mvsiz),
104 . vx2(mvsiz),vy2(mvsiz),vz2(mvsiz),
105 . vx3(mvsiz),vy3(mvsiz),vz3(mvsiz),
106 . vx4(mvsiz),vy4(mvsiz),vz4(mvsiz),
107 . vxh(4,mvsiz),vyh(4,mvsiz),vzh(4,mvsiz)
108C-----------------------------------------------
109 adrbuf=igrsurf(ksurf)%IAD_BUFR
110 DO i=1,9
111 rot(i)=bufsf(adrbuf+7+i-1)
112 END DO
113C-----------------------------------------------
114C Noeud main dans le repere de l'ellipsoide :
115C-----------------------------------------------
116 x1=bufsf(adrbuf+16)-bufsf(adrbuf+4)
117 x2=bufsf(adrbuf+17)-bufsf(adrbuf+5)
118 x3=bufsf(adrbuf+18)-bufsf(adrbuf+6)
119 xm=rot(1)*x1+rot(2)*x2+rot(3)*x3
120 ym=rot(4)*x1+rot(5)*x2+rot(6)*x3
121 zm=rot(7)*x1+rot(8)*x2+rot(9)*x3
122C-----------------------------------------------
123C Vitesse du noeud main en local.
124C-----------------------------------------------
125 v1=bufsf(adrbuf+19)
126 v2=bufsf(adrbuf+20)
127 v3=bufsf(adrbuf+21)
128 vxm=rot(1)*v1+rot(2)*v2+rot(3)*v3
129 vym=rot(4)*v1+rot(5)*v2+rot(6)*v3
130 vzm=rot(7)*v1+rot(8)*v2+rot(9)*v3
131C-----------------------------------------------
132C Vitesse de rotation du noeud main en local.
133C-----------------------------------------------
134 v1=bufsf(adrbuf+22)
135 v2=bufsf(adrbuf+23)
136 v3=bufsf(adrbuf+24)
137 vrx=rot(1)*v1+rot(2)*v2+rot(3)*v3
138 vry=rot(4)*v1+rot(5)*v2+rot(6)*v3
139 vrz=rot(7)*v1+rot(8)*v2+rot(9)*v3
140C----------------------------------------------------------
141C PASSAGE DES VITESSES EN LOCAL.
142C non optimise (plusieurs fois / noeud).
143C----------------------------------------------------------
144 DO 100 nls=ndeb+1,min(ndeb+mvsiz,nsc)
145 il =ksc(nls)
146 i =nls-ndeb
147C------
148 in1=ksi(1,il)
149 v1=v(1,in1)
150 v2=v(2,in1)
151 v3=v(3,in1)
152 vx1(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
153 vy1(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
154 vz1(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
155C------
156 in2=ksi(2,il)
157 v1=v(1,in2)
158 v2=v(2,in2)
159 v3=v(3,in2)
160 vx2(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
161 vy2(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
162 vz2(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
163C------
164 in3=ksi(3,il)
165 v1=v(1,in3)
166 v2=v(2,in3)
167 v3=v(3,in3)
168 vx3(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
169 vy3(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
170 vz3(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
171C------
172 in4=ksi(4,il)
173 v1=v(1,in4)
174 v2=v(2,in4)
175 v3=v(3,in4)
176 vx4(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
177 vy4(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
178 vz4(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
179C------
180 vxg(i)=fourth*(vx1(i)+vx2(i)+vx3(i)+vx4(i))
181 vyg(i)=fourth*(vy1(i)+vy2(i)+vy3(i)+vy4(i))
182 vzg(i)=fourth*(vz1(i)+vz2(i)+vz3(i)+vz4(i))
183 100 CONTINUE
184C----------------------------------------------------------
185C 1ER TRIANGLE / Force normale.
186C----------------------------------------------------------
187 DO 200 nls=ndeb+1,min(ndeb+mvsiz,nsc)
188 il =ksc(nls)
189 IF (iactiv(1,il)<=0) GOTO 200
190 i =nls-ndeb
191C STF(1,I)=STFAC*DEPTH(1,I)/MAX(DEPTH(1,I)-PENET(1,I),EM20)
192 stf(1,i)=stfac*depth(1,i)**2/max((depth(1,i)-penet(1,i))**2,em20)
193 fxn(1,i)=stf(1,i)*penet(1,i)*nxi(1,i)
194 fyn(1,i)=stf(1,i)*penet(1,i)*nyi(1,i)
195 fzn(1,i)=stf(1,i)*penet(1,i)*nzi(1,i)
196 200 CONTINUE
197C-------------------------------
198C 1ER TRIANGLE / Extraction du pt de contact au cycle precedent.
199C-------------------------------
200 DO 250 nls=ndeb+1,min(ndeb+mvsiz,nsc)
201 il =ksc(nls)
202 IF (iactiv(1,il)<=0) GOTO 250
203 i =nls-ndeb
204C Point sur le triangle tq D(P,L) etait maximum au cycle precedent.
205 l1(i)=iold(1,4*(il-1)+1)
206 l2(i)=iold(2,4*(il-1)+1)
207 l3(i)=iold(3,4*(il-1)+1)
208 250 CONTINUE
209C----------------------------------------------------------
210C 1ER TRIANGLE / Frottement :
211C----------------------------------------------------------
212#include "vectorize.inc"
213 DO 275 nls=ndeb+1,min(ndeb+mvsiz,nsc)
214 il =ksc(nls)
215 IF (iactiv(1,il)<=0) GOTO 275
216 i =nls-ndeb
217C------------------------------
218 fxt(1,i)=zero
219 fyt(1,i)=zero
220 fzt(1,i)=zero
221C-------------------------------
222C TOURNER DELTA(T-1) DANS LE PLAN TANGENT.
223C-------------------------------
224 IF (iactiv(1,il)>2) THEN
225 deltx=dold(1,4*(il-1)+1)
226 delty=dold(2,4*(il-1)+1)
227 deltz=dold(3,4*(il-1)+1)
228 nd =sqrt(deltx*deltx+delty*delty+deltz*deltz)
229 IF (nd/=zero) THEN
230 pn=deltx*nxi(1,i)+delty*nyi(1,i)+deltz*nzi(1,i)
231 deltx=deltx-pn*nxi(1,i)
232 delty=delty-pn*nyi(1,i)
233 deltz=deltz-pn*nzi(1,i)
234 scale=nd/sqrt(deltx*deltx+delty*delty+deltz*deltz)
235 deltx=scale*deltx
236 delty=scale*delty
237 deltz=scale*deltz
238 ENDIF
239 ELSE
240 deltx=zero
241 delty=zero
242 deltz=zero
243 ENDIF
244C-------------------------------
245C INCREMENT SUR DELTA (T-1 -> T) : VITESSE RELATIVE TANGENTE DE POLD * DT
246C-------------------------------
247 xh=hold(1,4*(il-1)+1)
248 yh=hold(2,4*(il-1)+1)
249 zh=hold(3,4*(il-1)+1)
250 v1x2=vrx*(yh-ym)
251 v2x1=vry*(xh-xm)
252 v2x3=vry*(zh-zm)
253 v3x2=vrz*(yh-ym)
254 v3x1=vrz*(xh-xm)
255 v1x3=vrx*(zh-zm)
256 v3 =v1x2-v2x1
257 v1 =v2x3-v3x2
258 v2 =v3x1-v1x3
259 vxh(1,i)=vxm+v1
260 vyh(1,i)=vym+v2
261 vzh(1,i)=vzm+v3
262 IF (iactiv(1,il)>=2) THEN
263 vxk=l1(i)*vxg(i)+l2(i)*vx1(i)+l3(i)*vx2(i)
264 vyk=l1(i)*vyg(i)+l2(i)*vy1(i)+l3(i)*vy2(i)
265 vzk=l1(i)*vzg(i)+l2(i)*vz1(i)+l3(i)*vz2(i)
266 dvx=vxh(1,i)-vxk
267 dvy=vyh(1,i)-vyk
268 dvz=vzh(1,i)-vzk
269 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
270 dvx=dvx-pn*nxi(1,i)
271 dvy=dvy-pn*nyi(1,i)
272 dvz=dvz-pn*nzi(1,i)
273 ELSE
274 dvx=zero
275 dvy=zero
276 dvz=zero
277 ENDIF
278C-----
279C Force tangente K.DELTA
280 dold(1,4*(il-1)+1)=deltx+dvx*dt1
281 dold(2,4*(il-1)+1)=delty+dvy*dt1
282 dold(3,4*(il-1)+1)=deltz+dvz*dt1
283 fxt(1,i)=stf(1,i)*dold(1,4*(il-1)+1)
284 fyt(1,i)=stf(1,i)*dold(2,4*(il-1)+1)
285 fzt(1,i)=stf(1,i)*dold(3,4*(il-1)+1)
286C-----
287 275 CONTINUE
288C-------------------------------
289C 2EME TRIANGLE / Force normale.
290C-------------------------------
291 DO 300 nls=ndeb+1,min(ndeb+mvsiz,nsc)
292 il =ksc(nls)
293 IF (iactiv(2,il)<=0) GOTO 300
294 i =nls-ndeb
295C STF(2,I)=STFAC*DEPTH(2,I)/MAX(DEPTH(2,I)-PENET(2,I),EM20)
296 stf(2,i)=stfac*depth(2,i)**2/max((depth(2,i)-penet(2,i))**2,em20)
297 fxn(2,i)=stf(2,i)*penet(2,i)*nxi(2,i)
298 fyn(2,i)=stf(2,i)*penet(2,i)*nyi(2,i)
299 fzn(2,i)=stf(2,i)*penet(2,i)*nzi(2,i)
300 300 CONTINUE
301C-------------------------------
302C 2EME TRIANGLE / Extraction du pt de contact au cycle precedent.
303C-------------------------------
304 DO 350 nls=ndeb+1,min(ndeb+mvsiz,nsc)
305 il =ksc(nls)
306 IF (iactiv(2,il)<=0) GOTO 350
307 i =nls-ndeb
308C Point sur le triangle tq D(P,L) etait maximum au cycle precedent.
309 l1(i)=iold(1,4*(il-1)+2)
310 l2(i)=iold(2,4*(il-1)+2)
311 l3(i)=iold(3,4*(il-1)+2)
312 350 CONTINUE
313C----------------------------------------------------------
314C 2EME TRIANGLE / Frottement :
315C----------------------------------------------------------
316#include "vectorize.inc"
317 DO 375 nls=ndeb+1,min(ndeb+mvsiz,nsc)
318 il =ksc(nls)
319 IF (iactiv(2,il)<=0) GOTO 375
320 i =nls-ndeb
321 fxt(2,i)=zero
322 fyt(2,i)=zero
323 fzt(2,i)=zero
324C-------------------------------
325C TOURNER DELTA(T-1) DANS LE PLAN TANGENT.
326C-------------------------------
327 IF (iactiv(2,il)>2) THEN
328 deltx=dold(1,4*(il-1)+2)
329 delty=dold(2,4*(il-1)+2)
330 deltz=dold(3,4*(il-1)+2)
331 nd =sqrt(deltx*deltx+delty*delty+deltz*deltz)
332 IF (nd/=0.) THEN
333 pn=deltx*nxi(2,i)+delty*nyi(2,i)+deltz*nzi(2,i)
334 deltx=deltx-pn*nxi(2,i)
335 delty=delty-pn*nyi(2,i)
336 deltz=deltz-pn*nzi(2,i)
337 scale=nd/sqrt(deltx*deltx+delty*delty+deltz*deltz)
338 deltx=scale*deltx
339 delty=scale*delty
340 deltz=scale*deltz
341 ENDIF
342 ELSE
343 deltx=zero
344 delty=zero
345 deltz=zero
346 ENDIF
347C-------------------------------
348C INCREMENT SUR DELTA (T-1 -> T) : VITESSE RELATIVE TANGENTE DE POLD * DT
349C-------------------------------
350 xh=hold(1,4*(il-1)+2)
351 yh=hold(2,4*(il-1)+2)
352 zh=hold(3,4*(il-1)+2)
353 v1x2=vrx*(yh-ym)
354 v2x1=vry*(xh-xm)
355 v2x3=vry*(zh-zm)
356 v3x2=vrz*(yh-ym)
357 v3x1=vrz*(xh-xm)
358 v1x3=vrx*(zh-zm)
359 v3 =v1x2-v2x1
360 v1 =v2x3-v3x2
361 v2 =v3x1-v1x3
362 vxh(2,i)=vxm+v1
363 vyh(2,i)=vym+v2
364 vzh(2,i)=vzm+v3
365 IF (iactiv(2,il)>=2) THEN
366 vxk=l1(i)*vxg(i)+l2(i)*vx2(i)+l3(i)*vx3(i)
367 vyk=l1(i)*vyg(i)+l2(i)*vy2(i)+l3(i)*vy3(i)
368 vzk=l1(i)*vzg(i)+l2(i)*vz2(i)+l3(i)*vz3(i)
369 dvx=vxh(2,i)-vxk
370 dvy=vyh(2,i)-vyk
371 dvz=vzh(2,i)-vzk
372 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
373 dvx=dvx-pn*nxi(2,i)
374 dvy=dvy-pn*nyi(2,i)
375 dvz=dvz-pn*nzi(2,i)
376 ELSE
377 dvx=zero
378 dvy=zero
379 dvz=zero
380 ENDIF
381C-----
382C Force tangente K.DELTA, K=1. et Norme.
383 dold(1,4*(il-1)+2)=deltx+dvx*dt1
384 dold(2,4*(il-1)+2)=delty+dvy*dt1
385 dold(3,4*(il-1)+2)=deltz+dvz*dt1
386 fxt(2,i)=stf(2,i)*dold(1,4*(il-1)+2)
387 fyt(2,i)=stf(2,i)*dold(2,4*(il-1)+2)
388 fzt(2,i)=stf(2,i)*dold(3,4*(il-1)+2)
389C-----
390 375 CONTINUE
391C-------------------------------
392C 3EME TRIANGLE / Force normale.
393C-------------------------------
394 DO 400 nls=ndeb+1,min(ndeb+mvsiz,nsc)
395 il =ksc(nls)
396 IF (iactiv(3,il)<=0) GOTO 400
397 i =nls-ndeb
398C STF(3,I)=STFAC*DEPTH(3,I)/MAX(DEPTH(3,I)-PENET(3,I),EM20)
399 stf(3,i)=stfac*depth(3,i)**2/max((depth(3,i)-penet(3,i))**2,em20)
400 fxn(3,i)=stf(3,i)*penet(3,i)*nxi(3,i)
401 fyn(3,i)=stf(3,i)*penet(3,i)*nyi(3,i)
402 fzn(3,i)=stf(3,i)*penet(3,i)*nzi(3,i)
403 400 CONTINUE
404C-------------------------------
405C 3EME TRIANGLE / Extraction du pt de contact au cycle precedent.
406C-------------------------------
407 DO 450 nls=ndeb+1,min(ndeb+mvsiz,nsc)
408 il =ksc(nls)
409 IF (iactiv(3,il)<=0) GOTO 450
410 i =nls-ndeb
411C Point sur le triangle tq D(P,L) etait maximum au cycle precedent.
412 l1(i)=iold(1,4*(il-1)+3)
413 l2(i)=iold(2,4*(il-1)+3)
414 l3(i)=iold(3,4*(il-1)+3)
415 450 CONTINUE
416C----------------------------------------------------------
417C 3EME TRIANGLE / Frottement :
418C----------------------------------------------------------
419#include "vectorize.inc"
420 DO 475 nls=ndeb+1,min(ndeb+mvsiz,nsc)
421 il =ksc(nls)
422 IF (iactiv(3,il)<=0) GOTO 475
423 i =nls-ndeb
424C------------------------------
425 fxt(3,i)=zero
426 fyt(3,i)=zero
427 fzt(3,i)=zero
428C-------------------------------
429C TOURNER DELTA(T-1) DANS LE PLAN TANGENT.
430C-------------------------------
431 IF (iactiv(3,il)>2) THEN
432 deltx=dold(1,4*(il-1)+3)
433 delty=dold(2,4*(il-1)+3)
434 deltz=dold(3,4*(il-1)+3)
435 nd =sqrt(deltx*deltx+delty*delty+deltz*deltz)
436 IF (nd/=zero) THEN
437 pn=deltx*nxi(3,i)+delty*nyi(3,i)+deltz*nzi(3,i)
438 deltx=deltx-pn*nxi(3,i)
439 delty=delty-pn*nyi(3,i)
440 deltz=deltz-pn*nzi(3,i)
441 scale=nd/sqrt(deltx*deltx+delty*delty+deltz*deltz)
442 deltx=scale*deltx
443 delty=scale*delty
444 deltz=scale*deltz
445 ENDIF
446 ELSE
447 deltx=zero
448 delty=zero
449 deltz=zero
450 ENDIF
451C-------------------------------
452C INCREMENT SUR DELTA (T-1 -> T) : VITESSE RELATIVE TANGENTE DE POLD * DT
453C-------------------------------
454 xh=hold(1,4*(il-1)+3)
455 yh=hold(2,4*(il-1)+3)
456 zh=hold(3,4*(il-1)+3)
457 v1x2=vrx*(yh-ym)
458 v2x1=vry*(xh-xm)
459 v2x3=vry*(zh-zm)
460 v3x2=vrz*(yh-ym)
461 v3x1=vrz*(xh-xm)
462 v1x3=vrx*(zh-zm)
463 v3 =v1x2-v2x1
464 v1 =v2x3-v3x2
465 v2 =v3x1-v1x3
466 vxh(3,i)=vxm+v1
467 vyh(3,i)=vym+v2
468 vzh(3,i)=vzm+v3
469 IF (iactiv(3,il)>=2) THEN
470 vxk=l1(i)*vxg(i)+l2(i)*vx3(i)+l3(i)*vx4(i)
471 vyk=l1(i)*vyg(i)+l2(i)*vy3(i)+l3(i)*vy4(i)
472 vzk=l1(i)*vzg(i)+l2(i)*vz3(i)+l3(i)*vz4(i)
473 dvx=vxh(3,i)-vxk
474 dvy=vyh(3,i)-vyk
475 dvz=vzh(3,i)-vzk
476 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
477 dvx=dvx-pn*nxi(3,i)
478 dvy=dvy-pn*nyi(3,i)
479 dvz=dvz-pn*nzi(3,i)
480 ELSE
481 dvx=zero
482 dvy=zero
483 dvz=zero
484 ENDIF
485C-----
486C Force tangente K.DELTA, K=1. et Norme.
487 dold(1,4*(il-1)+3)=deltx+dvx*dt1
488 dold(2,4*(il-1)+3)=delty+dvy*dt1
489 dold(3,4*(il-1)+3)=deltz+dvz*dt1
490 fxt(3,i)=stf(3,i)*dold(1,4*(il-1)+3)
491 fyt(3,i)=stf(3,i)*dold(2,4*(il-1)+3)
492 fzt(3,i)=stf(3,i)*dold(3,4*(il-1)+3)
493C-----
494 475 CONTINUE
495C-------------------------------
496C 4EME TRIANGLE / Force normale.
497C-------------------------------
498 DO 500 nls=ndeb+1,min(ndeb+mvsiz,nsc)
499 il =ksc(nls)
500 IF (iactiv(4,il)<=0) GOTO 500
501 i =nls-ndeb
502C STF(4,I)=STFAC*DEPTH(4,I)/MAX(DEPTH(4,I)-PENET(4,I),EM20)
503 stf(4,i)=stfac*depth(4,i)**2/max((depth(4,i)-penet(4,i))**2,em20)
504 fxn(4,i)=stf(4,i)*penet(4,i)*nxi(4,i)
505 fyn(4,i)=stf(4,i)*penet(4,i)*nyi(4,i)
506 fzn(4,i)=stf(4,i)*penet(4,i)*nzi(4,i)
507 500 CONTINUE
508C-------------------------------
509C 4EME TRIANGLE / Extraction du pt de contact au cycle precedent.
510C-------------------------------
511 DO 550 nls=ndeb+1,min(ndeb+mvsiz,nsc)
512 il =ksc(nls)
513 IF (iactiv(4,il)<=0) GOTO 550
514 i =nls-ndeb
515C Point sur le triangle tq D(P,L) etait maximum au cycle precedent.
516 l1(i)=iold(1,4*(il-1)+4)
517 l2(i)=iold(2,4*(il-1)+4)
518 l3(i)=iold(3,4*(il-1)+4)
519 550 CONTINUE
520C----------------------------------------------------------
521C 4EME TRIANGLE / Frottement :
522C----------------------------------------------------------
523#include "vectorize.inc"
524 DO 575 nls=ndeb+1,min(ndeb+mvsiz,nsc)
525 il =ksc(nls)
526 IF (iactiv(4,il)<=0) GOTO 575
527 i =nls-ndeb
528C------------------------------
529 fxt(4,i)=zero
530 fyt(4,i)=zero
531 fzt(4,i)=zero
532C-------------------------------
533C TOURNER DELTA(T-1) DANS LE PLAN TANGENT.
534C-------------------------------
535 IF (iactiv(4,il)>2) THEN
536 deltx=dold(1,4*(il-1)+4)
537 delty=dold(2,4*(il-1)+4)
538 deltz=dold(3,4*(il-1)+4)
539 nd =sqrt(deltx*deltx+delty*delty+deltz*deltz)
540 IF (nd/=zero) THEN
541 pn=deltx*nxi(4,i)+delty*nyi(4,i)+deltz*nzi(4,i)
542 deltx=deltx-pn*nxi(4,i)
543 delty=delty-pn*nyi(4,i)
544 deltz=deltz-pn*nzi(4,i)
545 scale=nd/sqrt(deltx*deltx+delty*delty+deltz*deltz)
546 deltx=scale*deltx
547 delty=scale*delty
548 deltz=scale*deltz
549 ENDIF
550 ELSE
551 deltx=zero
552 delty=zero
553 deltz=zero
554 ENDIF
555C-------------------------------
556C INCREMENT SUR DELTA (T-1 -> T) : VITESSE RELATIVE TANGENTE DE POLD * DT
557C-------------------------------
558 xh=hold(1,4*(il-1)+4)
559 yh=hold(2,4*(il-1)+4)
560 zh=hold(3,4*(il-1)+4)
561 v1x2=vrx*(yh-ym)
562 v2x1=vry*(xh-xm)
563 v2x3=vry*(zh-zm)
564 v3x2=vrz*(yh-ym)
565 v3x1=vrz*(xh-xm)
566 v1x3=vrx*(zh-zm)
567 v3 =v1x2-v2x1
568 v1 =v2x3-v3x2
569 v2 =v3x1-v1x3
570 vxh(4,i)=vxm+v1
571 vyh(4,i)=vym+v2
572 vzh(4,i)=vzm+v3
573 IF (iactiv(4,il)>=2) THEN
574 vxk=l1(i)*vxg(i)+l2(i)*vx4(i)+l3(i)*vx1(i)
575 vyk=l1(i)*vyg(i)+l2(i)*vy4(i)+l3(i)*vy1(i)
576 vzk=l1(i)*vzg(i)+l2(i)*vz4(i)+l3(i)*vz1(i)
577 dvx=vxh(4,i)-vxk
578 dvy=vyh(4,i)-vyk
579 dvz=vzh(4,i)-vzk
580 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
581 dvx=dvx-pn*nxi(4,i)
582 dvy=dvy-pn*nyi(4,i)
583 dvz=dvz-pn*nzi(4,i)
584 ELSE
585 dvx=zero
586 dvy=zero
587 dvz=zero
588 ENDIF
589C-----
590C Force tangente K.DELTA, K=1. et Norme.
591 dold(1,4*(il-1)+4)=deltx+dvx*dt1
592 dold(2,4*(il-1)+4)=delty+dvy*dt1
593 dold(3,4*(il-1)+4)=deltz+dvz*dt1
594 fxt(4,i)=stf(4,i)*dold(1,4*(il-1)+4)
595 fyt(4,i)=stf(4,i)*dold(2,4*(il-1)+4)
596 fzt(4,i)=stf(4,i)*dold(3,4*(il-1)+4)
597C-----
598 575 CONTINUE
599C------------------------------------------------------------
600C FRICTION (LOCAL COULOMB FRICTION).
601C------------------------------------------------------------
602C-----------------------------------------------
603C 1ER TRIANGLE.
604C-----------------------------------------------
605#include "vectorize.inc"
606 DO 580 nls=ndeb+1,min(ndeb+mvsiz,nsc)
607 il =ksc(nls)
608 IF (iactiv(1,il)<=0) GOTO 580
609 i =nls-ndeb
610C-----
611 ftx=fxt(1,i)
612 fty=fyt(1,i)
613 ftz=fzt(1,i)
614 fnx=fxn(1,i)
615 fny=fyn(1,i)
616 fnz=fzn(1,i)
617 ftn=sqrt(ftx*ftx+fty*fty+ftz*ftz)
618 fn =sqrt(fnx*fnx+fny*fny+fnz*fnz)
619 fmax= max(vfric*fn,zero)
620 IF (ftn>fmax) THEN
621 scale =fmax/ftn
622 ftx=scale*ftx
623 fty=scale*fty
624 ftz=scale*ftz
625 ELSE
626 scale=one
627 ENDIF
628 IF (iactiv(1,il)>1) THEN
629C Glissement DELTA=scale.DELTA
630 deltx=scale*dold(1,4*(il-1)+1)
631 delty=scale*dold(2,4*(il-1)+1)
632 deltz=scale*dold(3,4*(il-1)+1)
633C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
634 dold(1,4*(il-1)+1)=deltx
635 dold(2,4*(il-1)+1)=delty
636 dold(3,4*(il-1)+1)=deltz
637 ENDIF
638 fxt(1,i)=ftx
639 fyt(1,i)=fty
640 fzt(1,i)=ftz
641 580 CONTINUE
642C-----------------------------------------------
643C 2EME TRIANGLE.
644C-----------------------------------------------
645#include "vectorize.inc"
646 DO 585 nls=ndeb+1,min(ndeb+mvsiz,nsc)
647 il =ksc(nls)
648 IF (iactiv(2,il)<=0) GOTO 585
649 i =nls-ndeb
650C-----
651 ftx=fxt(2,i)
652 fty=fyt(2,i)
653 ftz=fzt(2,i)
654 fnx=fxn(2,i)
655 fny=fyn(2,i)
656 fnz=fzn(2,i)
657 ftn=sqrt(ftx*ftx+fty*fty+ftz*ftz)
658 fn =sqrt(fnx*fnx+fny*fny+fnz*fnz)
659 fmax= max(vfric*fn,zero)
660 IF (ftn>fmax) THEN
661 scale =fmax/ftn
662 ftx=scale*ftx
663 fty=scale*fty
664 ftz=scale*ftz
665 ELSE
666 scale=one
667 ENDIF
668 IF (iactiv(2,il)>1) THEN
669C Glissement DELTA=scale.DELTA
670 deltx=scale*dold(1,4*(il-1)+2)
671 delty=scale*dold(2,4*(il-1)+2)
672 deltz=scale*dold(3,4*(il-1)+2)
673C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
674 dold(1,4*(il-1)+2)=deltx
675 dold(2,4*(il-1)+2)=delty
676 dold(3,4*(il-1)+2)=deltz
677 ENDIF
678 fxt(2,i)=ftx
679 fyt(2,i)=fty
680 fzt(2,i)=ftz
681 585 CONTINUE
682C-----------------------------------------------
683C 3EME TRIANGLE.
684C-----------------------------------------------
685#include "vectorize.inc"
686 DO 590 nls=ndeb+1,min(ndeb+mvsiz,nsc)
687 il =ksc(nls)
688 IF (iactiv(3,il)<=0) GOTO 590
689 i =nls-ndeb
690C-----
691 ftx=fxt(3,i)
692 fty=fyt(3,i)
693 ftz=fzt(3,i)
694 fnx=fxn(3,i)
695 fny=fyn(3,i)
696 fnz=fzn(3,i)
697 ftn=sqrt(ftx*ftx+fty*fty+ftz*ftz)
698 fn =sqrt(fnx*fnx+fny*fny+fnz*fnz)
699 fmax= max(vfric*fn,zero)
700 IF (ftn>fmax) THEN
701 scale =fmax/ftn
702 ftx=scale*ftx
703 fty=scale*fty
704 ftz=scale*ftz
705 ELSE
706 scale=one
707 ENDIF
708 IF (iactiv(3,il)>1) THEN
709C Glissement DELTA=scale.DELTA
710 deltx=scale*dold(1,4*(il-1)+3)
711 delty=scale*dold(2,4*(il-1)+3)
712 deltz=scale*dold(3,4*(il-1)+3)
713C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
714 dold(1,4*(il-1)+3)=deltx
715 dold(2,4*(il-1)+3)=delty
716 dold(3,4*(il-1)+3)=deltz
717 ENDIF
718 fxt(3,i)=ftx
719 fyt(3,i)=fty
720 fzt(3,i)=ftz
721 590 CONTINUE
722C-----------------------------------------------
723C 4EME TRIANGLE.
724C-----------------------------------------------
725#include "vectorize.inc"
726 DO 595 nls=ndeb+1,min(ndeb+mvsiz,nsc)
727 il =ksc(nls)
728 IF (iactiv(4,il)<=0) GOTO 595
729 i =nls-ndeb
730C-----
731 ftx=fxt(4,i)
732 fty=fyt(4,i)
733 ftz=fzt(4,i)
734 fnx=fxn(4,i)
735 fny=fyn(4,i)
736 fnz=fzn(4,i)
737 ftn=sqrt(ftx*ftx+fty*fty+ftz*ftz)
738 fn =sqrt(fnx*fnx+fny*fny+fnz*fnz)
739 fmax= max(vfric*fn,zero)
740 IF (ftn>fmax) THEN
741 scale =fmax/ftn
742 ftx=scale*ftx
743 fty=scale*fty
744 ftz=scale*ftz
745 ELSE
746 scale=one
747 ENDIF
748 IF (iactiv(4,il)>1) THEN
749C Glissement DELTA=scale.DELTA
750 deltx=scale*dold(1,4*(il-1)+4)
751 delty=scale*dold(2,4*(il-1)+4)
752 deltz=scale*dold(3,4*(il-1)+4)
753C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
754 dold(1,4*(il-1)+4)=deltx
755 dold(2,4*(il-1)+4)=delty
756 dold(3,4*(il-1)+4)=deltz
757 ENDIF
758 fxt(4,i)=ftx
759 fyt(4,i)=fty
760 fzt(4,i)=ftz
761 595 CONTINUE
762C------------------------------------------------------------
763C LOCAL COORDINATES OF CONTACT POINT WITH RESPECT TO ELEMENT.
764C------------------------------------------------------------
765C#include "vectorize.inc"
766 DO 650 nls=ndeb+1,min(ndeb+mvsiz,nsc)
767 il =ksc(nls)
768 IF (iactiv(1,il)<=0) GOTO 650
769 i =nls-ndeb
770C-----
771 dx2=xp1(1,i)-xtk(1,i)
772 dy2=xp1(2,i)-ytk(1,i)
773 dz2=xp1(3,i)-ztk(1,i)
774 dx3=xp2(1,i)-xtk(1,i)
775 dy3=xp2(2,i)-ytk(1,i)
776 dz3=xp2(3,i)-ztk(1,i)
777 dx=dy2*dz3-dz2*dy3
778 dy=dx3*dz2-dz3*dx2
779 dz=dx2*dy3-dy2*dx3
780 s1=half*sqrt(dx*dx+dy*dy+dz*dz)
781C-----
782 dx1=gx(1,i) -xtk(1,i)
783 dy1=gx(2,i) -ytk(1,i)
784 dz1=gx(3,i) -ztk(1,i)
785 dx=dy1*dz3-dz1*dy3
786 dy=dx3*dz1-dz3*dx1
787 dz=dx1*dy3-dy1*dx3
788 s2=half*sqrt(dx*dx+dy*dy+dz*dz)
789C-----
790 dx=dy1*dz2-dz1*dy2
791 dy=dx2*dz1-dz2*dx1
792 dz=dx1*dy2-dy1*dx2
793 s3=half*sqrt(dx*dx+dy*dy+dz*dz)
794C-----
795 s=one/(s1+s2+s3)
796C-----
797 iold(1,4*(il-1)+1)=s1*s
798 iold(2,4*(il-1)+1)=s2*s
799 iold(3,4*(il-1)+1)=s3*s
800 650 CONTINUE
801C------------------------------------------------------------
802C#include "vectorize.inc"
803 DO 660 nls=ndeb+1,min(ndeb+mvsiz,nsc)
804 il =ksc(nls)
805 IF (iactiv(2,il)<=0) GOTO 660
806 i =nls-ndeb
807C-----
808 dx2=xp2(1,i)-xtk(2,i)
809 dy2=xp2(2,i)-ytk(2,i)
810 dz2=xp2(3,i)-ztk(2,i)
811 dx3=xp3(1,i)-xtk(2,i)
812 dy3=xp3(2,i)-ytk(2,i)
813 dz3=xp3(3,i)-ztk(2,i)
814 dx=dy2*dz3-dz2*dy3
815 dy=dx3*dz2-dz3*dx2
816 dz=dx2*dy3-dy2*dx3
817 s1=half*sqrt(dx*dx+dy*dy+dz*dz)
818C-----
819 dx1=gx(1,i) -xtk(2,i)
820 dy1=gx(2,i) -ytk(2,i)
821 dz1=gx(3,i) -ztk(2,i)
822 dx=dy1*dz3-dz1*dy3
823 dy=dx3*dz1-dz3*dx1
824 dz=dx1*dy3-dy1*dx3
825 s2=half*sqrt(dx*dx+dy*dy+dz*dz)
826C-----
827 dx=dy1*dz2-dz1*dy2
828 dy=dx2*dz1-dz2*dx1
829 dz=dx1*dy2-dy1*dx2
830 s3=half*sqrt(dx*dx+dy*dy+dz*dz)
831C-----
832 s=one/(s1+s2+s3)
833C-----
834 iold(1,4*(il-1)+2)=s1*s
835 iold(2,4*(il-1)+2)=s2*s
836 iold(3,4*(il-1)+2)=s3*s
837 660 CONTINUE
838C------------------------------------------------------------
839C#include "vectorize.inc"
840 DO 670 nls=ndeb+1,min(ndeb+mvsiz,nsc)
841 il =ksc(nls)
842 IF (iactiv(3,il)<=0) GOTO 670
843 i =nls-ndeb
844C-----
845 dx2=xp3(1,i)-xtk(3,i)
846 dy2=xp3(2,i)-ytk(3,i)
847 dz2=xp3(3,i)-ztk(3,i)
848 dx3=xp4(1,i)-xtk(3,i)
849 dy3=xp4(2,i)-ytk(3,i)
850 dz3=xp4(3,i)-ztk(3,i)
851 dx=dy2*dz3-dz2*dy3
852 dy=dx3*dz2-dz3*dx2
853 dz=dx2*dy3-dy2*dx3
854 s1=half*sqrt(dx*dx+dy*dy+dz*dz)
855C-----
856 dx1=gx(1,i) -xtk(3,i)
857 dy1=gx(2,i) -ytk(3,i)
858 dz1=gx(3,i) -ztk(3,i)
859 dx=dy1*dz3-dz1*dy3
860 dy=dx3*dz1-dz3*dx1
861 dz=dx1*dy3-dy1*dx3
862 s2=half*sqrt(dx*dx+dy*dy+dz*dz)
863C-----
864 dx=dy1*dz2-dz1*dy2
865 dy=dx2*dz1-dz2*dx1
866 dz=dx1*dy2-dy1*dx2
867 s3=half*sqrt(dx*dx+dy*dy+dz*dz)
868C-----
869 s=one/(s1+s2+s3)
870C-----
871 iold(1,4*(il-1)+3)=s1*s
872 iold(2,4*(il-1)+3)=s2*s
873 iold(3,4*(il-1)+3)=s3*s
874 670 CONTINUE
875C------------------------------------------------------------
876C#include "vectorize.inc"
877 DO 680 nls=ndeb+1,min(ndeb+mvsiz,nsc)
878 il =ksc(nls)
879 IF (iactiv(4,il)<=0) GOTO 680
880 i =nls-ndeb
881C-----
882 dx2=xp4(1,i)-xtk(4,i)
883 dy2=xp4(2,i)-ytk(4,i)
884 dz2=xp4(3,i)-ztk(4,i)
885 dx3=xp1(1,i)-xtk(4,i)
886 dy3=xp1(2,i)-ytk(4,i)
887 dz3=xp1(3,i)-ztk(4,i)
888 dx=dy2*dz3-dz2*dy3
889 dy=dx3*dz2-dz3*dx2
890 dz=dx2*dy3-dy2*dx3
891 s1=half*sqrt(dx*dx+dy*dy+dz*dz)
892C-----
893 dx1=gx(1,i) -xtk(4,i)
894 dy1=gx(2,i) -ytk(4,i)
895 dz1=gx(3,i) -ztk(4,i)
896 dx=dy1*dz3-dz1*dy3
897 dy=dx3*dz1-dz3*dx1
898 dz=dx1*dy3-dy1*dx3
899 s2=half*sqrt(dx*dx+dy*dy+dz*dz)
900C-----
901 dx=dy1*dz2-dz1*dy2
902 dy=dx2*dz1-dz2*dx1
903 dz=dx1*dy2-dy1*dx2
904 s3=half*sqrt(dx*dx+dy*dy+dz*dz)
905C-----
906 s=one/(s1+s2+s3)
907C-----
908 iold(1,4*(il-1)+4)=s1*s
909 iold(2,4*(il-1)+4)=s2*s
910 iold(3,4*(il-1)+4)=s3*s
911 680 CONTINUE
912C------------------------------------------------------------
913C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
914C------------------------------------------------------------
915#include "vectorize.inc"
916 DO 700 nls=ndeb+1,min(ndeb+mvsiz,nsc)
917 il =ksc(nls)
918 i =nls-ndeb
919C-----
920 IF (iactiv(1,il)>0) THEN
921 hold(1,4*(il-1)+1)=xi(1,i)
922 hold(2,4*(il-1)+1)=yi(1,i)
923 hold(3,4*(il-1)+1)=zi(1,i)
924 nold(1,4*(il-1)+1)=nxi(1,i)
925 nold(2,4*(il-1)+1)=nyi(1,i)
926 nold(3,4*(il-1)+1)=nzi(1,i)
927 ENDIF
928 IF (iactiv(2,il)>0) THEN
929 hold(1,4*(il-1)+2)=xi(2,i)
930 hold(2,4*(il-1)+2)=yi(2,i)
931 hold(3,4*(il-1)+2)=zi(2,i)
932 nold(1,4*(il-1)+2)=nxi(2,i)
933 nold(2,4*(il-1)+2)=nyi(2,i)
934 nold(3,4*(il-1)+2)=nzi(2,i)
935 ENDIF
936 IF (iactiv(3,il)>0) THEN
937 hold(1,4*(il-1)+3)=xi(3,i)
938 hold(2,4*(il-1)+3)=yi(3,i)
939 hold(3,4*(il-1)+3)=zi(3,i)
940 nold(1,4*(il-1)+3)=nxi(3,i)
941 nold(2,4*(il-1)+3)=nyi(3,i)
942 nold(3,4*(il-1)+3)=nzi(3,i)
943 ENDIF
944 IF (iactiv(4,il)>0) THEN
945 hold(1,4*(il-1)+4)=xi(4,i)
946 hold(2,4*(il-1)+4)=yi(4,i)
947 hold(3,4*(il-1)+4)=zi(4,i)
948 nold(1,4*(il-1)+4)=nxi(4,i)
949 nold(2,4*(il-1)+4)=nyi(4,i)
950 nold(3,4*(il-1)+4)=nzi(4,i)
951 ENDIF
952 700 CONTINUE
953C------------------------------------------------------------
954C MISE EN MEMOIRE (VECTEURS DE TRAVAIL -> ASSEMBLAGE).
955C------------------------------------------------------------
956C non vectoriel.
957 DO 600 nls=ndeb+1,min(ndeb+mvsiz,nsc)
958 il =ksc(nls)
959 i =nls-ndeb
960C-----
961 in1=ksi(1,il)
962 in2=ksi(2,il)
963 in3=ksi(3,il)
964 in4=ksi(4,il)
965C-----
966 IF (iactiv(1,il)>0) THEN
967 st1=stf(1,i)
968 ELSE
969 st1=zero
970 ENDIF
971 IF (iactiv(2,il)>0) THEN
972 st2=stf(2,i)
973 ELSE
974 st2=zero
975 ENDIF
976 IF (iactiv(3,il)>0) THEN
977 st3=stf(3,i)
978 ELSE
979 st3=zero
980 ENDIF
981 IF (iactiv(4,il)>0) THEN
982 st4=stf(4,i)
983 ELSE
984 st4=zero
985 ENDIF
986 st1g=iold(1,4*(il-1)+1)*st1
987 st11=iold(2,4*(il-1)+1)*st1
988 st12=iold(3,4*(il-1)+1)*st1
989 st2g=iold(1,4*(il-1)+2)*st2
990 st22=iold(2,4*(il-1)+2)*st2
991 st23=iold(3,4*(il-1)+2)*st2
992 st3g=iold(1,4*(il-1)+3)*st3
993 st33=iold(2,4*(il-1)+3)*st3
994 st34=iold(3,4*(il-1)+3)*st3
995 st4g=iold(1,4*(il-1)+4)*st4
996 st44=iold(2,4*(il-1)+4)*st4
997 st41=iold(3,4*(il-1)+4)*st4
998 stg =fourth*(st1g+st2g+st3g+st4g)
999 wns(in1)=wns(in1)+st11+st41+stg
1000 wns(in2)=wns(in2)+st12+st22+stg
1001 wns(in3)=wns(in3)+st23+st33+stg
1002 wns(in4)=wns(in4)+st34+st44+stg
1003C-----
1004 600 CONTINUE
1005C------------------------------------------------------------
1006C non vectoriel.
1007 DO 800 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1008 il =ksc(nls)
1009 i =nls-ndeb
1010C-----
1011 in1=ksi(1,il)
1012 in2=ksi(2,il)
1013 in3=ksi(3,il)
1014 in4=ksi(4,il)
1015C-----
1016 IF (iactiv(1,il)>0) THEN
1017 f1=fxn(1,i)
1018 ELSE
1019 f1=zero
1020 ENDIF
1021 IF (iactiv(2,il)>0) THEN
1022 f2=fxn(2,i)
1023 ELSE
1024 f2=zero
1025 ENDIF
1026 IF (iactiv(3,il)>0) THEN
1027 f3=fxn(3,i)
1028 ELSE
1029 f3=zero
1030 ENDIF
1031 IF (iactiv(4,il)>0) THEN
1032 f4=fxn(4,i)
1033 ELSE
1034 f4=zero
1035 ENDIF
1036C-----
1037 fnormx=fnormx+f1+f2+f3+f4
1038C-----
1039 f1g=iold(1,4*(il-1)+1)*f1
1040 f11=iold(2,4*(il-1)+1)*f1
1041 f12=iold(3,4*(il-1)+1)*f1
1042 f2g=iold(1,4*(il-1)+2)*f2
1043 f22=iold(2,4*(il-1)+2)*f2
1044 f23=iold(3,4*(il-1)+2)*f2
1045 f3g=iold(1,4*(il-1)+3)*f3
1046 f33=iold(2,4*(il-1)+3)*f3
1047 f34=iold(3,4*(il-1)+3)*f3
1048 f4g=iold(1,4*(il-1)+4)*f4
1049 f44=iold(2,4*(il-1)+4)*f4
1050 f41=iold(3,4*(il-1)+4)*f4
1051 fg =fourth*(f1g+f2g+f3g+f4g)
1052 wnf(1,in1)=wnf(1,in1)+f11+f41+fg
1053 wnf(1,in2)=wnf(1,in2)+f12+f22+fg
1054 wnf(1,in3)=wnf(1,in3)+f23+f33+fg
1055 wnf(1,in4)=wnf(1,in4)+f34+f44+fg
1056C-----
1057 IF (iactiv(1,il)>0) THEN
1058 f1=fyn(1,i)
1059 ELSE
1060 f1=zero
1061 ENDIF
1062 IF (iactiv(2,il)>0) THEN
1063 f2=fyn(2,i)
1064 ELSE
1065 f2=zero
1066 ENDIF
1067 IF (iactiv(3,il)>0) THEN
1068 f3=fyn(3,i)
1069 ELSE
1070 f3=zero
1071 ENDIF
1072 IF (iactiv(4,il)>0) THEN
1073 f4=fyn(4,i)
1074 ELSE
1075 f4=zero
1076 ENDIF
1077C-----
1078 fnormy=fnormy+f1+f2+f3+f4
1079C-----
1080 f1g=iold(1,4*(il-1)+1)*f1
1081 f11=iold(2,4*(il-1)+1)*f1
1082 f12=iold(3,4*(il-1)+1)*f1
1083 f2g=iold(1,4*(il-1)+2)*f2
1084 f22=iold(2,4*(il-1)+2)*f2
1085 f23=iold(3,4*(il-1)+2)*f2
1086 f3g=iold(1,4*(il-1)+3)*f3
1087 f33=iold(2,4*(il-1)+3)*f3
1088 f34=iold(3,4*(il-1)+3)*f3
1089 f4g=iold(1,4*(il-1)+4)*f4
1090 f44=iold(2,4*(il-1)+4)*f4
1091 f41=iold(3,4*(il-1)+4)*f4
1092 fg =fourth*(f1g+f2g+f3g+f4g)
1093 wnf(2,in1)=wnf(2,in1)+f11+f41+fg
1094 wnf(2,in2)=wnf(2,in2)+f12+f22+fg
1095 wnf(2,in3)=wnf(2,in3)+f23+f33+fg
1096 wnf(2,in4)=wnf(2,in4)+f34+f44+fg
1097C-----
1098 IF (iactiv(1,il)>0) THEN
1099 f1=fzn(1,i)
1100 ELSE
1101 f1=zero
1102 ENDIF
1103 IF (iactiv(2,il)>0) THEN
1104 f2=fzn(2,i)
1105 ELSE
1106 f2=zero
1107 ENDIF
1108 IF (iactiv(3,il)>0) THEN
1109 f3=fzn(3,i)
1110 ELSE
1111 f3=zero
1112 ENDIF
1113 IF (iactiv(4,il)>0) THEN
1114 f4=fzn(4,i)
1115 ELSE
1116 f4=zero
1117 ENDIF
1118C-----
1119 fnormz=fnormz+f1+f2+f3+f4
1120C-----
1121 f1g=iold(1,4*(il-1)+1)*f1
1122 f11=iold(2,4*(il-1)+1)*f1
1123 f12=iold(3,4*(il-1)+1)*f1
1124 f2g=iold(1,4*(il-1)+2)*f2
1125 f22=iold(2,4*(il-1)+2)*f2
1126 f23=iold(3,4*(il-1)+2)*f2
1127 f3g=iold(1,4*(il-1)+3)*f3
1128 f33=iold(2,4*(il-1)+3)*f3
1129 f34=iold(3,4*(il-1)+3)*f3
1130 f4g=iold(1,4*(il-1)+4)*f4
1131 f44=iold(2,4*(il-1)+4)*f4
1132 f41=iold(3,4*(il-1)+4)*f4
1133 fg =fourth*(f1g+f2g+f3g+f4g)
1134 wnf(3,in1)=wnf(3,in1)+f11+f41+fg
1135 wnf(3,in2)=wnf(3,in2)+f12+f22+fg
1136 wnf(3,in3)=wnf(3,in3)+f23+f33+fg
1137 wnf(3,in4)=wnf(3,in4)+f34+f44+fg
1138C-----
1139 800 CONTINUE
1140C------------------------------------------------------------
1141C non vectoriel.
1142 DO 825 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1143 il =ksc(nls)
1144 i =nls-ndeb
1145C-----
1146 in1=ksi(1,il)
1147 in2=ksi(2,il)
1148 in3=ksi(3,il)
1149 in4=ksi(4,il)
1150C-----
1151 IF (iactiv(1,il)>0) THEN
1152 f1=fxt(1,i)
1153 ELSE
1154 f1=zero
1155 ENDIF
1156 IF (iactiv(2,il)>0) THEN
1157 f2=fxt(2,i)
1158 ELSE
1159 f2=zero
1160 ENDIF
1161 IF (iactiv(3,il)>0) THEN
1162 f3=fxt(3,i)
1163 ELSE
1164 f3=zero
1165 ENDIF
1166 IF (iactiv(4,il)>0) THEN
1167 f4=fxt(4,i)
1168 ELSE
1169 f4=zero
1170 ENDIF
1171C-----
1172 ftangx=ftangx+f1+f2+f3+f4
1173C-----
1174 f1g=iold(1,4*(il-1)+1)*f1
1175 f11=iold(2,4*(il-1)+1)*f1
1176 f12=iold(3,4*(il-1)+1)*f1
1177 f2g=iold(1,4*(il-1)+2)*f2
1178 f22=iold(2,4*(il-1)+2)*f2
1179 f23=iold(3,4*(il-1)+2)*f2
1180 f3g=iold(1,4*(il-1)+3)*f3
1181 f33=iold(2,4*(il-1)+3)*f3
1182 f34=iold(3,4*(il-1)+3)*f3
1183 f4g=iold(1,4*(il-1)+4)*f4
1184 f44=iold(2,4*(il-1)+4)*f4
1185 f41=iold(3,4*(il-1)+4)*f4
1186 fg =fourth*(f1g+f2g+f3g+f4g)
1187 wtf(1,in1)=wtf(1,in1)+f11+f41+fg
1188 wtf(1,in2)=wtf(1,in2)+f12+f22+fg
1189 wtf(1,in3)=wtf(1,in3)+f23+f33+fg
1190 wtf(1,in4)=wtf(1,in4)+f34+f44+fg
1191C-----
1192 IF (iactiv(1,il)>0) THEN
1193 f1=fyt(1,i)
1194 ELSE
1195 f1=zero
1196 ENDIF
1197 IF (iactiv(2,il)>0) THEN
1198 f2=fyt(2,i)
1199 ELSE
1200 f2=zero
1201 ENDIF
1202 IF (iactiv(3,il)>0) THEN
1203 f3=fyt(3,i)
1204 ELSE
1205 f3=zero
1206 ENDIF
1207 IF (iactiv(4,il)>0) THEN
1208 f4=fyt(4,i)
1209 ELSE
1210 f4=zero
1211 ENDIF
1212C-----
1213 ftangy=ftangy+f1+f2+f3+f4
1214C-----
1215 f1g=iold(1,4*(il-1)+1)*f1
1216 f11=iold(2,4*(il-1)+1)*f1
1217 f12=iold(3,4*(il-1)+1)*f1
1218 f2g=iold(1,4*(il-1)+2)*f2
1219 f22=iold(2,4*(il-1)+2)*f2
1220 f23=iold(3,4*(il-1)+2)*f2
1221 f3g=iold(1,4*(il-1)+3)*f3
1222 f33=iold(2,4*(il-1)+3)*f3
1223 f34=iold(3,4*(il-1)+3)*f3
1224 f4g=iold(1,4*(il-1)+4)*f4
1225 f44=iold(2,4*(il-1)+4)*f4
1226 f41=iold(3,4*(il-1)+4)*f4
1227 fg =fourth*(f1g+f2g+f3g+f4g)
1228 wtf(2,in1)=wtf(2,in1)+f11+f41+fg
1229 wtf(2,in2)=wtf(2,in2)+f12+f22+fg
1230 wtf(2,in3)=wtf(2,in3)+f23+f33+fg
1231 wtf(2,in4)=wtf(2,in4)+f34+f44+fg
1232C-----
1233 IF (iactiv(1,il)>0) THEN
1234 f1=fzt(1,i)
1235 ELSE
1236 f1=zero
1237 ENDIF
1238 IF (iactiv(2,il)>0) THEN
1239 f2=fzt(2,i)
1240 ELSE
1241 f2=zero
1242 ENDIF
1243 IF (iactiv(3,il)>0) THEN
1244 f3=fzt(3,i)
1245 ELSE
1246 f3=zero
1247 ENDIF
1248 IF (iactiv(4,il)>0) THEN
1249 f4=fzt(4,i)
1250 ELSE
1251 f4=zero
1252 ENDIF
1253C-----
1254 ftangz=ftangz+f1+f2+f3+f4
1255C-----
1256 f1g=iold(1,4*(il-1)+1)*f1
1257 f11=iold(2,4*(il-1)+1)*f1
1258 f12=iold(3,4*(il-1)+1)*f1
1259 f2g=iold(1,4*(il-1)+2)*f2
1260 f22=iold(2,4*(il-1)+2)*f2
1261 f23=iold(3,4*(il-1)+2)*f2
1262 f3g=iold(1,4*(il-1)+3)*f3
1263 f33=iold(2,4*(il-1)+3)*f3
1264 f34=iold(3,4*(il-1)+3)*f3
1265 f4g=iold(1,4*(il-1)+4)*f4
1266 f44=iold(2,4*(il-1)+4)*f4
1267 f41=iold(3,4*(il-1)+4)*f4
1268 fg =fourth*(f1g+f2g+f3g+f4g)
1269 wtf(3,in1)=wtf(3,in1)+f11+f41+fg
1270 wtf(3,in2)=wtf(3,in2)+f12+f22+fg
1271 wtf(3,in3)=wtf(3,in3)+f23+f33+fg
1272 wtf(3,in4)=wtf(3,in4)+f34+f44+fg
1273C-----
1274 825 CONTINUE
1275C------------------------------------------------------------
1276C KINEMATIC TIME STEP.
1277C------------------------------------------------------------
1278 dti=ep20
1279 DO 900 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1280 il =ksc(nls)
1281 i =nls-ndeb
1282 IF (iactiv(1,il)>0) THEN
1283 dvx=vx1(i)-vxh(1,i)
1284 dvy=vy1(i)-vyh(1,i)
1285 dvz=vz1(i)-vzh(1,i)
1286 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
1287 IF (pn<zero) THEN
1288 h =xp1(1,i)*nxi(1,i)+xp1(2,i)*nyi(1,i)+xp1(3,i)*nzi(1,i)
1289 dti=min(dti,half*h/max(em20,-pn))
1290 ENDIF
1291 dvx=vx2(i)-vxh(1,i)
1292 dvy=vy2(i)-vyh(1,i)
1293 dvz=vz2(i)-vzh(1,i)
1294 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
1295 IF (pn<zero) THEN
1296 h =xp2(1,i)*nxi(1,i)+xp2(2,i)*nyi(1,i)+xp2(3,i)*nzi(1,i)
1297 dti=min(dti,half*h/max(em20,-pn))
1298 ENDIF
1299 dvx=vxg(i)-vxh(1,i)
1300 dvy=vyg(i)-vyh(1,i)
1301 dvz=vzg(i)-vzh(1,i)
1302 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
1303 IF (pn<zero) THEN
1304 h =gx(1,i)*nxi(1,i)+gx(2,i)*nyi(1,i)+gx(3,i)*nzi(1,i)
1305 dti=min(dti,half*h/max(em20,-pn))
1306 ENDIF
1307 ENDIF
1308 IF (iactiv(2,il)>0) THEN
1309 dvx=vx2(i)-vxh(2,i)
1310 dvy=vy2(i)-vyh(2,i)
1311 dvz=vz2(i)-vzh(2,i)
1312 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
1313 IF (pn<zero) THEN
1314 h =xp2(1,i)*nxi(2,i)+xp2(2,i)*nyi(2,i)+xp2(3,i)*nzi(2,i)
1315 dti=min(dti,half*h/max(em20,-pn))
1316 ENDIF
1317 dvx=vx3(i)-vxh(2,i)
1318 dvy=vy3(i)-vyh(2,i)
1319 dvz=vz3(i)-vzh(2,i)
1320 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
1321 IF (pn<zero) THEN
1322 h =xp3(1,i)*nxi(2,i)+xp3(2,i)*nyi(2,i)+xp3(3,i)*nzi(2,i)
1323 dti=min(dti,half*h/max(em20,-pn))
1324 ENDIF
1325 dvx=vxg(i)-vxh(2,i)
1326 dvy=vyg(i)-vyh(2,i)
1327 dvz=vzg(i)-vzh(2,i)
1328 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
1329 IF (pn<zero) THEN
1330 h =gx(1,i)*nxi(2,i)+gx(2,i)*nyi(2,i)+gx(3,i)*nzi(2,i)
1331 dti=min(dti,half*h/max(em20,-pn))
1332 ENDIF
1333 ENDIF
1334 IF (iactiv(3,il)>0) THEN
1335 dvx=vx3(i)-vxh(3,i)
1336 dvy=vy3(i)-vyh(3,i)
1337 dvz=vz3(i)-vzh(3,i)
1338 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
1339 IF (pn<zero) THEN
1340 h =xp3(1,i)*nxi(3,i)+xp3(2,i)*nyi(3,i)+xp3(3,i)*nzi(3,i)
1341 dti=min(dti,half*h/max(em20,-pn))
1342 ENDIF
1343 dvx=vx4(i)-vxh(3,i)
1344 dvy=vy4(i)-vyh(3,i)
1345 dvz=vz4(i)-vzh(3,i)
1346 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
1347 IF (pn<zero) THEN
1348 h =xp4(1,i)*nxi(3,i)+xp4(2,i)*nyi(3,i)+xp4(3,i)*nzi(3,i)
1349 dti=min(dti,half*h/max(em20,-pn))
1350 ENDIF
1351 dvx=vxg(i)-vxh(3,i)
1352 dvy=vyg(i)-vyh(3,i)
1353 dvz=vzg(i)-vzh(3,i)
1354 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
1355 IF (pn<zero) THEN
1356 h =gx(1,i)*nxi(3,i)+gx(2,i)*nyi(3,i)+gx(3,i)*nzi(3,i)
1357 dti=min(dti,half*h/max(em20,-pn))
1358 ENDIF
1359 ENDIF
1360 IF (iactiv(4,il)>0) THEN
1361 dvx=vx1(i)-vxh(4,i)
1362 dvy=vy1(i)-vyh(4,i)
1363 dvz=vz1(i)-vzh(4,i)
1364 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
1365 IF (pn<zero) THEN
1366 h =xp1(1,i)*nxi(4,i)+xp1(2,i)*nyi(4,i)+xp1(3,i)*nzi(4,i)
1367 dti=min(dti,half*h/max(em20,-pn))
1368 ENDIF
1369 dvx=vx4(i)-vxh(4,i)
1370 dvy=vy4(i)-vyh(4,i)
1371 dvz=vz4(i)-vzh(4,i)
1372 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
1373 IF (pn<zero) THEN
1374 h =xp4(1,i)*nxi(4,i)+xp4(2,i)*nyi(4,i)+xp4(3,i)*nzi(4,i)
1375 dti=min(dti,half*h/max(em20,-pn))
1376 ENDIF
1377 dvx=vxg(i)-vxh(4,i)
1378 dvy=vyg(i)-vyh(4,i)
1379 dvz=vzg(i)-vzh(4,i)
1380 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
1381 IF (pn<zero) THEN
1382 h =gx(1,i)*nxi(4,i)+gx(2,i)*nyi(4,i)+gx(3,i)*nzi(4,i)
1383 dti=min(dti,half*h/max(em20,-pn))
1384 ENDIF
1385 ENDIF
1386 IF(dti<dt2t)THEN
1387 dt2t = dti
1388 neltst = noint
1389 ityptst = 10
1390 ENDIF
1391C-----
1392 900 CONTINUE
1393C----------------------------------
1394 IF (dti<=zero) THEN
1395 dti=ep20
1396 DO 950 nls=ndeb+1,min(ndeb+mvsiz,nsc)
1397 il =ksc(nls)
1398 i =nls-ndeb
1399 IF (iactiv(1,il)>0) THEN
1400 dvx=vx1(i)-vxh(1,i)
1401 dvy=vy1(i)-vyh(1,i)
1402 dvz=vz1(i)-vzh(1,i)
1403 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
1404 IF (pn<zero) THEN
1405 h =xp1(1,i)*nxi(1,i)+xp1(2,i)*nyi(1,i)+xp1(3,i)*nzi(1,i)
1406 dti=min(dti,half*h/max(em20,-pn))
1407 ENDIF
1408 dvx=vx2(i)-vxh(1,i)
1409 dvy=vy2(i)-vyh(1,i)
1410 dvz=vz2(i)-vzh(1,i)
1411 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
1412 IF (pn<zero) THEN
1413 h =xp2(1,i)*nxi(1,i)+xp2(2,i)*nyi(1,i)+xp2(3,i)*nzi(1,i)
1414 dti=min(dti,half*h/max(em20,-pn))
1415 ENDIF
1416 dvx=vxg(i)-vxh(1,i)
1417 dvy=vyg(i)-vyh(1,i)
1418 dvz=vzg(i)-vzh(1,i)
1419 pn=dvx*nxi(1,i)+dvy*nyi(1,i)+dvz*nzi(1,i)
1420 IF (pn<zero) THEN
1421 h =gx(1,i)*nxi(1,i)+gx(2,i)*nyi(1,i)+gx(3,i)*nzi(1,i)
1422 dti=min(dti,half*h/max(em20,-pn))
1423 ENDIF
1424 ENDIF
1425 IF (iactiv(2,il)>0) THEN
1426 dvx=vx2(i)-vxh(2,i)
1427 dvy=vy2(i)-vyh(2,i)
1428 dvz=vz2(i)-vzh(2,i)
1429 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
1430 IF (pn<zero) THEN
1431 h =xp2(1,i)*nxi(2,i)+xp2(2,i)*nyi(2,i)+xp2(3,i)*nzi(2,i)
1432 dti=min(dti,half*h/max(em20,-pn))
1433 ENDIF
1434 dvx=vx3(i)-vxh(2,i)
1435 dvy=vy3(i)-vyh(2,i)
1436 dvz=vz3(i)-vzh(2,i)
1437 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
1438 IF (pn<zero) THEN
1439 h =xp3(1,i)*nxi(2,i)+xp3(2,i)*nyi(2,i)+xp3(3,i)*nzi(2,i)
1440 dti=min(dti,half*h/max(em20,-pn))
1441 ENDIF
1442 dvx=vxg(i)-vxh(2,i)
1443 dvy=vyg(i)-vyh(2,i)
1444 dvz=vzg(i)-vzh(2,i)
1445 pn=dvx*nxi(2,i)+dvy*nyi(2,i)+dvz*nzi(2,i)
1446 IF (pn<zero) THEN
1447 h =gx(1,i)*nxi(2,i)+gx(2,i)*nyi(2,i)+gx(3,i)*nzi(2,i)
1448 dti=min(dti,half*h/max(em20,-pn))
1449 ENDIF
1450 ENDIF
1451 IF (iactiv(3,il)>0) THEN
1452 dvx=vx3(i)-vxh(3,i)
1453 dvy=vy3(i)-vyh(3,i)
1454 dvz=vz3(i)-vzh(3,i)
1455 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
1456 IF (pn<zero) THEN
1457 h =xp3(1,i)*nxi(3,i)+xp3(2,i)*nyi(3,i)+xp3(3,i)*nzi(3,i)
1458 dti=min(dti,half*h/max(em20,-pn))
1459 ENDIF
1460 dvx=vx4(i)-vxh(3,i)
1461 dvy=vy4(i)-vyh(3,i)
1462 dvz=vz4(i)-vzh(3,i)
1463 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
1464 IF (pn<zero) THEN
1465 h =xp4(1,i)*nxi(3,i)+xp4(2,i)*nyi(3,i)+xp4(3,i)*nzi(3,i)
1466 dti=min(dti,half*h/max(em20,-pn))
1467 ENDIF
1468 dvx=vxg(i)-vxh(3,i)
1469 dvy=vyg(i)-vyh(3,i)
1470 dvz=vzg(i)-vzh(3,i)
1471 pn=dvx*nxi(3,i)+dvy*nyi(3,i)+dvz*nzi(3,i)
1472 IF (pn<zero) THEN
1473 h =gx(1,i)*nxi(3,i)+gx(2,i)*nyi(3,i)+gx(3,i)*nzi(3,i)
1474 dti=min(dti,half*h/max(em20,-pn))
1475 ENDIF
1476 ENDIF
1477 IF (iactiv(4,il)>0) THEN
1478 dvx=vx1(i)-vxh(4,i)
1479 dvy=vy1(i)-vyh(4,i)
1480 dvz=vz1(i)-vzh(4,i)
1481 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
1482 IF (pn<zero) THEN
1483 h =xp1(1,i)*nxi(4,i)+xp1(2,i)*nyi(4,i)+xp1(3,i)*nzi(4,i)
1484 dti=min(dti,half*h/max(em20,-pn))
1485 ENDIF
1486 dvx=vx4(i)-vxh(4,i)
1487 dvy=vy4(i)-vyh(4,i)
1488 dvz=vz4(i)-vzh(4,i)
1489 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
1490 IF (pn<zero) THEN
1491 h =xp4(1,i)*nxi(4,i)+xp4(2,i)*nyi(4,i)+xp4(3,i)*nzi(4,i)
1492 dti=min(dti,half*h/max(em20,-pn))
1493 ENDIF
1494 dvx=vxg(i)-vxh(4,i)
1495 dvy=vyg(i)-vyh(4,i)
1496 dvz=vzg(i)-vzh(4,i)
1497 pn=dvx*nxi(4,i)+dvy*nyi(4,i)+dvz*nzi(4,i)
1498 IF (pn<zero) THEN
1499 h =gx(1,i)*nxi(4,i)+gx(2,i)*nyi(4,i)+gx(3,i)*nzi(4,i)
1500 dti=min(dti,half*h/max(em20,-pn))
1501 ENDIF
1502 ENDIF
1503 IF(dti<=zero)THEN
1504 in1=ksi(1,il)
1505 in2=ksi(2,il)
1506 in3=ksi(3,il)
1507 in4=ksi(4,il)
1508#include "lockon.inc"
1509 WRITE(istdo,'(A,E12.4,A,I8)')
1510 . ' **WARNING MINIMUM TIME STEP ',dti,' IN INTERFACE ',noint
1511 WRITE(iout ,'(A,E12.4,A,I8)')
1512 . ' **WARNING MINIMUM TIME STEP ',dti,' IN INTERFACE ',noint
1513 WRITE(iout,'(A,4I8)') ' ELEMENT/SEGMENT NODES :',
1514 . in1,in2,in3,in4
1515#include "lockoff.inc"
1516 tstop = tt
1517 ENDIF
1518C-----
1519 950 CONTINUE
1520 ENDIF
1521C----------------------------------
1522 RETURN
1523 END
i15for1
subroutine i15for1(ndeb, nsc, stfac, x, v, ksurf, igrsurf, bufsf, ksc, ksi, iactiv, iold, hold, nold, dold, xp1, xp2, xp3, xp4, gx, xtk, ytk, ztk, ntx, nty, ntz, penet, depth, xi, yi, zi, nxi, nyi, nzi, ms, de, npc, pld, wnf, wtf, wns, fnormx, fnormy, fnormz, ftangx, ftangy, ftangz, dt2t, noint, neltst, ityptst, vfric)
Definition i15for1.F:41
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