OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25dst3_21.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| i25dst3_21 ../engine/source/interfaces/int25/i25dst3_21.F
25!||--- called by ------------------------------------------------------
26!|| i25comp_2 ../engine/source/interfaces/int25/i25comp_2.F
27!||--- uses -----------------------------------------------------
28!|| tri7box ../engine/share/modules/tri7box.F
29!||====================================================================
30 SUBROUTINE i25dst3_21(
31 1 JLT ,CAND_N ,CAND_E ,NRTM ,XX ,
32 2 YY ,ZZ ,XI ,YI ,ZI ,
33 3 NIN ,NSN ,IX1 ,IX2 ,IX3 ,
34 4 IX4 ,NSVG ,STIF ,INACTI ,MSEGLO ,
35 5 GAPS ,GAPM ,IRECT ,IRTLM ,TIME_S ,
36 6 GAP_NM ,ITAB ,NNX ,NNY ,NNZ ,
37 7 FAR ,PENT ,DIST ,LB ,LC ,
38 8 LBP ,LCP ,KSLIDE ,MVOISN ,GAPMXL ,
39 9 IBOUND ,VTX_BISECTOR ,ETYP ,ICODT ,ISKEW,
40 A DRAD ,DGAPLOAD)
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
843 END
i25dst3_21
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)
Definition i25dst3_21.F:41
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