OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i15fort1.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!|| i15fort1 ../engine/source/interfaces/int15/i15fort1.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 i15fort1(NDEB , NTC ,STFAC ,X ,V ,
31 2 KSURF ,IGRSURF ,BUFSF ,KTC ,KSI ,
32 3 IACTIV ,IOLD ,HOLD ,NOLD ,DOLD ,
33 4 XP1 ,XP2 ,XP3 ,XTK ,YTK ,
34 5 ZTK ,NTX ,NTY ,NTZ ,PENET,
35 6 DEPTH ,XI ,YI ,ZI ,NXI ,
36 7 NYI ,NZI ,MS ,DE ,NPC ,
37 8 PLD ,WNF ,WTF ,WNS ,FNORMX,
38 9 FNORMY,FNORMZ,FTANGX,FTANGY,FTANGZ ,
39 A DT2T ,NOINT , NELTST ,ITYPTST ,VFRIC )
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE groupdef_mod
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "units_c.inc"
59C-----------------------------------------------
60C D u m m y A r g u m e n t s
61C-----------------------------------------------
62 INTEGER NDEB, NTC
63 INTEGER KSURF,
64 . KSI(4,*),IACTIV(4,*),NPC(*),KTC(*),
65 . NOINT,NELTST,ITYPTST
66C REAL
67 my_real STFAC, BUFSF(*),X(3,*) , IOLD(3,*),
68 . HOLD(3,*) , NOLD(3,*) ,DOLD(3,*),MS(*) , V(3,*),
69 . DE, PLD(*), XTK(*) ,YTK(*) ,ZTK(*) ,
70 . NTX(*) ,NTY(*) ,NTZ(*) ,
71 . penet(*) ,depth(*),
72 . xi(*) ,yi(*) ,zi(*) ,nxi(*) ,nyi(*) ,nzi(*) ,
73 . wnf(3,*) ,wtf(3,*) ,wns(*) ,xp1(3,*) ,xp2(3,*) ,xp3(3,*) ,
74 . fnormx,fnormy,fnormz,ftangx,ftangy,ftangz,
75 . dt2t,vfric
76 TYPE (SURF_) , DIMENSION(NSURF) :: IGRSURF
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER ADRBUF, I, IL, IN1, IN2, IN3, IN4, NLS
81 my_real
82 . A, B, C, ROT(9), X1, X2, X3, XM, YM, ZM,
83 . V1, V2, V3, VXM, VYM, VZM, VRX, VRY, VRZ,
84 . V1X2, V2X1, V1X3, V3X1, V2X3, V3X2,
85 . deltx, delty, deltz, nd, pn, scale, nv,
86 . vxk, vyk, vzk, dvx, dvy, dvz, xh, yh, zh,
87 . dx1, dy1, dz1, dx2, dy2, dz2, dx3, dy3, dz3,
88 . dx, dy, dz,
89 . s1, s2, s3, s,
90 . f1, f11, f12, f13,
91 . ftx, fty, ftz, fnx, fny, fnz, ftn, fmax,fn,
92 . h, dti
93 my_real
94 . fxn(mvsiz),fyn(mvsiz),fzn(mvsiz),stf(mvsiz),
95 . fxt(mvsiz),fyt(mvsiz),fzt(mvsiz),
96 . l1(mvsiz),l2(mvsiz),l3(mvsiz),
97 . vx1(mvsiz),vy1(mvsiz),vz1(mvsiz),
98 . vx2(mvsiz),vy2(mvsiz),vz2(mvsiz),
99 . vx3(mvsiz),vy3(mvsiz),vz3(mvsiz),
100 . vxh(mvsiz),vyh(mvsiz),vzh(mvsiz)
101C-----------------------------------------------
102 adrbuf=igrsurf(ksurf)%IAD_BUFR
103 a =bufsf(adrbuf+1)
104 b =bufsf(adrbuf+2)
105 c =bufsf(adrbuf+3)
106 DO i=1,9
107 rot(i)=bufsf(adrbuf+7+i-1)
108 END DO
109C-----------------------------------------------
110C Noeud main dans le repere de l'ellipsoide :
111C-----------------------------------------------
112 x1=bufsf(adrbuf+16)-bufsf(adrbuf+4)
113 x2=bufsf(adrbuf+17)-bufsf(adrbuf+5)
114 x3=bufsf(adrbuf+18)-bufsf(adrbuf+6)
115 xm=rot(1)*x1+rot(2)*x2+rot(3)*x3
116 ym=rot(4)*x1+rot(5)*x2+rot(6)*x3
117 zm=rot(7)*x1+rot(8)*x2+rot(9)*x3
118C-----------------------------------------------
119C Vitesse du noeud main en local.
120C-----------------------------------------------
121 v1=bufsf(adrbuf+19)
122 v2=bufsf(adrbuf+20)
123 v3=bufsf(adrbuf+21)
124 vxm=rot(1)*v1+rot(2)*v2+rot(3)*v3
125 vym=rot(4)*v1+rot(5)*v2+rot(6)*v3
126 vzm=rot(7)*v1+rot(8)*v2+rot(9)*v3
127C-----------------------------------------------
128C Vitesse de rotation du noeud main en local.
129C-----------------------------------------------
130 v1=bufsf(adrbuf+22)
131 v2=bufsf(adrbuf+23)
132 v3=bufsf(adrbuf+24)
133 vrx=rot(1)*v1+rot(2)*v2+rot(3)*v3
134 vry=rot(4)*v1+rot(5)*v2+rot(6)*v3
135 vrz=rot(7)*v1+rot(8)*v2+rot(9)*v3
136C----------------------------------------------------------
137C PASSAGE DES VITESSES EN LOCAL.
138C non optimise (plusieurs fois / noeud).
139C----------------------------------------------------------
140 DO 100 nls=ndeb+1,min(ndeb+mvsiz,ntc)
141 il =ktc(nls)
142 IF (iactiv(1,il)<=0) GOTO 100
143 i =nls-ndeb
144C------
145 in1=ksi(1,il)
146 v1=v(1,in1)
147 v2=v(2,in1)
148 v3=v(3,in1)
149 vx1(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
150 vy1(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
151 vz1(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
152C------
153 in2=ksi(2,il)
154 v1=v(1,in2)
155 v2=v(2,in2)
156 v3=v(3,in2)
157 vx2(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
158 vy2(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
159 vz2(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
160C------
161 in3=ksi(3,il)
162 v1=v(1,in3)
163 v2=v(2,in3)
164 v3=v(3,in3)
165 vx3(i)=rot(1)*v1+rot(2)*v2+rot(3)*v3
166 vy3(i)=rot(4)*v1+rot(5)*v2+rot(6)*v3
167 vz3(i)=rot(7)*v1+rot(8)*v2+rot(9)*v3
168 100 CONTINUE
169C----------------------------------------------------------
170C Force normale.
171C----------------------------------------------------------
172 DO 200 nls=ndeb+1,min(ndeb+mvsiz,ntc)
173 il =ktc(nls)
174 IF (iactiv(1,il)<=0) GOTO 200
175 i =nls-ndeb
176C STF(I)=STFAC*DEPTH(I)/MAX(DEPTH(I)-PENET(I),EM20)
177 stf(i)=stfac*depth(i)**2/max((depth(i)-penet(i))**2,em20)
178 fxn(i)=stf(i)*penet(i)*nxi(i)
179 fyn(i)=stf(i)*penet(i)*nyi(i)
180 fzn(i)=stf(i)*penet(i)*nzi(i)
181 200 CONTINUE
182C-------------------------------
183C Extraction du pt de contact au cycle precedent.
184C-------------------------------
185 DO 250 nls=ndeb+1,min(ndeb+mvsiz,ntc)
186 il =ktc(nls)
187 IF (iactiv(1,il)<=0) GOTO 250
188 i =nls-ndeb
189C Point sur le triangle tq D(P,L) etait maximum au cycle precedent.
190 l1(i)=iold(1,4*(il-1)+1)
191 l2(i)=iold(2,4*(il-1)+1)
192 l3(i)=iold(3,4*(il-1)+1)
193 250 CONTINUE
194C----------------------------------------------------------
195C Frottement :
196C----------------------------------------------------------
197#include "vectorize.inc"
198 DO 275 nls=ndeb+1,min(ndeb+mvsiz,ntc)
199 il =ktc(nls)
200 IF (iactiv(1,il)<=0) GOTO 275
201 i =nls-ndeb
202C------------------------------
203 fxt(i)=zero
204 fyt(i)=zero
205 fzt(i)=zero
206C-------------------------------
207C TOURNER DELTA(T-1) DANS LE PLAN TANGENT.
208C-------------------------------
209 IF (iactiv(1,il)>2) THEN
210 deltx=dold(1,4*(il-1)+1)
211 delty=dold(2,4*(il-1)+1)
212 deltz=dold(3,4*(il-1)+1)
213 nd =sqrt(deltx*deltx+delty*delty+deltz*deltz)
214 IF (nd/=zero) THEN
215 pn=deltx*nxi(i)+delty*nyi(i)+deltz*nzi(i)
216 deltx=deltx-pn*nxi(i)
217 delty=delty-pn*nyi(i)
218 deltz=deltz-pn*nzi(i)
219 scale=nd/sqrt(deltx*deltx+delty*delty+deltz*deltz)
220 deltx=scale*deltx
221 delty=scale*delty
222 deltz=scale*deltz
223 ENDIF
224 ELSE
225 deltx=zero
226 delty=zero
227 deltz=zero
228 ENDIF
229C-------------------------------
230C INCREMENT SUR DELTA (T-1 -> T) : VITESSE RELATIVE TANGENTE DE POLD * DT
231C-------------------------------
232 xh=hold(1,4*(il-1)+1)
233 yh=hold(2,4*(il-1)+1)
234 zh=hold(3,4*(il-1)+1)
235 v1x2=vrx*(yh-ym)
236 v2x1=vry*(xh-xm)
237 v2x3=vry*(zh-zm)
238 v3x2=vrz*(yh-ym)
239 v3x1=vrz*(xh-xm)
240 v1x3=vrx*(zh-zm)
241 v3 =v1x2-v2x1
242 v1 =v2x3-v3x2
243 v2 =v3x1-v1x3
244 vxh(i)=vxm+v1
245 vyh(i)=vym+v2
246 vzh(i)=vzm+v3
247 IF (iactiv(1,il)>=2) THEN
248 vxk=l1(i)*vx1(i)+l2(i)*vx2(i)+l3(i)*vx3(i)
249 vyk=l1(i)*vy1(i)+l2(i)*vy2(i)+l3(i)*vy3(i)
250 vzk=l1(i)*vz1(i)+l2(i)*vz2(i)+l3(i)*vz3(i)
251 dvx=vxh(i)-vxk
252 dvy=vyh(i)-vyk
253 dvz=vzh(i)-vzk
254 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
255 dvx=dvx-pn*nxi(i)
256 dvy=dvy-pn*nyi(i)
257 dvz=dvz-pn*nzi(i)
258 ELSE
259 dvx=zero
260 dvy=zero
261 dvz=zero
262 ENDIF
263C-----
264C Force tangente K.DELTA, K=1. et Norme.
265 dold(1,4*(il-1)+1)=deltx+dvx*dt1
266 dold(2,4*(il-1)+1)=delty+dvy*dt1
267 dold(3,4*(il-1)+1)=deltz+dvz*dt1
268 fxt(i)=stf(i)*dold(1,4*(il-1)+1)
269 fyt(i)=stf(i)*dold(2,4*(il-1)+1)
270 fzt(i)=stf(i)*dold(3,4*(il-1)+1)
271C-----
272 275 CONTINUE
273C------------------------------------------------------------
274C FRICTION (LOCAL COULOMB FRICTION).
275C------------------------------------------------------------
276#include "vectorize.inc"
277 DO 625 nls=ndeb+1,min(ndeb+mvsiz,ntc)
278 il =ktc(nls)
279 IF (iactiv(1,il)<=0) GOTO 625
280 i =nls-ndeb
281C-----
282 ftx=fxt(i)
283 fty=fyt(i)
284 ftz=fzt(i)
285 fnx=fxn(i)
286 fny=fyn(i)
287 fnz=fzn(i)
288 ftn=sqrt(ftx*ftx+fty*fty+ftz*ftz)
289 fn =sqrt(fnx*fnx+fny*fny+fnz*fnz)
290 fmax= max(vfric*fn,zero)
291 IF (ftn>fmax) THEN
292 scale =fmax/ftn
293 ftx=scale*ftx
294 fty=scale*fty
295 ftz=scale*ftz
296 ELSE
297 scale=1.
298 ENDIF
299 IF (iactiv(1,il)>1) THEN
300C Glissement DELTA=scale.DELTA
301 deltx=scale*dold(1,4*(il-1)+1)
302 delty=scale*dold(2,4*(il-1)+1)
303 deltz=scale*dold(3,4*(il-1)+1)
304C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
305 dold(1,4*(il-1)+1)=deltx
306 dold(2,4*(il-1)+1)=delty
307 dold(3,4*(il-1)+1)=deltz
308 ENDIF
309 fxt(i)=ftx
310 fyt(i)=fty
311 fzt(i)=ftz
312 625 CONTINUE
313C------------------------------------------------------------
314C LOCAL COORDINATES OF CONTACT POINT WITH RESPECT TO ELEMENT.
315C------------------------------------------------------------
316C#include "vectorize.inc"
317 DO 650 nls=ndeb+1,min(ndeb+mvsiz,ntc)
318 il =ktc(nls)
319 IF (iactiv(1,il)<=0) GOTO 650
320 i =nls-ndeb
321C-----
322 dx2=xp2(1,i)-xtk(i)
323 dy2=xp2(2,i)-ytk(i)
324 dz2=xp2(3,i)-ztk(i)
325 dx3=xp3(1,i)-xtk(i)
326 dy3=xp3(2,i)-ytk(i)
327 dz3=xp3(3,i)-ztk(i)
328 dx=dy2*dz3-dz2*dy3
329 dy=dx3*dz2-dz3*dx2
330 dz=dx2*dy3-dy2*dx3
331 s1=half*sqrt(dx*dx+dy*dy+dz*dz)
332C-----
333 dx1=xp1(1,i) -xtk(i)
334 dy1=xp1(2,i) -ytk(i)
335 dz1=xp1(3,i) -ztk(i)
336 dx=dy1*dz3-dz1*dy3
337 dy=dx3*dz1-dz3*dx1
338 dz=dx1*dy3-dy1*dx3
339 s2=half*sqrt(dx*dx+dy*dy+dz*dz)
340C-----
341 dx=dy1*dz2-dz1*dy2
342 dy=dx2*dz1-dz2*dx1
343 dz=dx1*dy2-dy1*dx2
344 s3=half*sqrt(dx*dx+dy*dy+dz*dz)
345C-----
346 s=one/(s1+s2+s3)
347C-----
348 iold(1,4*(il-1)+1)=s1*s
349 iold(2,4*(il-1)+1)=s2*s
350 iold(3,4*(il-1)+1)=s3*s
351 650 CONTINUE
352C------------------------------------------------------------
353C MISE EN MEMOIRE BUFFERS D'INTERFACE (-> PROCHAIN CYCLE).
354C------------------------------------------------------------
355#include "vectorize.inc"
356 DO 700 nls=ndeb+1,min(ndeb+mvsiz,ntc)
357 il =ktc(nls)
358 IF (iactiv(1,il)<=0) GOTO 700
359 i =nls-ndeb
360C-----
361 hold(1,4*(il-1)+1)=xi(i)
362 hold(2,4*(il-1)+1)=yi(i)
363 hold(3,4*(il-1)+1)=zi(i)
364C-----
365 nold(1,4*(il-1)+1)=nxi(i)
366 nold(2,4*(il-1)+1)=nyi(i)
367 nold(3,4*(il-1)+1)=nzi(i)
368 700 CONTINUE
369C------------------------------------------------------------
370C MISE EN MEMOIRE (VECTEURS DE TRAVAIL -> ASSEMBLAGE).
371C------------------------------------------------------------
372C non vectoriel.
373 DO 600 nls=ndeb+1,min(ndeb+mvsiz,ntc)
374 il =ktc(nls)
375 IF (iactiv(1,il)<=0) GOTO 600
376 i =nls-ndeb
377C-----
378 in1=ksi(1,il)
379 in2=ksi(2,il)
380 in3=ksi(3,il)
381C-----
382 wns(in1) =wns(in1)+iold(1,4*(il-1)+1)*stf(i)
383 wns(in2) =wns(in2)+iold(2,4*(il-1)+1)*stf(i)
384 wns(in3) =wns(in3)+iold(3,4*(il-1)+1)*stf(i)
385C-----
386 600 CONTINUE
387C------------------------------------------------------------
388C non vectoriel.
389 DO 800 nls=ndeb+1,min(ndeb+mvsiz,ntc)
390 il =ktc(nls)
391 IF (iactiv(1,il)<=0) GOTO 800
392 i =nls-ndeb
393C-----
394 in1=ksi(1,il)
395 in2=ksi(2,il)
396 in3=ksi(3,il)
397C-----
398 f1=fxn(i)
399 fnormx=fnormx+f1
400 f11=iold(1,4*(il-1)+1)*f1
401 f12=iold(2,4*(il-1)+1)*f1
402 f13=iold(3,4*(il-1)+1)*f1
403 wnf(1,in1)=wnf(1,in1)+f11
404 wnf(1,in2)=wnf(1,in2)+f12
405 wnf(1,in3)=wnf(1,in3)+f13
406C-----
407 f1=fyn(i)
408 fnormy=fnormy+f1
409 f11=iold(1,4*(il-1)+1)*f1
410 f12=iold(2,4*(il-1)+1)*f1
411 f13=iold(3,4*(il-1)+1)*f1
412 wnf(2,in1)=wnf(2,in1)+f11
413 wnf(2,in2)=wnf(2,in2)+f12
414 wnf(2,in3)=wnf(2,in3)+f13
415C-----
416 f1=fzn(i)
417 fnormz=fnormz+f1
418 f11=iold(1,4*(il-1)+1)*f1
419 f12=iold(2,4*(il-1)+1)*f1
420 f13=iold(3,4*(il-1)+1)*f1
421 wnf(3,in1)=wnf(3,in1)+f11
422 wnf(3,in2)=wnf(3,in2)+f12
423 wnf(3,in3)=wnf(3,in3)+f13
424C-----
425 800 CONTINUE
426C------------------------------------------------------------
427C non vectoriel.
428 DO 825 nls=ndeb+1,min(ndeb+mvsiz,ntc)
429 il =ktc(nls)
430 IF (iactiv(1,il)<=0) GOTO 825
431 i =nls-ndeb
432C-----
433 in1=ksi(1,il)
434 in2=ksi(2,il)
435 in3=ksi(3,il)
436C-----
437 f1=fxt(i)
438 ftangx=ftangx+f1
439 f11=iold(1,4*(il-1)+1)*f1
440 f12=iold(2,4*(il-1)+1)*f1
441 f13=iold(3,4*(il-1)+1)*f1
442 wtf(1,in1)=wtf(1,in1)+f11
443 wtf(1,in2)=wtf(1,in2)+f12
444 wtf(1,in3)=wtf(1,in3)+f13
445C-----
446 f1=fyt(i)
447 ftangy=ftangy+f1
448 f11=iold(1,4*(il-1)+1)*f1
449 f12=iold(2,4*(il-1)+1)*f1
450 f13=iold(3,4*(il-1)+1)*f1
451 wtf(2,in1)=wtf(2,in1)+f11
452 wtf(2,in2)=wtf(2,in2)+f12
453 wtf(2,in3)=wtf(2,in3)+f13
454C-----
455 f1=fzt(i)
456 ftangz=ftangz+f1
457 f11=iold(1,4*(il-1)+1)*f1
458 f12=iold(2,4*(il-1)+1)*f1
459 f13=iold(3,4*(il-1)+1)*f1
460 wtf(3,in1)=wtf(3,in1)+f11
461 wtf(3,in2)=wtf(3,in2)+f12
462 wtf(3,in3)=wtf(3,in3)+f13
463C-----
464 825 CONTINUE
465C------------------------------------------------------------
466C KINEMATIC TIME STEP.
467C------------------------------------------------------------
468 dti=ep20
469 DO 900 nls=ndeb+1,min(ndeb+mvsiz,ntc)
470 il =ktc(nls)
471 IF (iactiv(1,il)<=0) GOTO 900
472 i =nls-ndeb
473 dvx=vx1(i)-vxh(i)
474 dvy=vy1(i)-vyh(i)
475 dvz=vz1(i)-vzh(i)
476 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
477 IF (pn<zero) THEN
478 h =xp1(1,i)*nxi(i)+xp1(2,i)*nyi(i)+xp1(3,i)*nzi(i)
479 dti=min(dti,half*h/max(em20,-pn))
480 ENDIF
481 dvx=vx2(i)-vxh(i)
482 dvy=vy2(i)-vyh(i)
483 dvz=vz2(i)-vzh(i)
484 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
485 IF (pn<zero) THEN
486 h =xp2(1,i)*nxi(i)+xp2(2,i)*nyi(i)+xp2(3,i)*nzi(i)
487 dti=min(dti,half*h/max(em20,-pn))
488 ENDIF
489 dvx=vx3(i)-vxh(i)
490 dvy=vy3(i)-vyh(i)
491 dvz=vz3(i)-vzh(i)
492 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
493 IF (pn<zero) THEN
494 h =xp3(1,i)*nxi(i)+xp3(2,i)*nyi(i)+xp3(3,i)*nzi(i)
495 dti=min(dti,half*h/max(em20,-pn))
496 ENDIF
497 IF(dti<dt2t)THEN
498 dt2t = dti
499 neltst = noint
500 ityptst = 10
501 ENDIF
502C-----
503 900 CONTINUE
504C----------------------------------
505 IF (dti<=0.) THEN
506 dti=ep20
507 DO 950 nls=ndeb+1,min(ndeb+mvsiz,ntc)
508 il =ktc(nls)
509 IF (iactiv(1,il)<=0) GOTO 950
510 i =nls-ndeb
511 dvx=vx1(i)-vxh(i)
512 dvy=vy1(i)-vyh(i)
513 dvz=vz1(i)-vzh(i)
514 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
515 IF (pn<zero) THEN
516 h =xp1(1,i)*nxi(i)+xp1(2,i)*nyi(i)+xp1(3,i)*nzi(i)
517 dti=min(dti,half*h/max(em20,-pn))
518 ENDIF
519 dvx=vx2(i)-vxh(i)
520 dvy=vy2(i)-vyh(i)
521 dvz=vz2(i)-vzh(i)
522 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
523 IF (pn<zero) THEN
524 h =xp2(1,i)*nxi(i)+xp2(2,i)*nyi(i)+xp2(3,i)*nzi(i)
525 dti=min(dti,half*h/max(em20,-pn))
526 ENDIF
527 dvx=vx3(i)-vxh(i)
528 dvy=vy3(i)-vyh(i)
529 dvz=vz3(i)-vzh(i)
530 pn=dvx*nxi(i)+dvy*nyi(i)+dvz*nzi(i)
531 IF (pn<zero) THEN
532 h =xp3(1,i)*nxi(i)+xp3(2,i)*nyi(i)+xp3(3,i)*nzi(i)
533 dti=min(dti,half*h/max(em20,-pn))
534 ENDIF
535 IF(dti<=zero)THEN
536 in1=ksi(1,il)
537 in2=ksi(2,il)
538 in3=ksi(3,il)
539#include "lockon.inc"
540 WRITE(istdo,'(A,E12.4,A,I8)')
541 . ' **WARNING MINIMUM TIME STEP ',dti,' IN INTERFACE ',noint
542 WRITE(iout ,'(A,E12.4,A,I8)')
543 . ' **WARNING MINIMUM TIME STEP ',dti,' IN INTERFACE ',noint
544 WRITE(iout,'(a,3i8)') ' element/segment nodes :',
545 . IN1,IN2,IN3
546#include "lockoff.inc"
547 TSTOP = TT
548 ENDIF
549C-----
550 950 CONTINUE
551 ENDIF
552C----------------------------------
553 RETURN
554 END
i15fort1
subroutine i15fort1(ndeb, ntc, stfac, x, v, ksurf, igrsurf, bufsf, ktc, ksi, iactiv, iold, hold, nold, dold, xp1, xp2, xp3, 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 i15fort1.F:40
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