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