OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25dst3_1.F File Reference
#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "comlock.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i25dst3_1 (jlt, cand_n, cand_e, xx, yy, zz, xi, yi, zi, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, mseglo, gaps, gapm, gapmxl, irect, irtlm, time_s, gap_nm, itab, icont_i, nnx, nny, nnz, far, pent, dist, lb, lc, lbp, lcp, subtria, mvoisn, ibound, vtx_bisector, drad, dgapload)
subroutine i25glob_1 (jlt, cand_n, cand_e, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, mseglo, irtlm, time_s, itab, far, pent, lb, lc, farm, penm, lbm, lcm)
subroutine i25glob (jlt, cand_n, cand_e, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, mseglo, irtlm, time_s, itab, far, pent, lb, lc, index, farm, penm, lbm, lcm)

Function/Subroutine Documentation

◆ i25dst3_1()

subroutine i25dst3_1 ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
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,
gapmxl,
integer, dimension(4,*) irect,
integer, dimension(4,nsn) irtlm,
time_s,
gap_nm,
integer, dimension(*) itab,
integer, dimension(*) icont_i,
nnx,
nny,
nnz,
integer, dimension(mvsiz,4) far,
pent,
dist,
lb,
lc,
lbp,
lcp,
integer, dimension(mvsiz) subtria,
integer, dimension(mvsiz,4) mvoisn,
integer, dimension(4,*) ibound,
real*4, dimension(3,2,*) vtx_bisector,
intent(in) drad,
intent(in) dgapload )

Definition at line 30 of file i25dst3_1.F.

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

◆ i25glob()

subroutine i25glob ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
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,
integer, dimension(4,nsn) irtlm,
time_s,
integer, dimension(*) itab,
integer, dimension(mvsiz,4) far,
pent,
lb,
lc,
integer, dimension(*) index,
integer, dimension(4,*) farm,
penm,
lbm,
lcm )

Definition at line 890 of file i25dst3_1.F.

898C-----------------------------------------------
899C M o d u l e s
900C-----------------------------------------------
901 USE tri7box
902C-----------------------------------------------
903C I m p l i c i t T y p e s
904C-----------------------------------------------
905#include "implicit_f.inc"
906C-----------------------------------------------
907C G l o b a l P a r a m e t e r s
908C-----------------------------------------------
909#include "mvsiz_p.inc"
910#include "comlock.inc"
911C-----------------------------------------------
912C D u m m y A r g u m e n t s
913C-----------------------------------------------
914 INTEGER JLT, NIN, NSN, INACTI,
915 . CAND_N(*),
916 . CAND_E(*),ITAB(*)
917 INTEGER NSVG(MVSIZ), IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
918 . FAR(MVSIZ,4), MSEGLO(*), IRTLM(4,NSN) ,INDEX(*), FARM(4,*)
919 my_real
920 . stif(*), time_s(*),
921 . pent(mvsiz,4), lb(mvsiz,4), lc(mvsiz,4),
922 . penm(4,*), lbm(4,*), lcm(4,*)
923C-----------------------------------------------
924C L o c a l V a r i a b l e s
925C-----------------------------------------------
926 INTEGER IT, I, J
927C--------------------------------------------------------
928C Back to global arrays
929C--------------------------------------------------------
930 DO i=1,jlt
931 j=index(i)
932 DO it=1,4
933 farm(it,j)=far(i,it)
934 penm(it,j)=pent(i,it)
935 lbm(it,j) =lb(i,it)
936 lcm(it,j) =lc(i,it)
937 END DO
938 END DO
939 RETURN

◆ i25glob_1()

subroutine i25glob_1 ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
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,
integer, dimension(4,nsn) irtlm,
time_s,
integer, dimension(*) itab,
integer, dimension(mvsiz,4) far,
pent,
lb,
lc,
integer, dimension(4,*) farm,
penm,
lbm,
lcm )

Definition at line 838 of file i25dst3_1.F.

845C-----------------------------------------------
846C I m p l i c i t T y p e s
847C-----------------------------------------------
848#include "implicit_f.inc"
849C-----------------------------------------------
850C G l o b a l P a r a m e t e r s
851C-----------------------------------------------
852#include "mvsiz_p.inc"
853#include "comlock.inc"
854C-----------------------------------------------
855C D u m m y A r g u m e n t s
856C-----------------------------------------------
857 INTEGER JLT, NIN, NSN, INACTI,
858 . CAND_N(*),
859 . CAND_E(*),ITAB(*)
860 INTEGER NSVG(MVSIZ), IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
861 . FAR(MVSIZ,4), MSEGLO(*), IRTLM(4,NSN), FARM(4,*)
862 my_real
863 . stif(*), time_s(*),
864 . pent(mvsiz,4), lb(mvsiz,4), lc(mvsiz,4),
865 . penm(4,*), lbm(4,*), lcm(4,*)
866C-----------------------------------------------
867C L o c a l V a r i a b l e s
868C-----------------------------------------------
869 INTEGER IT, I
870C--------------------------------------------------------
871C Back to global arrays
872C--------------------------------------------------------
873 DO i=1,jlt
874 DO it=1,4
875 farm(it,i)=far(i,it)
876 penm(it,i)=pent(i,it)
877 lbm(it,i) =lb(i,it)
878 lcm(it,i) =lc(i,it)
879 END DO
880 END DO
881 RETURN