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 . INTTH,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, ITQ, N4N, IKEEP, 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 . IBC21, IBC22, IBC23, IBCS, ISKS, IBCM(4), ISKM(4),
83 . IBCX(MVSIZ), IBCY(MVSIZ), IBCZ(MVSIZ)
84 my_real
85 . xij(mvsiz,4), xi0v(mvsiz), xi5,
86 . yij(mvsiz,4), yi0v(mvsiz), yi5,
87 . zij(mvsiz,4), zi0v(mvsiz), zi5,
88 . x51, x52, x53, x54,
89 . y51, y52, y53, y54,
90 . z51, z52, z53, z54,
91 . xo1, xo2, xo3, xo4, xo5, xoi,
92 . yo1, yo2, yo3, yo4, yo5, yoi,
93 . zo1, zo2, zo3, zo4, zo5, zoi,
94 . vo12, vo23, vo34, vo41, 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 . unp,zerom,eps,unpt,zeromt,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 INTEGER ISHEL(MVSIZ)
114 DATA itria/1,2,2,3,3,4,4,1/
115C-----------------------------------------------------------------------
116C
117C initialisation (cf triangles)
118 far(1:jlt,1:4) = 0
119 pent(1:jlt,1:4) = ep20
120 dd(1:jlt,1:4) = ep20
121 dist(1:jlt)= ep20
122 ld(1:jlt) = ep20
123C
124 DO i=1,jlt
125C
126C For computing LBP, LCP, FAR=2, etc
127C
128 IF(stif(i) <= zero)cycle
129C
130 x0n(i,1) = xx(i,1) - xx(i,5)
131 y0n(i,1) = yy(i,1) - yy(i,5)
132 z0n(i,1) = zz(i,1) - zz(i,5)
133C
134 x0n(i,2) = xx(i,2) - xx(i,5)
135 y0n(i,2) = yy(i,2) - yy(i,5)
136 z0n(i,2) = zz(i,2) - zz(i,5)
137C
138 x0n(i,3) = xx(i,3) - xx(i,5)
139 y0n(i,3) = yy(i,3) - yy(i,5)
140 z0n(i,3) = zz(i,3) - zz(i,5)
141C
142 x0n(i,4) = xx(i,4) - xx(i,5)
143 y0n(i,4) = yy(i,4) - yy(i,5)
144 z0n(i,4) = zz(i,4) - zz(i,5)
145C
146 IF(ix3(i)/=ix4(i))THEN
147 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
148 ELSE
149 gap_mm(i)=gap_nm(3,i)
150 END IF
151 ENDDO
152C--------------------------------------------------------
153#include "vectorize.inc"
154 DO i=1,jlt
155C
156 IF(stif(i) <= zero)cycle
157C
158 xi0v(i) = xx(i,5) - xi(i)
159 yi0v(i) = yy(i,5) - yi(i)
160 zi0v(i) = zz(i,5) - zi(i)
161C
162 xij(i,1) = xx(i,1) - xi(i)
163 yij(i,1) = yy(i,1) - yi(i)
164 zij(i,1) = zz(i,1) - zi(i)
165C
166 xij(i,2) = xx(i,2) - xi(i)
167 yij(i,2) = yy(i,2) - yi(i)
168 zij(i,2) = zz(i,2) - zi(i)
169C
170 xij(i,3) = xx(i,3) - xi(i)
171 yij(i,3) = yy(i,3) - yi(i)
172 zij(i,3) = zz(i,3) - zi(i)
173C
174 xij(i,4) = xx(i,4) - xi(i)
175 yij(i,4) = yy(i,4) - yi(i)
176 zij(i,4) = zz(i,4) - zi(i)
177C
178 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
179 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
180 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
181C
182 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
183 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
184 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
185C
186 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
187 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
188 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
189 nn = one/
190 . max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
191 xn(i,1)=xn(i,1)*nn
192 yn(i,1)=yn(i,1)*nn
193 zn(i,1)=zn(i,1)*nn
194 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
195 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
196C
197 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
198 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
199 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
200C
201 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
202 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
203 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
204 nn = one/
205 . max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i,2)+ zn(i,2)*zn(i,2)))
206 xn(i,2)=xn(i,2)*nn
207 yn(i,2)=yn(i,2)*nn
208 zn(i,2)=zn(i,2)*nn
209 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn(i,2)*sz3)*nn
210 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
211C
212 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
213 sy4 = zi0v(i)*xij(i,4) - xi0v(i)*zij(i,4)
214 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
215C
216 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
217 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
218 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
219 nn = one/
220 . max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
221 xn(i,3)=xn(i,3)*nn
222 yn(i,3)=yn(i,3)*nn
223 zn(i,3)=zn(i,3)*nn
224 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
225 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
226C
227 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
228 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
229 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
230 nn = one/
231 . max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn(i,4)))
232 xn(i,4)=xn(i,4)*nn
233 yn(i,4)=yn(i,4)*nn
234 zn(i,4)=zn(i,4)*nn
235 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
236 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
237C
238 END DO
239C--------------------------------------------------------
240C
241#include "vectorize.inc"
242 DO i=1,jlt
243C
244 IF(stif(i) <= zero)cycle
245C
246 aaa = one/max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
247 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
248 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
249 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
250 al(i,1) = max(zero,min(one,al(i,1)))
251 aaa = one/max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
252 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
253 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
254 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
255 al(i,2) = max(zero,min(one,al(i,2)))
256 aaa = one/max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
257 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
258 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
259 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
260 al(i,3) = max(zero,min(one,al(i,3)))
261 aaa = one/max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
262 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
263 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
264 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
265 al(i,4) = max(zero,min(one,al(i,4)))
266C
267 END DO
268C--------------------------------------------------------
269 ingap(1:jlt,1:4) = 0
270C--------------------------------------------------------
271 it=1
272 i1=itria(1,it)
273 i2=itria(2,it)
274#include "vectorize.inc"
275 DO i=1,jlt
276C
277 IF(stif(i) <= zero)cycle
278C
279 x12(i) = xx(i,i2) - xx(i,i1)
280 y12(i) = yy(i,i2) - yy(i,i1)
281 z12(i) = zz(i,i2) - zz(i,i1)
282C
283 lbp(i,it) = lb(i,it)
284 lcp(i,it) = lc(i,it)
285 lap = one-lbp(i,it)-lcp(i,it)
286C HLA, HLB, HLC necessaires pour triangle angle obtu
287 aaa = one / max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
288 hla= lap*abs(lap) * aaa
289 IF(lap<zero.AND.
290 + hla<=hlb(i,it).AND.hla<=hlc(i,it))THEN
291 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
292 lbp(i,it) = max(zero,min(one,lbp(i,it)))
293 lcp(i,it) = one - lbp(i,it)
294 ELSEIF(lbp(i,it)<zero.AND.
295 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)THEN
296 lbp(i,it) = zero
297 lcp(i,it) = al(i,i2)
298 ELSEIF(lcp(i,it)<zero.AND.
299 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))THEN
300 lcp(i,it) = zero
301 lbp(i,it) = al(i,i1)
302 ENDIF
303
304 lap = one-lbp(i,it)-lcp(i,it)
305 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
306 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
307 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
308 dx =xi(i)-xp
309 dy =yi(i)-yp
310 dz =zi(i)-zp
311 dd(i,it)=dx*dx+dy*dy+dz*dz
312
313 END DO
314C--------------------------------------------------------
315C Verifie si penetre vs gap cylindrique ...
316C--------------------------------------------------------
317#include "vectorize.inc"
318 DO i =1,jlt
319C
320 IF (stif(i) <= zero)cycle
321C
322 lap = one-lbp(i,it)-lcp(i,it)
323C
324 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),
325 . max(drad,gapmxl(i)+dgapload))
326 gap2 = gap**2
327C
328 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))
329C
330 IF(dd(i,it) <= gap2 .AND. bb(i,it) <= zero) ingap(i,it)=1
331
332 IF(bb(i,it) > zero)THEN
333C
334C penetration approximee (cf zone limite interpolation des normales)
335 pent(i,it)=max(zero,gap+bb(i,it))
336 ELSE
337 pent(i,it)=max(zero,gap-sqrt(dd(i,it)))
338 END IF
339C
340 ENDDO
341C--------------------------------------------------------
342 n4n=0
343 DO i=1,jlt
344C
345 IF(stif(i) <= zero)cycle
346C
347 IF(ix3(i)/=ix4(i))THEN
348 n4n = n4n+1
349 i4n(n4n)=i
350 ENDIF
351 ENDDO
352C--------------------------------------------------------
353 DO it=2,4
354
355 i1=itria(1,it)
356 i2=itria(2,it)
357
358#include "vectorize.inc"
359 DO k=1,n4n
360 i=i4n(k)
361C
362 x12(i) = xx(i,i2) - xx(i,i1)
363 y12(i) = yy(i,i2) - yy(i,i1)
364 z12(i) = zz(i,i2) - zz(i,i1)
365C
366 lbp(i,it) = lb(i,it)
367 lcp(i,it) = lc(i,it)
368 lap = one-lbp(i,it)-lcp(i,it)
369C HLA, HLB, HLC necessaires pour triangle angle obtu
370 aaa = one / max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
371 hla= lap*abs(lap) * aaa
372 IF(lap<zero.AND.
373 + hla<=hlb(i,it).AND.hla<=hlc(i,it))THEN
374 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
375 lbp(i,it) = max(zero,min(one,lbp(i,it)))
376 lcp(i,it) = one - lbp(i,it)
377 ELSEIF(lbp(i,it)<zero.AND.
378 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)THEN
379 lbp(i,it) = zero
380 lcp(i,it) = al(i,i2)
381 ELSEIF(lcp(i,it)<zero.AND.
382 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))THEN
383 lcp(i,it) = zero
384 lbp(i,it) = al(i,i1)
385 ENDIF
386
387 lap = one-lbp(i,it)-lcp(i,it)
388 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
389 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
390 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
391 dx =xi(i)-xp
392 dy =yi(i)-yp
393 dz =zi(i)-zp
394 dd(i,it)=dx*dx+dy*dy+dz*dz
395
396 END DO
397C--------------------------------------------------------
398C Verifie si penetre vs gap cylindrique ...
399C--------------------------------------------------------
400#include "vectorize.inc"
401 DO k=1,n4n
402 i=i4n(k)
403C
404 lap = one-lbp(i,it)-lcp(i,it)
405C
406 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),
407 . max(drad,gapmxl(i)+dgapload))
408 gap2 = gap**2
409C
410 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))
411C
412 IF(dd(i,it) <= gap2 .AND. bb(i,it) <= zero) ingap(i,it)=1
413
414 IF(bb(i,it) > zero)THEN
415C
416C penetration approximee (cf zone limite interpolation des normales)
417 pent(i,it)=max(zero,gap+bb(i,it))
418 ELSE
419 pent(i,it)=max(zero,gap-sqrt(dd(i,it)))
420 END IF
421C
422 ENDDO
423 END DO ! DO IT=2,4
424C--------------------------------------------------------
425C Look for closest sub-triangle
426C--------------------------------------------------------
427 DO i=1,jlt
428C
429 IF(stif(i) <= zero)cycle
430C
431C Erase Subtria (cf intersections)
432 subtria(i)=0
433C
434 IF(ix3(i)/=ix4(i))THEN
435 dmin=min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
436
437 DO jj=1,4
438 IF(dd(i,jj) <= onep03*dmin)THEN
439 lbx = lb(i,jj)
440 lcx = lc(i,jj)
441 lax = one-lb(i,jj)-lc(i,jj)
442C
443C Privilegier secteur dans lequel on se trouve.
444 IF(lbx >= zero .AND. lcx >= zero)THEN
445 lx=zero
446 ELSE
447 ! point le plus proche dans le secteur angulaire == centre
448 lx = max(zero,dd(i,jj)-bb(i,jj)*bb(i,jj))
449 END IF
450
451 IF(lx < ld(i)) THEN
452 subtria(i)= jj
453 dist(i) = dd(i,jj)
454 ld(i) = lx
455 END IF
456
457 END IF
458 END DO
459 ELSE
460 IF(dd(i,1) <= dist(i))THEN
461 subtria(i)= 1
462 dist(i) = dd(i,1)
463 END IF
464 END IF
465C
466 it = subtria(i)
467 DO j=1,4
468 IF(j /= it) pent(i,j)=zero
469 END DO
470 END DO
471C--------------------------------------------------------
472C PREVENT IMPACT IN A GIVEN DIRECTION IF BOUNDARY CONDITION
473C IN THE SAME DIRECTION FOR BOTH SECONDARY AND PRIMARY NODES
474C--------------------------------------------------------
475 ibcx(1:jlt)=0
476 ibcy(1:jlt)=0
477 ibcz(1:jlt)=0
478 DO i=1,jlt
479C
480 IF(stif(i) <= zero)cycle
481C
482 IF(.NOT.(etyp(i)==0 .OR. etyp(i) > nrtm))cycle ! only solids including coating shell.
483C
484 IF(nsvg(i) > 0)THEN
485 ibcs = icodt(nsvg(i))
486 isks = iskew(nsvg(i))
487 ELSE
488 n = cand_n(i)-nsn
489 ibcs = icodt_fi(nin)%P(n)
490 isks = iskew_fi(nin)%P(n)
491 END IF
492
493 ibcm(1)=icodt(ix1(i))
494 iskm(1)=iskew(ix1(i))
495
496 ibcm(2)=icodt(ix2(i))
497 iskm(2)=iskew(ix2(i))
498
499 ibcm(3)=icodt(ix3(i))
500 iskm(3)=iskew(ix3(i))
501
502 ibcm(4)=icodt(ix4(i))
503 iskm(4)=iskew(ix4(i))
504
505 IF(isks==1)THEN
506 IF((ibcs ==1.OR.ibcs ==3.OR.ibcs ==5.OR.ibcs ==7).AND.
507 . (ibcm(1)==1.OR.ibcm(1)==3.OR.ibcm(1)==5.OR.ibcm(1)==7).AND.
508 . (ibcm(2)==1.OR.ibcm(2)==3.OR.ibcm(2)==5.OR.ibcm(2)==7).AND.
509 . (ibcm(3)==1.OR.ibcm(3)==3.OR.ibcm(3)==5.OR.ibcm(3)==7).AND.
510 . (ibcm(4)==1.OR.ibcm(4)==3.OR.ibcm(4)==5.OR.ibcm(4)==7))THEN
511 ibcz(i)=1
512 END IF
513 IF((ibcs ==2.OR.ibcs ==3.OR.ibcs ==6.OR.ibcs ==7).AND.
514 . (ibcm(1)==2.OR.ibcm(1)==3.OR.ibcm(1)==6.OR.ibcm(1)==7).AND.
515 . (ibcm(2)==2.OR.ibcm(2)==3.OR.ibcm(2)==6.OR.ibcm(2)==7).AND.
516 . (ibcm(3)==2.OR.ibcm(3)==3.OR.ibcm(3)==6.OR.ibcm(3)==7).AND.
517 . (ibcm(4)==2.OR.ibcm(4)==3.OR.ibcm(4)==6.OR.ibcm(4)==7))THEN
518 ibcy(i)=1
519 END IF
520 IF((ibcs ==4.OR.ibcs ==5.OR.ibcs ==6.OR.ibcs ==7).AND.
521 . (ibcm(1)==4.OR.ibcm(1)==5.OR.ibcm(1)==6.OR.ibcm(1)==7).AND.
522 . (ibcm(2)==4.OR.ibcm(2)==5.OR.ibcm(2)==6.OR.ibcm(2)==7).AND.
523 . (ibcm(3)==4.OR.ibcm(3)==5.OR.ibcm(3)==6.OR.ibcm(3)==7).AND.
524 . (ibcm(4)==4.OR.ibcm(4)==5.OR.ibcm(4)==6.OR.ibcm(4)==7))THEN
525 ibcx(i)=1
526 END IF
527 END IF
528 END DO
529C-----
530 prec=one-em04
531 DO i=1,jlt
532C
533 IF(stif(i) <= zero)cycle
534C
535 it = subtria(i)
536 IF(pent(i,it) == zero)cycle
537C
538 IF((ibcz(i)==1.AND.abs(zn(i,it)) > prec).OR.
539 . (ibcy(i)==1.AND.abs(yn(i,it)) > prec).OR.
540 . (ibcx(i)==1.AND.abs(xn(i,it)) > prec))THEN
541 pent(i,it)=zero
542 END IF
543C
544 END DO
545C--------------------------------------------------------
546C Check vs bissectors ...
547C--------------------------------------------------------
548#include "vectorize.inc"
549 DO i =1,jlt
550C
551 IF (stif(i) <= zero)cycle
552C
553 it=subtria(i)
554 IF(pent(i,it)==zero) cycle
555C
556 i1=itria(1,it)
557 i2=itria(2,it)
558C
559C Projection en dehors du 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 => tjrs regarder si l'on est dans le cone (permet de determiner dans quel cone on se trouve)
617C
618C Glissement => utiliser le cone pour choisir entre les segments voisins
619C mais garder le contact si le nd a glisse de plus d'un elt (cf 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 => tjrs regarder si l'on est dans le cone (permet de determiner dans quel cone on se trouve)
710C
711C Glissement => utiliser le cone pour choisir entre les segments voisins
712C mais garder le contact si le nd a glisse de plus d'un elt (cf 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 => tjrs regarder si l'on est dans le cone (permet de determiner dans quel cone on se trouve)
734C
735C Glissement => utiliser le cone pour choisir entre les segments voisins
736C mais garder le contact si le nd a glisse de plus d'un elt (cf 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 => tjrs regarder si l'on est dans le cone (permet de determiner dans quel cone on se trouve)
758C
759C Glissement => utiliser le cone pour choisir entre les segments voisins
760C mais garder le contact si le nd a glisse de plus d'un elt (cf 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