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

Go to the source code of this file.

Functions/Subroutines

subroutine i25dst3_21 (jlt, cand_n, cand_e, nrtm, xx, yy, zz, xi, yi, zi, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, mseglo, gaps, gapm, irect, irtlm, time_s, gap_nm, itab, nnx, nny, nnz, far, pent, dist, lb, lc, lbp, lcp, kslide, mvoisn, gapmxl, ibound, vtx_bisector, etyp, icodt, iskew, drad, dgapload)

Function/Subroutine Documentation

◆ i25dst3_21()

subroutine i25dst3_21 ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
integer nrtm,
xx,
yy,
zz,
xi,
yi,
zi,
integer nin,
integer nsn,
integer, dimension(mvsiz) ix1,
integer, dimension(mvsiz) ix2,
integer, dimension(mvsiz) ix3,
integer, dimension(mvsiz) ix4,
integer, dimension(mvsiz) nsvg,
stif,
integer inacti,
integer, dimension(*) mseglo,
gaps,
gapm,
integer, dimension(4,*) irect,
integer, dimension(4,nsn) irtlm,
time_s,
gap_nm,
integer, dimension(*) itab,
nnx,
nny,
nnz,
integer, dimension(mvsiz,4) far,
pent,
dist,
lb,
lc,
lbp,
lcp,
integer, dimension(mvsiz,4) kslide,
integer, dimension(mvsiz,4) mvoisn,
gapmxl,
integer, dimension(4,mvsiz) ibound,
real*4, dimension(3,2,*) vtx_bisector,
integer, dimension(*) etyp,
integer, dimension(*) icodt,
integer, dimension(*) iskew,
intent(in) drad,
intent(in) dgapload )

Definition at line 30 of file i25dst3_21.F.

41C-----------------------------------------------
42C M o d u l e s
43C-----------------------------------------------
44 USE tri7box
45C-----------------------------------------------
46C I m p l i c i t T y p e s
47C-----------------------------------------------
48#include "implicit_f.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"
53#include "task_c.inc"
54#include "comlock.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER JLT, NIN, NSN, INACTI, NRTM,
59 . CAND_N(*),
60 . CAND_E(*),NSVG(MVSIZ), ETYP(*), ICODT(*), ISKEW(*)
61 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
62 . MSEGLO(*),IRTLM(4,NSN), KSLIDE(MVSIZ,4), MVOISN(MVSIZ,4),
63 . IBOUND(4,MVSIZ)
64 my_real , INTENT(IN) :: dgapload ,drad
65 my_real
66 . time_s(2,nsn), gaps(*), gapm(*)
67 my_real
68 . xx(mvsiz,5),yy(mvsiz,5),zz(mvsiz,5),
69 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
70 . nnx(mvsiz,5), nny(mvsiz,5), nnz(mvsiz,5),
71 . pent(mvsiz,4), dist(mvsiz), lb(mvsiz,4), lc(mvsiz,4),
72 . lbp(mvsiz,4), lcp(mvsiz,4), gap_nm(4,mvsiz), gapmxl(mvsiz)
73 INTEGER IRECT(4,*),ITAB(*),FAR(MVSIZ,4)
74 real*4 vtx_bisector(3,2,*)
75C-----------------------------------------------
76C L o c a l V a r i a b l e s
77C-----------------------------------------------
78 INTEGER I, J, K, L, N, N4N, MGLOB, IB1, IB2, IB3, IBX, IX, IY, IZ
79 INTEGER I4N(MVSIZ),
80 . INGAP(MVSIZ,4),
81 . IT, JJ, I1, I2, ITRIA(2,4), SUBTRIA(MVSIZ),
82 . IBCS, ISKS, IBCM(4), ISKM(4),
83 . IBCX(MVSIZ), IBCY(MVSIZ), IBCZ(MVSIZ)
84 my_real
85 . xij(mvsiz,4), xi0v(mvsiz),
86 . yij(mvsiz,4), yi0v(mvsiz),
87 . zij(mvsiz,4), zi0v(mvsiz),
88 .
89 .
90 .
91 .
92 .
93 .
94 . pene
95 my_real
96 . gap_mm(mvsiz), gap, gap2,
97 . xp, yp, zp,
98 . dx, dy, dz, dmin, dd(mvsiz,4),
99 . lap, hla, hlb(mvsiz,4), hlc(mvsiz,4),
100 . al(mvsiz,4)
101 my_real
102 . aaa,bb(mvsiz,4),
103 . sx1,sx2,sx3,sx4,
104 . sy1,sy2,sy3,sy4,
105 . sz1,sz2,sz3,sz4,
106 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
107 . xn(mvsiz,4),yn(mvsiz,4),zn(mvsiz,4),
108 . side(mvsiz,4),
109 . x12(mvsiz), y12(mvsiz), z12(mvsiz),
110 . px, py, pz, pp, p1, p2, xh, yh, zh, d1, d2, d3, vx, vy, vz,
111 . ll ,nn, pn, ld(mvsiz), lx, lax, lbx, lcx,
112 . prec
113 DATA itria/1,2,2,3,3,4,4,1/
114C-----------------------------------------------------------------------
115C
116C initialisation (cf triangles)
117 far(1:jlt,1:4) = 0
118 pent(1:jlt,1:4) = ep20
119 dd(1:jlt,1:4) = ep20
120 dist(1:jlt)= ep20
121 ld(1:jlt) = ep20
122C
123 DO i=1,jlt
124C
125C For computing LBP, LCP, FAR=2, etc
126C
127 IF(stif(i) <= zero)cycle
128C
129 x0n(i,1) = xx(i,1) - xx(i,5)
130 y0n(i,1) = yy(i,1) - yy(i,5)
131 z0n(i,1) = zz(i,1) - zz(i,5)
132C
133 x0n(i,2) = xx(i,2) - xx(i,5)
134 y0n(i,2) = yy(i,2) - yy(i,5)
135 z0n(i,2) = zz(i,2) - zz(i,5)
136C
137 x0n(i,3) = xx(i,3) - xx(i,5)
138 y0n(i,3) = yy(i,3) - yy(i,5)
139 z0n(i,3) = zz(i,3) - zz(i,5)
140C
141 x0n(i,4) = xx(i,4) - xx(i,5)
142 y0n(i,4) = yy(i,4) - yy(i,5)
143 z0n(i,4) = zz(i,4) - zz(i,5)
144C
145 IF(ix3(i)/=ix4(i))THEN
146 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
147 ELSE
148 gap_mm(i)=gap_nm(3,i)
149 END IF
150 ENDDO
151C--------------------------------------------------------
152#include "vectorize.inc"
153 DO i=1,jlt
154C
155 IF(stif(i) <= zero)cycle
156C
157 xi0v(i) = xx(i,5) - xi(i)
158 yi0v(i) = yy(i,5) - yi(i)
159 zi0v(i) = zz(i,5) - zi(i)
160C
161 xij(i,1) = xx(i,1) - xi(i)
162 yij(i,1) = yy(i,1) - yi(i)
163 zij(i,1) = zz(i,1) - zi(i)
164C
165 xij(i,2) = xx(i,2) - xi(i)
166 yij(i,2) = yy(i,2) - yi(i)
167 zij(i,2) = zz(i,2) - zi(i)
168C
169 xij(i,3) = xx(i,3) - xi(i)
170 yij(i,3) = yy(i,3) - yi(i)
171 zij(i,3) = zz(i,3) - zi(i)
172C
173 xij(i,4) = xx(i,4) - xi(i)
174 yij(i,4) = yy(i,4) - yi(i)
175 zij(i,4) = zz(i,4) - zi(i)
176C
177 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
178 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
179 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
180C
181 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
182 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
183 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
184C
185 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
186 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
187 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
188 nn = one/
189 . max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
190 xn(i,1)=xn(i,1)*nn
191 yn(i,1)=yn(i,1)*nn
192 zn(i,1)=zn(i,1)*nn
193 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
194 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
195C
196 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
197 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
198 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
199C
200 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
201 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
202 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
203 nn = one/
204 . max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i,2)+ zn(i,2)*zn(i,2)))
205 xn(i,2)=xn(i,2)*nn
206 yn(i,2)=yn(i,2)*nn
207 zn(i,2)=zn(i,2)*nn
208 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn(i,2)*sz3)*nn
209 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
210C
211 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
212 sy4 = zi0v(i)*xij(i,4) - xi0v(i)*zij(i,4)
213 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
214C
215 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
216 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
217 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
218 nn = one/
219 . max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
220 xn(i,3)=xn(i,3)*nn
221 yn(i,3)=yn(i,3)*nn
222 zn(i,3)=zn(i,3)*nn
223 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
224 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
225C
226 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
227 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
228 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
229 nn = one/
230 . max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn(i,4)))
231 xn(i,4)=xn(i,4)*nn
232 yn(i,4)=yn(i,4)*nn
233 zn(i,4)=zn(i,4)*nn
234 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
235 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
236C
237 END DO
238C--------------------------------------------------------
239C
240#include "vectorize.inc"
241 DO i=1,jlt
242C
243 IF(stif(i) <= zero)cycle
244C
245 aaa = one/max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
246 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
247 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
248 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
249 al(i,1) = max(zero,min(one,al(i,1)))
250 aaa = one/max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
251 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
252 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
253 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
254 al(i,2) = max(zero,min(one,al(i,2)))
255 aaa = one/max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
256 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
257 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
258 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
259 al(i,3) = max(zero,min(one,al(i,3)))
260 aaa = one/max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
261 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
262 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
263 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
264 al(i,4) = max(zero,min(one,al(i,4)))
265C
266 END DO
267C--------------------------------------------------------
268 ingap(1:jlt,1:4) = 0
269C--------------------------------------------------------
270 it=1
271 i1=itria(1,it)
272 i2=itria(2,it)
273#include "vectorize.inc"
274 DO i=1,jlt
275C
276 IF(stif(i) <= zero)cycle
277C
278 x12(i) = xx(i,i2) - xx(i,i1)
279 y12(i) = yy(i,i2) - yy(i,i1)
280 z12(i) = zz(i,i2) - zz(i,i1)
281C
282 lbp(i,it) = lb(i,it)
283 lcp(i,it) = lc(i,it)
284 lap = one-lbp(i,it)-lcp(i,it)
285C HLA, HLB, HLC required for triangle sharp angle
286 aaa = one / max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
287 hla= lap*abs(lap) * aaa
288 IF(lap<zero.AND.
289 + hla<=hlb(i,it).AND.hla<=hlc(i,it))THEN
290 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
291 lbp(i,it) = max(zero,min(one,lbp(i,it)))
292 lcp(i,it) = one - lbp(i,it)
293 ELSEIF(lbp(i,it)<zero.AND.
294 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)THEN
295 lbp(i,it) = zero
296 lcp(i,it) = al(i,i2)
297 ELSEIF(lcp(i,it)<zero.AND.
298 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))THEN
299 lcp(i,it) = zero
300 lbp(i,it) = al(i,i1)
301 ENDIF
302
303 lap = one-lbp(i,it)-lcp(i,it)
304 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
305 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
306 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
307 dx =xi(i)-xp
308 dy =yi(i)-yp
309 dz =zi(i)-zp
310 dd(i,it)=dx*dx+dy*dy+dz*dz
311
312 END DO
313C--------------------------------------------------------
314C Check if penetra vs cylindrical gap ...
315C--------------------------------------------------------
316#include "vectorize.inc"
317 DO i =1,jlt
318C
319 IF (stif(i) <= zero)cycle
320C
321 lap = one-lbp(i,it)-lcp(i,it)
322C
323 gap = min(max(drad,gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i)+dgapload),
324 . max(drad,gapmxl(i)+dgapload))
325 gap2 = gap**2
326C
327 bb(i,it) =((xx(i,5)-xi(i))*xn(i,it)+(yy(i,5)-yi(i))*yn(i,it)+(zz(i,5)-zi(i))*zn(i,it))
328C
329 IF(dd(i,it) <= gap2 .AND. bb(i,it) <= zero) ingap(i,it)=1
330
331 IF(bb(i,it) > zero)THEN
332C
333C penetration approximee (cf zone limite interpolation des normales)
334 pent(i,it)=max(zero,gap+bb(i,it))
335 ELSE
336 pent(i,it)=max(zero,gap-sqrt(dd(i,it)))
337 END IF
338C
339 ENDDO
340C--------------------------------------------------------
341 n4n=0
342 DO i=1,jlt
343C
344 IF(stif(i) <= zero)cycle
345C
346 IF(ix3(i)/=ix4(i))THEN
347 n4n = n4n+1
348 i4n(n4n)=i
349 ENDIF
350 ENDDO
351C--------------------------------------------------------
352 DO it=2,4
353
354 i1=itria(1,it)
355 i2=itria(2,it)
356
357#include "vectorize.inc"
358 DO k=1,n4n
359 i=i4n(k)
360C
361 x12(i) = xx(i,i2) - xx(i,i1)
362 y12(i) = yy(i,i2) - yy(i,i1)
363 z12(i) = zz(i,i2) - zz(i,i1)
364C
365 lbp(i,it) = lb(i,it)
366 lcp(i,it) = lc(i,it)
367 lap = one-lbp(i,it)-lcp(i,it)
368C HLA, HLB, HLC required for triangle sharp angle
369 aaa = one / max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
370 hla= lap*abs(lap) * aaa
371 IF(lap<zero.AND.
372 + hla<=hlb(i,it).AND.hla<=hlc(i,it))THEN
373 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
374 lbp(i,it) = max(zero,min(one,lbp(i,it)))
375 lcp(i,it) = one - lbp(i,it)
376 ELSEIF(lbp(i,it)<zero.AND.
377 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)THEN
378 lbp(i,it) = zero
379 lcp(i,it) = al(i,i2)
380 ELSEIF(lcp(i,it)<zero.AND.
381 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))THEN
382 lcp(i,it) = zero
383 lbp(i,it) = al(i,i1)
384 ENDIF
385
386 lap = one-lbp(i,it)-lcp(i,it)
387 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
388 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
389 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
390 dx =xi(i)-xp
391 dy =yi(i)-yp
392 dz =zi(i)-zp
393 dd(i,it)=dx*dx+dy*dy+dz*dz
394
395 END DO
396C--------------------------------------------------------
397C Check if penetra vs cylindrical gap ...
398C--------------------------------------------------------
399#include "vectorize.inc"
400 DO k=1,n4n
401 i=i4n(k)
402C
403 lap = one-lbp(i,it)-lcp(i,it)
404C
405 gap = min(max(drad,gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i)+dgapload),
406 . max(drad,gapmxl(i)+dgapload))
407 gap2 = gap**2
408C
409 bb(i,it) =((xx(i,5)-xi(i))*xn(i,it)+(yy(i,5)-yi(i))*yn(i,it)+(zz(i,5)-zi(i))*zn(i,it))
410C
411 IF(dd(i,it) <= gap2 .AND. bb(i,it) <= zero) ingap(i,it)=1
412
413 IF(bb(i,it) > zero)THEN
414C
415C penetration approximee (cf zone limite interpolation des normales)
416 pent(i,it)=max(zero,gap+bb(i,it))
417 ELSE
418 pent(i,it)=max(zero,gap-sqrt(dd(i,it)))
419 END IF
420C
421 ENDDO
422 END DO ! DO IT=2,4
423C--------------------------------------------------------
424C Look for closest sub-triangle
425C--------------------------------------------------------
426 DO i=1,jlt
427C
428 IF(stif(i) <= zero)cycle
429C
430C Erase Subtria (cf intersections)
431 subtria(i)=0
432C
433 IF(ix3(i)/=ix4(i))THEN
434 dmin=min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
435
436 DO jj=1,4
437 IF(dd(i,jj) <= onep03*dmin)THEN
438 lbx = lb(i,jj)
439 lcx = lc(i,jj)
440 lax = one-lb(i,jj)-lc(i,jj)
441C
442C Privilegier sector in which we are.
443 IF(lbx >= zero .AND. lcx >= zero)THEN
444 lx=zero
445 ELSE
446 ! closest point in the angular sector == center
447 lx = max(zero,dd(i,jj)-bb(i,jj)*bb(i,jj))
448 END IF
449
450 IF(lx < ld(i)) THEN
451 subtria(i)= jj
452 dist(i) = dd(i,jj)
453 ld(i) = lx
454 END IF
455
456 END IF
457 END DO
458 ELSE
459 IF(dd(i,1) <= dist(i))THEN
460 subtria(i)= 1
461 dist(i) = dd(i,1)
462 END IF
463 END IF
464C
465 it = subtria(i)
466 DO j=1,4
467 IF(j /= it) pent(i,j)=zero
468 END DO
469 END DO
470C--------------------------------------------------------
471C PREVENT IMPACT IN A GIVEN DIRECTION IF BOUNDARY CONDITION
472C IN THE SAME DIRECTION FOR BOTH SECONDARY AND PRIMARY NODES
473C--------------------------------------------------------
474 ibcx(1:jlt)=0
475 ibcy(1:jlt)=0
476 ibcz(1:jlt)=0
477 DO i=1,jlt
478C
479 IF(stif(i) <= zero)cycle
480C
481 IF(.NOT.(etyp(i)==0 .OR. etyp(i) > nrtm))cycle ! only solids including coating shell.
482C
483 IF(nsvg(i) > 0)THEN
484 ibcs = icodt(nsvg(i))
485 isks = iskew(nsvg(i))
486 ELSE
487 n = cand_n(i)-nsn
488 ibcs = icodt_fi(nin)%P(n)
489 isks = iskew_fi(nin)%P(n)
490 END IF
491
492 ibcm(1)=icodt(ix1(i))
493 iskm(1)=iskew(ix1(i))
494
495 ibcm(2)=icodt(ix2(i))
496 iskm(2)=iskew(ix2(i))
497
498 ibcm(3)=icodt(ix3(i))
499 iskm(3)=iskew(ix3(i))
500
501 ibcm(4)=icodt(ix4(i))
502 iskm(4)=iskew(ix4(i))
503
504 IF(isks==1)THEN
505 IF((ibcs ==1.OR.ibcs ==3.OR.ibcs ==5.OR.ibcs ==7).AND.
506 . (ibcm(1)==1.OR.ibcm(1)==3.OR.ibcm(1)==5.OR.ibcm(1)==7).AND.
507 . (ibcm(2)==1.OR.ibcm(2)==3.OR.ibcm(2)==5.OR.ibcm(2)==7).AND.
508 . (ibcm(3)==1.OR.ibcm(3)==3.OR.ibcm(3)==5.OR.ibcm(3)==7).AND.
509 . (ibcm(4)==1.OR.ibcm(4)==3.OR.ibcm(4)==5.OR.ibcm(4)==7))THEN
510 ibcz(i)=1
511 END IF
512 IF((ibcs ==2.OR.ibcs ==3.OR.ibcs ==6.OR.ibcs ==7).AND.
513 . (ibcm(1)==2.OR.ibcm(1)==3.OR.ibcm(1)==6.OR.ibcm(1)==7).AND.
514 . (ibcm(2)==2.OR.ibcm(2)==3.OR.ibcm(2)==6.OR.ibcm(2)==7).AND.
515 . (ibcm(3)==2.OR.ibcm(3)==3.OR.ibcm(3)==6.OR.ibcm(3)==7).AND.
516 . (ibcm(4)==2.OR.ibcm(4)==3.OR.ibcm(4)==6.OR.ibcm(4)==7))THEN
517 ibcy(i)=1
518 END IF
519 IF((ibcs ==4.OR.ibcs ==5.OR.ibcs ==6.OR.ibcs ==7).AND.
520 . (ibcm(1)==4.OR.ibcm(1)==5.OR.ibcm(1)==6.OR.ibcm(1)==7).AND.
521 . (ibcm(2)==4.OR.ibcm(2)==5.OR.ibcm(2)==6.OR.ibcm(2)==7).AND.
522 . (ibcm(3)==4.OR.ibcm(3)==5.OR.ibcm(3)==6.OR.ibcm(3)==7).AND.
523 . (ibcm(4)==4.OR.ibcm(4)==5.OR.ibcm(4)==6.OR.ibcm(4)==7))THEN
524 ibcx(i)=1
525 END IF
526 END IF
527 END DO
528C-----
529 prec=one-em04
530 DO i=1,jlt
531C
532 IF(stif(i) <= zero)cycle
533C
534 it = subtria(i)
535 IF(pent(i,it) == zero)cycle
536C
537 IF((ibcz(i)==1.AND.abs(zn(i,it)) > prec).OR.
538 . (ibcy(i)==1.AND.abs(yn(i,it)) > prec).OR.
539 . (ibcx(i)==1.AND.abs(xn(i,it)) > prec))THEN
540 pent(i,it)=zero
541 END IF
542C
543 END DO
544C--------------------------------------------------------
545C Check vs bissectors ...
546C--------------------------------------------------------
547#include "vectorize.inc"
548 DO i =1,jlt
549C
550 IF (stif(i) <= zero)cycle
551C
552 it=subtria(i)
553 IF(it==0)cycle
554 IF(pent(i,it)==zero) cycle
555C
556 i1=itria(1,it)
557 i2=itria(2,it)
558C
559C Projection outside the triangle
560 IF(lb(i,it) < -em03 .OR. lc(i,it) < -em03 .OR.
561 . lb(i,it)+lc(i,it) > one+em03) far(i,it)=1
562C
563 xh=xi(i)+bb(i,it)*xn(i,it)
564 yh=yi(i)+bb(i,it)*yn(i,it)
565 zh=zi(i)+bb(i,it)*zn(i,it)
566C
567 IF(ix3(i) /= ix4(i))THEN
568C
569 ib1=ibound(i1,i)
570 ib2=ibound(i2,i)
571 IF(mvoisn(i,it)==0)THEN
572C
573 IF( (xh-xx(i,i1))* nnx(i,it)+
574 . (yh-yy(i,i1))* nny(i,it)+
575 . (zh-zz(i,i1))* nnz(i,it) >= gaps(i)) far(i,it) =2
576C
577 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
578 . (ib2 /= 0 .AND. ib1 == 0))THEN
579C
580 ibx=max(ib1,ib2)
581 IF(ib1/=0)THEN
582 ix =i1
583 ELSEIF(ib2/=0)THEN
584 ix =i2
585 END IF
586C
587 IF(vtx_bisector(1,1,ibx)/=zero.OR.
588 . vtx_bisector(2,1,ibx)/=zero.OR.
589 . vtx_bisector(3,1,ibx)/=zero.OR.
590 . vtx_bisector(1,2,ibx)/=zero.OR.
591 . vtx_bisector(2,2,ibx)/=zero.OR.
592 . vtx_bisector(3,2,ibx)/=zero)THEN
593 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
594 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
595 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
596 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
597 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
598 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
599 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
600 ELSE
601 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
602 vy = y0n(i,ix)
603 vz = z0n(i,ix)
604 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
605 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
606 IF(pn >= gaps(i)) far(i,it) =2
607 END IF
608C
609 END IF
610C
611C Cone
612 IF(far(i,it)==1 .OR. bb(i,it) <= zero) THEN
613
614 IF(ingap(i,it) == 1 .OR. (kslide(i,i1)/=0.OR.kslide(i,i2)/=0))THEN
615C
616C INGAP = 1 => always check if we are in the cone (allows to determine which cone we are in)
617C
618C sliding => use the cone to choose between neighboring segments
619C but keep the contact if the node has slipped more than one element (see supersonic)
620C
621 x12(i)=xx(i,i2)-xx(i,i1)
622 y12(i)=yy(i,i2)-yy(i,i1)
623 z12(i)=zz(i,i2)-zz(i,i1)
624C
625C normal to the bisecting plane (pointing toward the inside)
626 px = z12(i)*nny(i,it)-y12(i)*nnz(i,it)
627 py = x12(i)*nnz(i,it)-z12(i)*nnx(i,it)
628 pz = y12(i)*nnx(i,it)-x12(i)*nny(i,it)
629 pp = px*px+py*py+pz*pz
630C
631 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
632C
633 side(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/max(em30,ll*pp))
634 IF(side(i,it) < -zep01) far(i,it) =2
635C
636 END IF
637
638 END IF
639
640 ELSE
641
642 ib1=ibound(1,i)
643 ib2=ibound(2,i)
644 ib3=ibound(3,i)
645C
646 d1=(xh-xx(i,1))* nnx(i,1)+
647 . (yh-yy(i,1))* nny(i,1)+
648 . (zh-zz(i,1))* nnz(i,1)
649 d2=(xh-xx(i,2))* nnx(i,2)+
650 . (yh-yy(i,2))* nny(i,2)+
651 . (zh-zz(i,2))* nnz(i,2)
652 d3=(xh-xx(i,3))* nnx(i,4)+
653 . (yh-yy(i,3))* nny(i,4)+
654 . (zh-zz(i,3))* nnz(i,4)
655C
656 IF( (mvoisn(i,1) == 0 .AND. d1 >= gaps(i)).OR.
657 . (mvoisn(i,2) == 0 .AND. d2 >= gaps(i)).OR.
658 . (mvoisn(i,4) == 0 .AND. d3 >= gaps(i)) )THEN
659C
660 far(i,it)=2
661C
662 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
663 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
664 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
665C
666 ibx=max(ib1,ib2,ib3)
667 IF(ib1/=0)THEN
668 ix =1
669 iy =2
670 iz =3
671 ELSEIF(ib2/=0)THEN
672 ix =2
673 iy =3
674 iz =1
675 ELSEIF(ib3/=0)THEN
676 ix =3
677 iy =1
678 iz =2
679 END IF
680C
681 IF(vtx_bisector(1,1,ibx)/=zero.OR.
682 . vtx_bisector(2,1,ibx)/=zero.OR.
683 . vtx_bisector(3,1,ibx)/=zero.OR.
684 . vtx_bisector(1,2,ibx)/=zero.OR.
685 . vtx_bisector(2,2,ibx)/=zero.OR.
686 . vtx_bisector(3,2,ibx)/=zero)THEN
687 p1 = (xh-xx(i,ix))* vtx_bisector(1,1,ibx)+
688 . (yh-yy(i,ix))* vtx_bisector(2,1,ibx)+
689 . (zh-zz(i,ix))* vtx_bisector(3,1,ibx)
690 p2 = (xh-xx(i,ix))* vtx_bisector(1,2,ibx)+
691 . (yh-yy(i,ix))* vtx_bisector(2,2,ibx)+
692 . (zh-zz(i,ix))* vtx_bisector(3,2,ibx)
693 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
694 ELSE
695 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
696 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
697 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
698 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
699 pn = ((xh-xx(i,ix))*vx+(yh-yy(i,ix))*vy+(zh-zz(i,ix))*vz)*nn
700 IF(pn >= gaps(i)) far(i,it) =2
701 END IF
702C
703 END IF
704C
705 IF(far(i,it)==1 .OR. bb(i,it) <= zero) THEN
706
707 IF(mvoisn(i,1) /= 0 .AND. (ingap(i,it) == 1 .OR. (kslide(i,1)/=0.OR.kslide(i,2)/=0)))THEN
708C
709C INGAP = 1 => always check if we are in the cone (allows to determine which cone we are in)
710C
711C sliding => use the cone to choose between neighboring segments
712C but keep the contact if the node has slipped more than one element (see supersonic)
713C
714 x12(i)=xx(i,2)-xx(i,1)
715 y12(i)=yy(i,2)-yy(i,1)
716 z12(i)=zz(i,2)-zz(i,1)
717C
718C normal to the bisecting plane (pointing toward the inside)
719 px = z12(i)*nny(i,1)-y12(i)*nnz(i,1)
720 py = x12(i)*nnz(i,1)-z12(i)*nnx(i,1)
721 pz = y12(i)*nnx(i,1)-x12(i)*nny(i,1)
722 pp = px*px+py*py+pz*pz
723C
724 ll = xij(i,1)*xij(i,1)+yij(i,1)*yij(i,1)+zij(i,1)*zij(i,1)
725C
726 side(i,1)=-(xij(i,1)*px+yij(i,1)*py+zij(i,1)*pz)*sqrt(one/max(em30,ll*pp))
727 IF(side(i,1) < -zep01) far(i,it) =2
728C
729 END IF
730
731 IF(mvoisn(i,2) /= 0 .AND. (ingap(i,it) == 1 .OR. (kslide(i,2)/=0.OR.kslide(i,3)/=0)))THEN
732C
733C INGAP = 1 => always check if we are in the cone (allows to determine which cone we are in)
734C
735C sliding => use the cone to choose between neighboring segments
736C but keep the contact if the node has slipped more than one element (see supersonic)
737C
738 x12(i)=xx(i,3)-xx(i,2)
739 y12(i)=yy(i,3)-yy(i,2)
740 z12(i)=zz(i,3)-zz(i,2)
741C
742C normal to the bisecting plane (pointing toward the inside)
743 px = z12(i)*nny(i,2)-y12(i)*nnz(i,2)
744 py = x12(i)*nnz(i,2)-z12(i)*nnx(i,2)
745 pz = y12(i)*nnx(i,2)-x12(i)*nny(i,2)
746 pp = px*px+py*py+pz*pz
747C
748 ll = xij(i,2)*xij(i,2)+yij(i,2)*yij(i,2)+zij(i,2)*zij(i,2)
749C
750 side(i,2)=-(xij(i,2)*px+yij(i,2)*py+zij(i,2)*pz)*sqrt(one/max(em30,ll*pp))
751 IF(side(i,2) < -zep01) far(i,it) =2
752C
753 END IF
754
755 IF(mvoisn(i,4) /= 0 .AND. (ingap(i,it) == 1 .OR. (kslide(i,3)/=0.OR.kslide(i,1)/=0)))THEN
756C
757C INGAP = 1 => always check if we are in the cone (allows to determine which cone we are in)
758C
759C sliding => use the cone to choose between neighboring segments
760C but keep the contact if the node has slipped more than one element (see supersonic)
761C
762 x12(i)=xx(i,1)-xx(i,3)
763 y12(i)=yy(i,1)-yy(i,3)
764 z12(i)=zz(i,1)-zz(i,3)
765C
766C normal to the bisecting plane (pointing toward the inside)
767 px = z12(i)*nny(i,4)-y12(i)*nnz(i,4)
768 py = x12(i)*nnz(i,4)-z12(i)*nnx(i,4)
769 pz = y12(i)*nnx(i,4)-x12(i)*nny(i,4)
770 pp = px*px+py*py+pz*pz
771C
772 ll = xij(i,3)*xij(i,3)+yij(i,3)*yij(i,3)+zij(i,3)*zij(i,3)
773C
774 side(i,4)=-(xij(i,3)*px+yij(i,3)*py+zij(i,3)*pz)*sqrt(one/max(em30,ll*pp))
775 IF(side(i,4) < -zep01) far(i,it) =2
776C
777 END IF
778 END IF
779 END IF
780C
781 IF(far(i,it)==2) pent(i,it)=zero
782C
783 ENDDO
784C--------------------------------------------------------
785 DO i=1,jlt
786C
787 IF(stif(i) <= zero)cycle
788C
789 it = subtria(i)
790 IF(it/=0.AND.pent(i,it)==zero) it=0
791
792 IF(it == 0)cycle
793
794 n = cand_n(i)
795 l = cand_e(i)
796
797 mglob=mseglo(l)
798
799 pene=pent(i,it) ! PENE /= ZERO
800 IF(n<=nsn)THEN
801#include "lockon.inc"
802
803c if(itab(nsvg(i))==27363)
804cc if(itab(nsvg(i))==27952.or.
805cc . itab(ix1(i))==27952.or.itab(ix2(i))==27952.or.itab(ix3(i))==27952.or.itab(ix4(i))==27952)
806c . print *,'dst21 nat',ispmd+1,it,irtlm(1,n),mglob,cand_e(i),pent(i,it),far(i,it),
807c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
808c . ibound(1:4,i),lb(i,it),lc(i,it),lb(i,it)+lc(i,it),
809c . BB(i,it),far(i,it),TIME_S(1,N),time_s(2,n),
810c . dist(i),ld(i)
811 IF(dist(i) < onep02*time_s(1,n) .AND.
812 * (dist(i) < time_s(2,n) .OR.
813 * (dist(i) == time_s(2,n) .AND. irtlm(1,n) < mglob)))THEN
814 irtlm(1,n) = mglob
815 irtlm(2,n) = it
816 irtlm(3,n) = cand_e(i)
817 irtlm(4,n) = ispmd+1
818 time_s(2,n) = dist(i)
819 END IF
820#include "lockoff.inc"
821 ELSE
822#include "lockon.inc"
823 IF(dist(i) < onep02*time_sfi(nin)%P(2*(n-nsn-1)+1) .AND.
824 * (dist(i) < time_sfi(nin)%P(2*(n-nsn-1)+2) .OR.
825 * (dist(i) == time_sfi(nin)%P(2*(n-nsn-1)+2) .AND. irtlm_fi(nin)%P(1,n-nsn) < mglob)))THEN
826 irtlm_fi(nin)%P(1,n-nsn) = mglob
827 irtlm_fi(nin)%P(2,n-nsn) = it
828 irtlm_fi(nin)%P(3,n-nsn) = cand_e(i)
829 irtlm_fi(nin)%P(4,n-nsn) = ispmd+1
830 time_sfi(nin)%P(2*(n-nsn-1)+2) = dist(i)
831 END IF
832c if(itafi(nin)%p(n-nsn)==3817238)
833cc if(itafi(nin)%p(n-nsn)==3817238.or.
834cc . itab(ix1(i))==3817238.or.itab(ix2(i))==3817238.or.itab(ix3(i))==3817238.or.itab(ix4(i))==3817238)
835c . print *,'dst21 rem',ispmd+1,it,IRTLM_FI(NIN)%P(1,n-nsn),mglob,cand_e(i),pent(i,it),far(i,it),
836c . TIME_SFI(NIN)%P(N-NSN),dist(i),BB(i,it),
837c . itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i))
838#include "lockoff.inc"
839 END IF
840 END DO
841
842 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
tri7box::iskew_fi
type(int_pointer), dimension(:), allocatable iskew_fi
Definition tri7box.F:550
tri7box::time_sfi
type(real_pointer), dimension(:), allocatable time_sfi
Definition tri7box.F:542
tri7box::irtlm_fi
type(int_pointer2), dimension(:), allocatable irtlm_fi
Definition tri7box.F:533
tri7box::icodt_fi
type(int_pointer), dimension(:), allocatable icodt_fi
Definition tri7box.F:551