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

907C-----------------------------------------------
908C M o d u l e s
909C-----------------------------------------------
910 USE tri7box
911C-----------------------------------------------
912C I m p l i c i t T y p e s
913C-----------------------------------------------
914#include "implicit_f.inc"
915C-----------------------------------------------
916C G l o b a l P a r a m e t e r s
917C-----------------------------------------------
918#include "mvsiz_p.inc"
919#include "comlock.inc"
920C-----------------------------------------------
921C D u m m y A r g u m e n t s
922C-----------------------------------------------
923 INTEGER JLT, NIN, NSN, INACTI,
924 . CAND_N(*),
925 . CAND_E(*),ITAB(*)
926 INTEGER NSVG(MVSIZ), IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
927 . FAR(MVSIZ,4), MSEGLO(*), IRTLM(4,NSN) ,INDEX(*), FARM(4,*)
928 my_real
929 . stif(*), time_s(*),
930 . pent(mvsiz,4), lb(mvsiz,4), lc(mvsiz,4),
931 . penm(4,*), lbm(4,*), lcm(4,*)
932C-----------------------------------------------
933C L o c a l V a r i a b l e s
934C-----------------------------------------------
935 INTEGER IT, I, J, L, N, MGLOB
936C--------------------------------------------------------
937C Back to global arrays
938C--------------------------------------------------------
939 DO i=1,jlt
940 j=index(i)
941 DO it=1,4
942 farm(it,j)=far(i,it)
943 penm(it,j)=pent(i,it)
944 lbm(it,j) =lb(i,it)
945 lcm(it,j) =lc(i,it)
946 END DO
947 END DO
948 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 847 of file i25dst3_1.F.

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