OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i24dst3e.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!|| i24dst3e ../engine/source/interfaces/int24/i24dst3e.F
25!||--- called by ------------------------------------------------------
26!|| i24mainf ../engine/source/interfaces/int24/i24main.F
27!||--- calls -----------------------------------------------------
28!|| bitget ../engine/source/interfaces/intsort/i20sto.F
29!||--- uses -----------------------------------------------------
30!|| h3d_mod ../engine/share/modules/h3d_mod.F
31!|| tri7box ../engine/share/modules/tri7box.F
32!||====================================================================
33 SUBROUTINE i24dst3e(
34 1 JLT ,A ,X ,CAND_N ,CAND_E ,
35 2 MBINFLG,ISEADD ,ISEDGE ,NSVG ,NIN ,
36 3 IXX ,STIF ,
37 4 JLT_NEW,INACTI ,XI ,YI ,ZI ,
38 5 XX ,YY ,ZZ ,PMAX_GAP,
39 6 FSKYI ,ISKY ,CAND_T ,FCONT ,H3D_DATA)
40C-----------------------------------------------
41C M o d u l e s
42C-----------------------------------------------
43 USE tri7box
44 USE h3d_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 D u m m y A r g u m e n t s
56C-----------------------------------------------
57 INTEGER JLT, JLT_NEW,NIN,INACTI, CAND_N(*),NSVG(MVSIZ),
58 . CAND_E(*),MBINFLG(*),ISEDGE(*),ISEADD(*),CAND_T(*)
59 INTEGER ISKY(*),IXX(MVSIZ,13)
60 my_real
61 . N1(MVSIZ), N2(MVSIZ), N3(MVSIZ), PMAX_GAP,
62 . XI(MVSIZ), YI(MVSIZ), ZI(MVSIZ), STIF(MVSIZ),
63 . fskyi(lskyi,nfskyi),x(3,*),a(3,*),fcont(3,*),
64 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17)
65 TYPE(h3d_database) :: H3D_DATA
66C-----------------------------------------------
67C C o m m o n B l o c k s
68C-----------------------------------------------
69#include "scr14_c.inc"
70#include "scr16_c.inc"
71#include "parit_c.inc"
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER I, J, K, IFSEG,IK1(4),IK2(4),IK14(4),IK5(4),ICONT,
76 . JM,JP,JG,JMG,JPG,NES,IAD,K1,K2,K5,K14,KJ,IFE,IX(4),
77 . NISKYL,M,IG
78 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ)
79 my_real
80 . XJM,YJM,ZJM,XJP,YJP,ZJP,FM(4,4),FSI(4),X1,Y1,Z1,X2,Y2,Z2,F,
81 . xijm,yijm,zijm,xijp,yijp,zijp,fx,fy,fz,
82 . x12,y12,z12,x14,y14,z14,x15,y15,z15,xij,yij,zij,aaa,bbb,
83 . sxm1,sym1,szm1,sxm2,sym2,szm2,sxm,sym,szm,xm2,xa,xb,det,
84 . sxsm,sysm,szsm,nxij,nyij,nzij,xsm,xs2m2,ys2m2,zs2m2,
85 . tx12,ty12,tz12,xj,yj,zj,dist,nxsm,nysm,nzsm,h1m,h2m,
86 . sxs1,sys1,szs1,sxs2,sys2,szs2,sxs,sys,szs,xs2,h1s,h2s,
87 . nx12,ny12,nz12,txij,tyij,tzij,f1,f2,unp,zerom,eps,fi,fj,fsj
88 INTEGER BITGET
89 EXTERNAL BITGET
90
91 DATA IK1 /1,2,3,4/
92 DATA IK2 /2,3,4,1/
93 DATA IK5 /3,1,0,2/ ! triangle only
94 DATA ik14 /14,15,16,17/
95
96 unp = one + em4
97 zerom = zero - em4
98 eps = em3
99c write(iout,*)'t=',tt,'i24dst3e 1'
100
101 DO i=1,jlt
102 IF(cand_t(i)==0)cycle ! no edge candidate
103 m=cand_e(i)
104 IF(bitget(mbinflg(m),8)==0)cycle !flag activ edge MAIN
105
106 DO j = 1,4
107 fsi(j) = zero
108 DO k = 1,4
109 fm(j,k) = zero
110 ENDDO
111 ENDDO
112
113 icont=0
114
115 IF(cand_n(i)<0)THEN
116c afaire !!!!!!!!!!!!!!!!!!!!!!!!!
117 stop 987
118 ENDIF
119 iad = iseadd(cand_n(i))
120 nes = isedge(iad)
121
122 ix1(i) = ixx(i,1)
123 ix2(i) = ixx(i,2)
124 ix3(i) = ixx(i,3)
125 ix4(i) = ixx(i,4)
126
127 k5=5
128
129 DO k = 1,4 ! edge loop on segment
130 IF(bitget(mbinflg(m),k+1)==0)cycle !flag edge MAIN
131
132 k1=ik1(k)
133 k2=ik2(k)
134 k14=ik14(k)
135 IF(ix3(i) == ix4(i)) k5=ik5(k)
136
137 x1 = xx(i,k1)
138 y1 = yy(i,k1)
139 z1 = zz(i,k1)
140
141 x2 = xx(i,k2)
142 y2 = yy(i,k2)
143 z2 = zz(i,k2)
144
145 x15 = xx(i,k5) - x1
146 y15 = yy(i,k5) - y1
147 z15 = zz(i,k5) - z1
148
149 x12 = x2 - x1
150 y12 = y2 - y1
151 z12 = z2 - z1
152
153 x14 = xx(i,k14) - x1
154 y14 = yy(i,k14) - y1
155 z14 = zz(i,k14) - z1
156
157c deux vecteurs surfaces du diedre MAIN
158 sxm1 = y12*z15 - z12*y15
159 sym1 = z12*x15 - x12*z15
160 szm1 = x12*y15 - y12*x15
161
162 sxm2 = y14*z12 - z14*y12
163 sym2 = z14*x12 - x14*z12
164 szm2 = x14*y12 - y14*x12
165
166c check if convex
167
168 aaa = (sym1*szm2 - szm1*sym2)*x12
169 . + (szm1*sxm2 - sxm1*szm2)*y12
170 . + (sxm1*sym2 - sym1*sxm2)*z12
171
172 IF(aaa<=zero)cycle
173
174 sxm = sxm1 + sxm2
175 sym = sym1 + sym2
176 szm = szm1 + szm2
177
178 xm2 = x12*x12+y12*y12+z12*z12
179 aaa = one/sqrt(xm2)
180 tx12 = x12*aaa
181 ty12 = y12*aaa
182 tz12 = z12*aaa
183
184
185 jm=nes
186 DO kj=1,nes ! loop on SECONDARY edge
187
188 jg = isedge(iad+kj)
189 jmg = isedge(iad+kj+nes)
190 jpg = isedge(iad+kj+nes+nes)
191
192 xj=x(1,jg)
193 yj=x(2,jg)
194 zj=x(3,jg)
195
196 xjm=x(1,jmg)
197 yjm=x(2,jmg)
198 zjm=x(3,jmg)
199
200 xjp=x(1,jpg)
201 yjp=x(2,jpg)
202 zjp=x(3,jpg)
203
204 xij = xj - xi(i)
205 yij = yj - yi(i)
206 zij = zj - zi(i)
207
208 xijm = xjm - xi(i)
209 yijm = yjm - yi(i)
210 zijm = zjm - zi(i)
211
212 xijp = xjp - xi(i)
213 yijp = yjp - yi(i)
214 zijp = zjp - zi(i)
215
216c two area vectors on SECONDARY dihedral
217 sxs1 = yij*zijm - zij*yijm
218 sys1 = zij*xijm - xij*zijm
219 szs1 = xij*yijm - yij*xijm
220
221 sxs2 = yijp*zij - zijp*yij
222 sys2 = zijp*xij - xijp*zij
223 szs2 = xijp*yij - yijp*xij
224
225c check if convex
226
227 aaa = (sys1*szs2 - szs1*sys2)*xij
228 . + (szs1*sxs2 - sxs1*szs2)*yij
229 . + (sxs1*sys2 - sys1*sxs2)*zij
230
231 IF(aaa<=zero)cycle
232
233 sxs = sxs1 + sxs2
234 sys = sys1 + sys2
235 szs = szs1 + szs2
236
237 aaa = sxs*sxm+sys*sym+szs*szm
238 IF(aaa > zero) cycle ! MAIN and SECONDARY dihedral are not locking each other
239
240c vecteur SECONDARY MAIN
241 sxsm = y12*zij - z12*yij
242 sysm = z12*xij - x12*zij
243 szsm = x12*yij - y12*xij
244
245 aaa = sxs*sxsm+sys*sysm+szs*szsm
246 IF(aaa < zero) THEN
247 sxsm = - sxsm
248 sysm = - sysm
249 szsm = - szsm
250 ENDIF
251 aaa = sxm*sxsm+sym*sysm+szm*szsm
252 IF(aaa > zero) cycle ! inconsistent
253
254 aaa = sxsm*(xi(i)-x1)
255 . + sysm*(yi(i)-y1)
256 . + szsm*(zi(i)-z1)
257 aaa = -aaa
258 bbb = sxsm*sxsm + sysm*sysm + szsm*szsm
259
260 IF(aaa >= zero)cycle ! no intersection
261 IF(bbb == zero)cycle ! // edges
262
263 bbb = one/sqrt(bbb)
264 dist = aaa*bbb
265
266 nxsm = sxsm*bbb
267 nysm = sysm*bbb
268 nzsm = szsm*bbb
269
270 xs2 = xij*xij+yij*yij+zij*zij
271
272 aaa = one/sqrt(nxsm*nxsm+nysm*nysm+nzsm*nzsm)
273 nxsm = nxsm*aaa
274 nysm = nysm*aaa
275 nzsm = nzsm*aaa
276
277 goto 123
278
279C=======================================================================
280c check if normal is inside MAIN dihedra
281
282 aaa = tx12*nxsm + ty12*nysm + tz12*nzsm
283 nx12 = nxsm*aaa
284 ny12 = nysm*aaa
285 nz12 = nzsm*aaa
286
287 nxsm = nxsm - nx12 ! 12 component is removed
288 nysm = nysm - ny12 ! and resored
289 nzsm = nzsm - nz12 ! later
290
291
292 aaa = (sym1*nzsm - szm1*nysm)*tx12
293 . + (szm1*nxsm - sxm1*nzsm)*ty12
294 . + (sxm1*nysm - sym1*nxsm)*tz12
295
296 IF(aaa<=zero)THEN
297c project on SM1
298 aaa = sxm1*sxm1 + sym1*sym1 + szm1*szm1
299 aaa = (sxm1*nxsm + sym1*nysm + szm1*nzsm)/sqrt(aaa)
300 nxsm = sxm1*aaa
301 nysm = sym1*aaa
302 nzsm = szm1*aaa
303 ENDIF
304
305 aaa = (nysm*szm2 - nzsm*sym2)*tx12
306 . + (nzsm*sxm2 - nxsm*szm2)*ty12
307 . + (nxsm*sym2 - nysm*sxm2)*tz12
308
309 IF(aaa<=zero)THEN
310c project on SM2
311 aaa = sxm2*sxm2 + sym2*sym2 + szm2*szm2
312 aaa = (sxm2*nxsm + sym2*nysm + szm2*nzsm)/sqrt(aaa)
313 nxsm = sxm2*aaa
314 nysm = sym2*aaa
315 nzsm = szm2*aaa
316 ENDIF
317
318 nxsm = -(nxsm + nx12) ! restore
319 nysm = -(nysm + ny12) ! and
320 nzsm = -(nzsm + nz12) ! invert sign
321
322c check if normal is inside SECONDARY dihedra
323
324 aaa = one/sqrt(xs2)
325 txij = xij*aaa
326 tyij = yij*aaa
327 tzij = zij*aaa
328
329 aaa = txij*nxsm + tyij*nysm + tzij*nzsm
330 nxij = nxsm*aaa
331 nyij = nysm*aaa
332 nzij = nzsm*aaa
333
334 nxsm = nxsm - nxij ! IJ component is removed
335 nysm = nysm - nyij ! and resored
336 nzsm = nzsm - nzij ! later
337
338
339 aaa = (sys1*nzsm - szs1*nysm)*txij
340 . + (szs1*nxsm - sxs1*nzsm)*tyij
341 . + (sxs1*nysm - sys1*nxsm)*tzij
342
343 IF(aaa<=zero)THEN
344c project on SM1
345 aaa = sxs1*sxs1 + sys1*sys1 + szs1*szs1
346 aaa = (sxs1*nxsm + sys1*nysm + szs1*nzsm)/sqrt(aaa)
347 nxsm = sxs1*aaa
348 nysm = sys1*aaa
349 nzsm = szs1*aaa
350 ENDIF
351
352 aaa = (nysm*szs2 - nzsm*sys2)*txij
353 . + (nzsm*sxs2 - nxsm*szs2)*tyij
354 . + (nxsm*sys2 - nysm*sxs2)*tzij
355
356 IF(aaa<=zero)THEN
357c project on SM2
358 aaa = sxs2*sxs2 + sys2*sys2 + szs2*szs2
359 aaa = (sxs2*nxsm + sys2*nysm + szs2*nzsm)/sqrt(aaa)
360 nxsm = sxs2*aaa
361 nysm = sys2*aaa
362 nzsm = szs2*aaa
363 ENDIF
364
365 nxsm = -(nxsm + nxij) ! restore
366 nysm = -(nysm + nyij) ! and
367 nzsm = -(nzsm + nzij) ! invert sign
368
369C=======================================================================
370
371 123 continue
372
373
374c compute relative location on MAIN and SECONDARY edge
375c same as int 11 (i11dst3.F)
376 xsm = xij*x12 + yij*y12 + zij*z12
377
378 xs2m2 = x2-xj
379 ys2m2 = y2-yj
380 zs2m2 = z2-zj
381
382 xa = xij*xs2m2 + yij*ys2m2 + zij*zs2m2
383 xb = -x12*xs2m2 - y12*ys2m2 - z12*zs2m2
384 det = xm2*xs2 - xsm*xsm
385 det = max(em20,det)
386C
387 h1m = (xa*xsm-xb*xs2) / det
388
389 IF(h1m < -em3)cycle ! no contact
390 IF(h1m > one+em3)cycle ! no contact
391C
392 xs2 = max(xs2,em20)
393 xm2 = max(xm2,em20)
394 h1m = min(one,max(zero,h1m ))
395 h1s = -(xa + h1m *xsm) / xs2
396
397 IF(h1s > one+em3)cycle ! no contact
398 IF(h1s < half-em3)cycle ! edge IJ will be solved as JI
399
400 h1s = min(one,max(zero,h1s))
401
402 h1m = -(xb + h1s *xsm) / xm2
403 h1m = min(one,max(zero,h1m ))
404C
405 h2s = one -h1s
406 h2m = one -h1m
407
408 icont=1
409c compute forces
410 f = dist*stif(i)
411 IF(h1s < half+em3)THEN
412 aaa=(h1s-half+em3)/(two*em3)
413 f = f*aaa
414 ELSEIF(h1s > one-em3)THEN
415 aaa=-(h1s-one-em3)/(two*em3)
416 f = f*aaa
417 ENDIF
418 IF(h1m < em3)THEN
419 aaa=(h1m+em3)/(two*em3)
420 f = f*aaa
421 ELSEIF(h1m > one-em3)THEN
422 aaa=-(h1m-one-em3)/(two*em3)
423 f = f*aaa
424 ENDIF
425 f1=-h1m*f
426 f2=-h2m*f
427 fi=h1s*f
428 fj=h2s*f
429
430 fm(1,k1) = fm(1,k1) + f1*nxsm
431 fm(2,k1) = fm(2,k1) + f1*nysm
432 fm(3,k1) = fm(3,k1) + f1*nzsm
433 fm(4,k1) = fm(4,k1) + stif(i)*abs(h1m)
434 fm(1,k2) = fm(1,k2) + f2*nxsm
435 fm(2,k2) = fm(2,k2) + f2*nysm
436 fm(3,k2) = fm(3,k2) + f2*nzsm
437 fm(4,k2) = fm(4,k2) + stif(i)*abs(h2m)
438 fsi(1) = fsi(1) + fi*nxsm
439 fsi(2) = fsi(2) + fi*nysm
440 fsi(3) = fsi(3) + fi*nzsm
441 fsi(4) = fsi(4) + stif(i)*abs(h1s)
442c save FJ
443 IF(fj/=zero)THEN
444#include "lockon.inc"
445 nisky = nisky + 1
446 niskyl = nisky
447#include "lockoff.inc"
448 fx = fj*nxsm
449 fy = fj*nysm
450 fz = fj*nzsm
451 fskyi(niskyl,1) = fx
452 fskyi(niskyl,2) = fy
453 fskyi(niskyl,3) = fz
454 fskyi(niskyl,4)=stif(i)*abs(h2s)
455 isky(niskyl) = jg
456 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT > 0)THEN
457 fcont(1,jg) =fcont(1,jg) + fx
458 fcont(2,jg) =fcont(2,jg) + fy
459 fcont(3,jg) =fcont(3,jg) + fz
460 ENDIF
461 ENDIF
462
463 ENDDO ! J=1,NES ! loop on SECONDARY edge
464 ENDDO ! K = 1,4 ! edge loop on segment
465c save F(Ki) and FSI
466 IF(icont==0)cycle
467
468 ix(1) = ix1(i)
469 ix(2) = ix2(i)
470 ix(3) = ix3(i)
471 ix(4) = ix4(i)
472 DO k=1,4
473 IF(fm(4,k)/=zero)THEN
474#include "lockon.inc"
475 nisky = nisky + 1
476 niskyl = nisky
477#include "lockoff.inc"
478 fskyi(niskyl,1)=fm(1,k)
479 fskyi(niskyl,2)=fm(2,k)
480 fskyi(niskyl,3)=fm(3,k)
481 fskyi(niskyl,4)=fm(4,k)
482 isky(niskyl) = ix(k)
483 ENDIF
484 ENDDO
485 IF(fsi(4)/=zero)THEN
486#include "lockon.inc"
487 nisky = nisky + 1
488 niskyl = nisky
489#include "lockoff.inc"
490 fskyi(niskyl,1)=fsi(1)
491 fskyi(niskyl,2)=fsi(2)
492 fskyi(niskyl,3)=fsi(3)
493 fskyi(niskyl,4)=fsi(4)
494 isky(niskyl) = nsvg(i)
495 ENDIF
496
497 IF(anim_v(4)+outp_v(4)+h3d_data%N_VECT_CONT > 0)THEN
498 DO k=1,4
499 IF(fm(4,k)/=zero)THEN
500 ig=ix(k)
501 fcont(1,ig) = fcont(1,ig) + fm(1,k)
502 fcont(2,ig) = fcont(2,ig) + fm(2,k)
503 fcont(3,ig) = fcont(3,ig) + fm(3,k)
504 ENDIF
505 ENDDO
506 IF(fsi(4)/=zero)THEN
507 jg=nsvg(i)
508 fcont(1,jg) = fcont(1,jg) + fsi(1)
509 fcont(2,jg) = fcont(2,jg) + fsi(2)
510 fcont(3,jg) = fcont(3,jg) + fsi(3)
511 ENDIF
512 ENDIF
513
514
515
516 ENDDO ! I=1,JLT
517
518c write(iout,*)'i24dst3e 10'
519
520
521 RETURN
522 END
i24dst3e
subroutine i24dst3e(jlt, a, x, cand_n, cand_e, mbinflg, iseadd, isedge, nsvg, nin, ixx, stif, jlt_new, inacti, xi, yi, zi, xx, yy, zz, pmax_gap, fskyi, isky, cand_t, fcont, h3d_data)
Definition i24dst3e.F:40
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
h3d_mod
Definition h3d_mod.F:308
h3d_mod::h3d_database
Definition h3d_mod.F:435