OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25dst3_3.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_3 ../engine/source/interfaces/int25/i25dst3_3.F
25!||--- called by ------------------------------------------------------
26!|| i25mainf ../engine/source/interfaces/int25/i25mainf.F
27!||--- uses -----------------------------------------------------
28!|| tri7box ../engine/share/modules/tri7box.F
29!||====================================================================
30 SUBROUTINE i25dst3_3(
31 1 JLT ,CAND_N ,CAND_E ,CN_LOC ,CE_LOC ,
32 2 IRTLM ,XX ,YY ,ZZ ,GAP_NM ,
33 3 XI ,YI ,ZI ,GAPS ,GAPMXL ,
34 4 ISHARP ,NNX ,NNY ,NNZ ,
35 5 N1 ,N2 ,N3 ,H1 ,H2 ,
36 5 H3 ,H4 ,NIN ,NSN ,IX1 ,
37 6 IX2 ,IX3 ,IX4 ,NSVG ,STIF ,
38 7 INACTI ,KINI ,ITAB ,LB ,LC ,
39 8 PENMIN ,EPS ,PENE ,PENE_OLD ,SUBTRIA ,
40 9 GAPV ,IVIS2 ,IF_ADH ,IFADHI ,BASE_ADH ,
41 A MVOISN ,IBOUND ,VTX_BISECTOR,DIST ,TIME )
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, IVIS2, ISHARP, ITAB(*),
59 . CAND_N(*),CAND_E(*),
60 . MVOISN(MVSIZ,4),IBOUND(4,*)
61 INTEGER NSVG(MVSIZ), KINI(MVSIZ), CN_LOC(MVSIZ), CE_LOC(MVSIZ),
62 . IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ), SUBTRIA(*), IRTLM(4,*),
63 . IF_ADH(*), IFADHI(MVSIZ)
64 my_real
65 . PENMIN, EPS, PENE_OLD(5,*)
66 my_real
67 . N1(MVSIZ), N2(MVSIZ), N3(MVSIZ),
68 . H1(MVSIZ), H2(MVSIZ), H3(MVSIZ), H4(MVSIZ),
69 . xx(mvsiz,5), yy(mvsiz,5), zz(mvsiz,5),
70 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
71 . nnx(mvsiz,5), nny(mvsiz,5), nnz(mvsiz,5),
72 . lb(mvsiz), lc(mvsiz), gap_nm(4,mvsiz), gaps(mvsiz),
73 . pene(mvsiz), gapmxl(mvsiz), gapv(mvsiz), base_adh(mvsiz)
74 real*4 vtx_bisector(3,2,*)
75 my_real , INTENT(INOUT) :: dist(mvsiz)
76 my_real , INTENT(IN) :: time
77C-----------------------------------------------
78C L o c a l V a r i a b l e s
79C-----------------------------------------------
80 INTEGER I, J, K, L, N, I1, I2, JG, IT, ITRIA(2,4), I3, I4,
81 . IB1, IB2, IB3, IBX, IX, IY, IZ, NBORD, KBORD(MVSIZ)
82 my_real
83 . AAA, NNI, NI2, H0, PENE_SHFT,
84 . NN, NNE(MVSIZ), XH(MVSIZ), YH(MVSIZ), ZH(MVSIZ), XC, YC, ZC, DC, P1, P2, GAPM,
85 . bb(mvsiz), la(mvsiz)
86 my_real
87 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
88 . nn1(mvsiz), nn2(mvsiz), nn3(mvsiz), dd, s,
89 . epseg, ax, bx, cx, gap_mm(mvsiz),
90 . lbs(mvsiz), lcs(mvsiz), xp(mvsiz), yp(mvsiz), zp(mvsiz),
91 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
92 . n1x,n1y,n1z,n1n,
93 . n2x,n2y,n2z,n2n,
94 . linterp, lbm, lcm, lam, xl, yl, zl,
95 . pn, vx, vy, vz,
96 . nne1,nne2,nne4
97 DATA itria/1,2,2,3,3,4,4,1/
98C--------------------------------------------------------
99C zone limite interpolation des normales
100 epseg = (two+half)/hundred
101C--------------------------------------------------------
102 DO i=1,jlt
103C
104 x0n(i,1) = xx(i,1) - xx(i,5)
105 y0n(i,1) = yy(i,1) - yy(i,5)
106 z0n(i,1) = zz(i,1) - zz(i,5)
107C
108 x0n(i,2) = xx(i,2) - xx(i,5)
109 y0n(i,2) = yy(i,2) - yy(i,5)
110 z0n(i,2) = zz(i,2) - zz(i,5)
111C
112 x0n(i,3) = xx(i,3) - xx(i,5)
113 y0n(i,3) = yy(i,3) - yy(i,5)
114 z0n(i,3) = zz(i,3) - zz(i,5)
115C
116 x0n(i,4) = xx(i,4) - xx(i,5)
117 y0n(i,4) = yy(i,4) - yy(i,5)
118 z0n(i,4) = zz(i,4) - zz(i,5)
119C
120 IF(ix3(i)/=ix4(i))THEN
121 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
122 ELSE
123 gap_mm(i)=gap_nm(3,i)
124 END IF
125C
126 ENDDO
127C--------------------------------------------------------
128C normales aux triangles (recalculees ici pour pas les stocker).
129C--------------------------------------------------------
130 DO i=1,jlt
131C
132 it = subtria(i)
133 i1=itria(1,it)
134 i2=itria(2,it)
135C
136 nx(i) = y0n(i,i1)*z0n(i,i2) - z0n(i,i1)*y0n(i,i2)
137 ny(i) = z0n(i,i1)*x0n(i,i2) - x0n(i,i1)*z0n(i,i2)
138 nz(i) = x0n(i,i1)*y0n(i,i2) - y0n(i,i1)*x0n(i,i2)
139C
140 nn=one/max(em30,sqrt(nx(i)*nx(i)+ ny(i)*ny(i)+ nz(i)*nz(i)))
141 nx(i)=nx(i)*nn
142 ny(i)=ny(i)*nn
143 nz(i)=nz(i)*nn
144C
145 ENDDO
146C--------------------------------------------------------
147C cas general
148C--------------------------------------------------------
149 DO i=1,jlt
150C
151 it = subtria(i)
152 i1=itria(1,it)
153 i2=itria(2,it)
154C
155 la(i)=one-lb(i)-lc(i)
156C
157 gapv(i)= min(gaps(i)+la(i)*gap_mm(i)+lb(i)*gap_nm(i1,i)+lc(i)*gap_nm(i2,i),gapmxl(i))
158 bb(i) = (xx(i,5)-xi(i))*nx(i)+(yy(i,5)-yi(i))*ny(i)+(zz(i,5)-zi(i))*nz(i)
159C
160 IF(ix3(i)/=ix4(i))THEN
161
162 IF(bb(i) <= zero)THEN
163
164 xp(i) =la(i)*xx(i,5)+lb(i)*xx(i,i1)+lc(i)*xx(i,i2)
165 yp(i) =la(i)*yy(i,5)+lb(i)*yy(i,i1)+lc(i)*yy(i,i2)
166 zp(i) =la(i)*zz(i,5)+lb(i)*zz(i,i1)+lc(i)*zz(i,i2)
167 nn1(i) =xi(i)-xp(i)
168 nn2(i) =yi(i)-yp(i)
169 nn3(i) =zi(i)-zp(i)
170 dd = sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i))
171 IF(dd > em03) THEN
172 nn = one/dd
173 n1(i) = nn1(i)*nn
174 n2(i) = nn2(i)*nn
175 n3(i) = nn3(i)*nn
176 ELSE
177 n1(i) = nx(i)
178 n2(i) = ny(i)
179 n3(i) = nz(i)
180 END IF
181 pene(i)=max(zero,gapv(i)-dd)
182 dist(i)=dd
183 ELSE
184
185 IF(bb(i) > zero .AND. la(i) < epseg .AND. mvoisn(i,it)/=0)THEN
186C
187C zone limite interpolation des normales
188C
189C rota solides works well ::
190 nn1(i)=nnx(i,it)
191 nn2(i)=nny(i,it)
192 nn3(i)=nnz(i,it)
193C
194 nni = nx(i)*nn1(i) + ny(i)*nn2(i) + nz(i)*nn3(i)
195 ni2 = nn1(i)*nn1(i) + nn2(i)*nn2(i) + nn3(i)*nn3(i)
196 IF(nni < zero .OR. two*nni*nni < ni2)THEN
197c scharp angle bound nodal normal to 45 from segment normal
198 aaa = sqrt(max(zero,ni2-nni*nni)) - nni
199 nn1(i) = nn1(i) + aaa*nx(i)
200 nn2(i) = nn2(i) + aaa*ny(i)
201 nn3(i) = nn3(i) + aaa*nz(i)
202 ENDIF
203 nn = one/
204 . max(em30,sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i)))
205 nn1(i)=nn1(i)*nn
206 nn2(i)=nn2(i)*nn
207 nn3(i)=nn3(i)*nn
208C
209 s = la(i)/epseg
210C
211C continuite de la normale a la frontiere de la zone limite d'interpolation
212 nn1(i)=(one-s)*nn1(i)+s*nx(i)
213 nn2(i)=(one-s)*nn2(i)+s*ny(i)
214 nn3(i)=(one-s)*nn3(i)+s*nz(i)
215 nn = one/
216 . max(em30,sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i)))
217 nn1(i)=nn1(i)*nn
218 nn2(i)=nn2(i)*nn
219 nn3(i)=nn3(i)*nn
220C
221 n1(i) = nn1(i)
222 n2(i) = nn2(i)
223 n3(i) = nn3(i)
224C
225 xp(i) =la(i)*xx(i,5)+lb(i)*xx(i,i1)+lc(i)*xx(i,i2)
226 yp(i) =la(i)*yy(i,5)+lb(i)*yy(i,i1)+lc(i)*yy(i,i2)
227 zp(i) =la(i)*zz(i,5)+lb(i)*zz(i,i1)+lc(i)*zz(i,i2)
228C
229C PENE(I)=MAX(ZERO,GAPV(I)+(XP(I)-XI(I))*N1(I)+(YP(I)-YI(I))*N2(I)+(ZP(I)-ZI(I))*N3(I))
230 pene(i)=max(zero,gapv(i)+bb(i)) ! The distance against the normal estimates the penetration
231C
232 dist(i)=bb(i)
233 ELSE
234C
235C All other cases : normal == normal to the segment
236 n1(i) = nx(i)
237 n2(i) = ny(i)
238 n3(i) = nz(i)
239 pene(i)=max(zero,gapv(i)+bb(i))
240 dist(i)=bb(i)
241 END IF
242
243 END IF
244
245 ELSEIF(ix3(i)==ix4(i))THEN
246
247 IF(bb(i) <= zero)THEN
248
249 xp(i) =la(i)*xx(i,5)+lb(i)*xx(i,i1)+lc(i)*xx(i,i2)
250 yp(i) =la(i)*yy(i,5)+lb(i)*yy(i,i1)+lc(i)*yy(i,i2)
251 zp(i) =la(i)*zz(i,5)+lb(i)*zz(i,i1)+lc(i)*zz(i,i2)
252 nn1(i) =xi(i)-xp(i)
253 nn2(i) =yi(i)-yp(i)
254 nn3(i) =zi(i)-zp(i)
255 dd = sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i))
256 IF(dd > em03) THEN
257 nn = one/dd
258 n1(i) = nn1(i)*nn
259 n2(i) = nn2(i)*nn
260 n3(i) = nn3(i)*nn
261 ELSE
262 n1(i) = nx(i)
263 n2(i) = ny(i)
264 n3(i) = nz(i)
265 END IF
266 pene(i)=max(zero,gapv(i)-dd)
267 dist(i)=dd
268
269 ELSEIF(bb(i) > zero .AND. ((la(i) < epseg .AND. mvoisn(i,1)/=0).OR.
270 . (lb(i) < epseg .AND. mvoisn(i,2)/=0).OR.
271 . (lc(i) < epseg .AND. mvoisn(i,4)/=0)))THEN
272C
273C zone limite interpolation des normales
274 IF(la(i) < epseg .AND. mvoisn(i,1)/=0)THEN
275 aaa=lb(i)+lc(i)
276 ax=zero
277 bx=lb(i)/aaa
278 cx=lc(i)/aaa
279 s = la(i)/epseg
280 ELSEIF(lb(i) < epseg .AND. mvoisn(i,2)/=0)THEN
281 aaa=la(i)+lc(i)
282 ax=la(i)/aaa
283 bx=zero
284 cx=lc(i)/aaa
285 s = lb(i)/epseg
286 ELSEIF(lc(i) < epseg .AND. mvoisn(i,4)/=0)THEN
287 aaa=la(i)+lb(i)
288 ax=la(i)/aaa
289 bx=lb(i)/aaa
290 cx=zero
291 s = lc(i)/epseg
292 END IF
293 nn1(i)=(bx+cx-ax)*nnx(i,i1)+(ax+cx-bx)*nnx(i,i2)+(ax+bx-cx)*nnx(i,5)
294 nn2(i)=(bx+cx-ax)*nny(i,i1)+(ax+cx-bx)*nny(i,i2)+(ax+bx-cx)*nny(i,5)
295 nn3(i)=(bx+cx-ax)*nnz(i,i1)+(ax+cx-bx)*nnz(i,i2)+(ax+bx-cx)*nnz(i,5)
296C
297 nni = nx(i)*nn1(i) + ny(i)*nn2(i) + nz(i)*nn3(i)
298 ni2 = nn1(i)*nn1(i) + nn2(i)*nn2(i) + nn3(i)*nn3(i)
299 IF(nni < zero .OR. two*nni*nni < ni2)THEN
300c scharp angle bound nodal normal to 45 from segment normal
301 aaa = sqrt(max(zero,ni2-nni*nni)) - nni
302 nn1(i) = nn1(i) + aaa*nx(i)
303 nn2(i) = nn2(i) + aaa*ny(i)
304 nn3(i) = nn3(i) + aaa*nz(i)
305 ENDIF
306 nn = one/
307 . max(em30,sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i)))
308 nn1(i)=nn1(i)*nn
309 nn2(i)=nn2(i)*nn
310 nn3(i)=nn3(i)*nn
311C
312C continuite de la normale a la frontiere de la zone limite d'interpolation
313 nn1(i)=(one-s)*nn1(i)+s*nx(i)
314 nn2(i)=(one-s)*nn2(i)+s*ny(i)
315 nn3(i)=(one-s)*nn3(i)+s*nz(i)
316 nn = one/
317 . max(em30,sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i)))
318 nn1(i)=nn1(i)*nn
319 nn2(i)=nn2(i)*nn
320 nn3(i)=nn3(i)*nn
321C
322 n1(i) = nn1(i)
323 n2(i) = nn2(i)
324 n3(i) = nn3(i)
325C
326 xp(i) =la(i)*xx(i,5)+lb(i)*xx(i,i1)+lc(i)*xx(i,i2)
327 yp(i) =la(i)*yy(i,5)+lb(i)*yy(i,i1)+lc(i)*yy(i,i2)
328 zp(i) =la(i)*zz(i,5)+lb(i)*zz(i,i1)+lc(i)*zz(i,i2)
329C
330C PENE(I)=MAX(ZERO,GAPV(I)+(XP(I)-XI(I))*N1(I)+(YP(I)-YI(I))*N2(I)+(ZP(I)-ZI(I))*N3(I))
331 pene(i)=max(zero,gapv(i)+bb(i)) ! The distance against the normal estimates the penetration
332 dist(i)=bb(i)
333C
334 ELSE
335C
336C All other cases : normal == normal to the segment
337 n1(i) = nx(i)
338 n2(i) = ny(i)
339 n3(i) = nz(i)
340 pene(i)=max(zero,gapv(i)+bb(i))
341 dist(i)=bb(i)
342C
343 END IF
344
345C
346 END IF
347
348C-------------------------------------------
349c if(itab(nsvg(i))==2810875.or.
350c . itab(ix1(i))==2810875.or.itab(ix2(i))==2810875.or.itab(ix3(i))==2810875.or.itab(ix4(i))==2810875)
351c . print *,'dst3-avant',itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),mvoisn(i,1:4),
352c . pene(i),pene_old(5,cand_n(i))
353C-------------------------------------------
354c if(itab(nsvg(i))==10105970)
355c . print *,'dst3-avant',itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),mvoisn(i,1:4),
356c . pene(i),pene_old(5,cand_n(i))
357 END DO
358C--------------------------------------------------------
359C cas particulier au bord
360C--------------------------------------------------------
361 nbord = 0
362 DO i=1,jlt
363C
364 it = subtria(i)
365 i1=itria(1,it)
366 i2=itria(2,it)
367C
368 IF(ix3(i)/=ix4(i))THEN
369
370 ib1=ibound(i1,i)
371 ib2=ibound(i2,i)
372C
373C Projection sur le plan et Distance signee au plan
374 xh(i)=xi(i)+bb(i)*nx(i)
375 yh(i)=yi(i)+bb(i)*ny(i)
376 zh(i)=zi(i)+bb(i)*nz(i)
377
378 IF(mvoisn(i,it)==0)THEN
379C
380C Distance signee a l'arete
381C
382C upper skin --------------------
383C |
384C I |
385C -BB: | |
386C | NNE |
387C neutral fiber - - - C - - H < - - >
388C <-------------->|
389C Gap |
390C |
391C --------------------
392C
393 nn1(i)=nnx(i,it)
394 nn2(i)=nny(i,it)
395 nn3(i)=nnz(i,it)
396 nne(i)= (xh(i)-xx(i,i1))*nn1(i)+ (yh(i)-yy(i,i1))*nn2(i)+ (zh(i)-zz(i,i1))*nn3(i)
397C
398 nbord=nbord+1
399 kbord(nbord)=i
400C
401 ELSEIF((ib1/=0.AND.ib2==0).OR.
402 . (ib2/=0.AND.ib1==0))THEN
403C
404 ibx=max(ib1,ib2)
405 IF(ib1/=0)THEN
406 ix =i1
407 ELSE !IF(IB2/=0)THEN
408 ix =i2
409 END IF
410
411 IF(vtx_bisector(1,1,ibx)/=zero.OR.
412 . vtx_bisector(2,1,ibx)/=zero.OR.
413 . vtx_bisector(3,1,ibx)/=zero.OR.
414 . vtx_bisector(1,2,ibx)/=zero.OR.
415 . vtx_bisector(2,2,ibx)/=zero.OR.
416 . vtx_bisector(3,2,ibx)/=zero)THEN
417C
418 p1 = (xh(i)-xx(i,ix))* vtx_bisector(1,1,ibx)+
419 . (yh(i)-yy(i,ix))* vtx_bisector(2,1,ibx)+
420 . (zh(i)-zz(i,ix))* vtx_bisector(3,1,ibx)
421 p2 = (xh(i)-xx(i,ix))* vtx_bisector(1,2,ibx)+
422 . (yh(i)-yy(i,ix))* vtx_bisector(2,2,ibx)+
423 . (zh(i)-zz(i,ix))* vtx_bisector(3,2,ibx)
424C
425 IF(p1 < gaps(i) .AND. p2 < gaps(i))THEN
426C
427 nn1(i)=vtx_bisector(1,1,ibx)+vtx_bisector(1,2,ibx)
428 nn2(i)=vtx_bisector(2,1,ibx)+vtx_bisector(2,2,ibx)
429 nn3(i)=vtx_bisector(3,1,ibx)+vtx_bisector(3,2,ibx)
430C
431 nn=sqrt(nn1(i)*nn1(i)+nn2(i)*nn2(i)+nn3(i)*nn3(i))
432 nn = one/max(em30,nn)
433 nn1(i)=nn1(i)*nn
434 nn2(i)=nn2(i)*nn
435 nn3(i)=nn3(i)*nn
436C
437 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
438C
439 nbord=nbord+1
440 kbord(nbord)=i
441C
442 ELSEIF(p1 < gaps(i))THEN
443
444 nn1(i)= vtx_bisector(1,1,ibx)
445 nn2(i)= vtx_bisector(2,1,ibx)
446 nn3(i)= vtx_bisector(3,1,ibx)
447 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
448
449 nbord=nbord+1
450 kbord(nbord)=i
451
452 ELSEIF(p2 < gaps(i))THEN
453
454 nn1(i)= vtx_bisector(1,2,ibx)
455 nn2(i)= vtx_bisector(2,2,ibx)
456 nn3(i)= vtx_bisector(3,2,ibx)
457 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
458
459 nbord=nbord+1
460 kbord(nbord)=i
461
462 ELSE
463
464c IF(NSVG(I) > 0)THEN
465c print *,' ** possible internal error wrt p1,p2 in i25dst3-3',ib1,ib2,
466c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),p1,p2
467c ELSE
468c print *,' ** possible internal error wrt p1,p2 in i25dst3-3',ib1,ib2,
469c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),p1,p2
470c END IF
471
472 END IF
473C
474 ELSE ! IF(VTX_BISECTOR(1,1,IBX)/=ZERO.OR.
475
476 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
477 vy = y0n(i,ix)
478 vz = z0n(i,ix)
479 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
480 pn = ((xh(i)-xx(i,ix))*vx+(yh(i)-yy(i,ix))*vy+(zh(i)-zz(i,ix))*vz)*nn
481 IF(pn < gaps(i))THEN
482
483 nn1(i)= vx*nn
484 nn2(i)= vy*nn
485 nn3(i)= vz*nn
486 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
487
488 nbord=nbord+1
489 kbord(nbord)=i
490
491 ELSE
492
493c IF(NSVG(I) > 0)THEN
494c print *,' ** possible internal error wrt pn in i25dst3-3',ib1,ib2,
495c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),pn
496c ELSE
497c print *,' ** possible internal error wrt pn in i25dst3-3',ib1,ib2,
498c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),pn
499c END IF
500
501 END IF
502C
503 END IF ! IF(VTX_BISECTOR(1,1,IBX)/=ZERO.OR.
504
505 END IF
506
507 ELSEIF(ix3(i)==ix4(i))THEN
508
509 ib1=ibound(1,i)
510 ib2=ibound(2,i)
511 ib3=ibound(3,i)
512C
513C Projection sur le plan et Distance signee au plan
514 xh(i)=xi(i)+bb(i)*nx(i)
515 yh(i)=yi(i)+bb(i)*ny(i)
516 zh(i)=zi(i)+bb(i)*nz(i)
517
518 IF(mvoisn(i,1)==0.OR.
519 . mvoisn(i,2)==0.OR.
520 . mvoisn(i,4)==0)THEN
521
522 nne(i)=gaps(i)
523 nne1 = (xh(i)-xx(i,i1))*nnx(i,i1)+(yh(i)-yy(i,i1))*nny(i,i1)+(zh(i)-zz(i,i1))*nnz(i,i1)
524 nne2 = (xh(i)-xx(i,i2))*nnx(i,i2)+(yh(i)-yy(i,i2))*nny(i,i2)+(zh(i)-zz(i,i2))*nnz(i,i2)
525 nne4 = (xh(i)-xx(i,5))*nnx(i,5)+(yh(i)-yy(i,5))*nny(i,5)+(zh(i)-zz(i,5))*nnz(i,5)
526
527
528 IF((mvoisn(i,1)==0 .AND. nne1 < nne(i)) .OR.
529 . (mvoisn(i,2)==0 .AND. nne2 < nne(i)) .OR.
530 . (mvoisn(i,4)==0 .AND. nne4 < nne(i))) THEN
531
532 nbord=nbord+1
533 kbord(nbord)=i
534
535 IF(mvoisn(i,1) == 0 .AND. nne1 < nne(i)) THEN
536 nne(i)=nne1
537 nn1(i)=nnx(i,i1)
538 nn2(i)=nny(i,i1)
539 nn3(i)=nnz(i,i1)
540 ENDIF
541 IF(mvoisn(i,2) == 0 .AND. nne2 < nne(i))THEN
542C Distance signee a l'arete
543 nne(i)=nne2
544 nn1(i)=nnx(i,i2)
545 nn2(i)=nny(i,i2)
546 nn3(i)=nnz(i,i2)
547 ENDIF
548 IF(mvoisn(i,4) == 0 .AND. nne4 < nne(i)) THEN
549 nne(i)=nne4
550 nn1(i)=nnx(i,5)
551 nn2(i)=nny(i,5)
552 nn3(i)=nnz(i,5)
553 END IF
554 ENDIF
555C
556 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
557 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
558 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
559C
560 ibx=max(ib1,ib2,ib3)
561 IF(ib1/=0)THEN
562 ix =1
563 iy =2
564 iz =3
565 ELSEIF(ib2/=0)THEN
566 ix =2
567 iy =3
568 iz =1
569 ELSE !IF(IB3/=0)THEN
570 ix =3
571 iy =1
572 iz =2
573 END IF
574C
575 IF(vtx_bisector(1,1,ibx)/=zero.OR.
576 . vtx_bisector(2,1,ibx)/=zero.OR.
577 . vtx_bisector(3,1,ibx)/=zero.OR.
578 . vtx_bisector(1,2,ibx)/=zero.OR.
579 . vtx_bisector(2,2,ibx)/=zero.OR.
580 . vtx_bisector(3,2,ibx)/=zero)THEN
581 p1 = (xh(i)-xx(i,ix))* vtx_bisector(1,1,ibx)+
582 . (yh(i)-yy(i,ix))* vtx_bisector(2,1,ibx)+
583 . (zh(i)-zz(i,ix))* vtx_bisector(3,1,ibx)
584 p2 = (xh(i)-xx(i,ix))* vtx_bisector(1,2,ibx)+
585 . (yh(i)-yy(i,ix))* vtx_bisector(2,2,ibx)+
586 . (zh(i)-zz(i,ix))* vtx_bisector(3,2,ibx)
587
588 IF(p1 < gaps(i) .AND. p2 < gaps(i))THEN
589C
590 nn1(i)=vtx_bisector(1,1,ibx)+vtx_bisector(1,2,ibx)
591 nn2(i)=vtx_bisector(2,1,ibx)+vtx_bisector(2,2,ibx)
592 nn3(i)=vtx_bisector(3,1,ibx)+vtx_bisector(3,2,ibx)
593C
594 nn=sqrt(nn1(i)*nn1(i)+nn2(i)*nn2(i)+nn3(i)*nn3(i))
595 nn = one/max(em30,nn)
596 nn1(i)=nn1(i)*nn
597 nn2(i)=nn2(i)*nn
598 nn3(i)=nn3(i)*nn
599C
600 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
601C
602 nbord=nbord+1
603 kbord(nbord)=i
604C
605 ELSEIF(p1 < gaps(i))THEN
606
607 nn1(i)= vtx_bisector(1,1,ibx)
608 nn2(i)= vtx_bisector(2,1,ibx)
609 nn3(i)= vtx_bisector(3,1,ibx)
610 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
611
612 nbord=nbord+1
613 kbord(nbord)=i
614
615 ELSEIF(p2 < gaps(i))THEN
616
617 nn1(i)= vtx_bisector(1,2,ibx)
618 nn2(i)= vtx_bisector(2,2,ibx)
619 nn3(i)= vtx_bisector(3,2,ibx)
620 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
621
622 nbord=nbord+1
623 kbord(nbord)=i
624
625 ELSE
626
627c IF(NSVG(I) > 0)THEN
628c print *,' ** possible internal error wrt p1,p2 in i25dst3-3',ib1,ib2,
629c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),p1,p2
630c ELSE
631c print *,' ** possible internal error wrt p1,p2 in i25dst3-3',ib1,ib2,
632c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),p1,p2
633c END IF
634
635 END IF
636
637 ELSE ! IF(VTX_BISECTOR(1,1,IBX)/=ZERO.OR.
638
639 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
640 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
641 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
642 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
643 pn = ((xh(i)-xx(i,ix))*vx+(yh(i)-yy(i,ix))*vy+(zh(i)-zz(i,ix))*vz)*nn
644
645 IF(pn < gaps(i))THEN
646
647 nn1(i)= vx*nn
648 nn2(i)= vy*nn
649 nn3(i)= vz*nn
650 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
651
652 nbord=nbord+1
653 kbord(nbord)=i
654
655 ELSE
656
657c IF(NSVG(I) > 0)THEN
658c print *,' ** possible internal error wrt pn in i25dst3-3',ib1,ib2,
659c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),pn
660c ELSE
661c print *,' ** possible internal error wrt pn in i25dst3-3',ib1,ib2,
662c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),pn
663c END IF
664
665 END IF
666
667 END IF ! IF(VTX_BISECTOR(1,1,IBX)/=ZERO.OR.
668
669 END IF
670C
671 END IF
672
673 END DO
674C-------------------------------------------
675 IF(isharp==3)THEN
676
677 DO k=1,nbord
678 i=kbord(k)
679
680 IF(nne(i) > zero)THEN
681 pene(i)=zero
682 END IF
683
684 END DO
685 ELSEIF(isharp==1)THEN
686
687 DO k=1,nbord
688 i=kbord(k)
689
690 gapm = gapv(i)-gaps(i)
691 IF(nne(i) > zero .AND. bb(i)+gapm < zero)THEN
692C
693C ___________________________ xxx
694C |xxxxxx
695C GapS |xxxxxxxx <=> Upper Round Corner
696C ___________________________|xxxxxxxx
697C | |
698C GapM | | <=> Straight shape of the edge
699C ---------------------------M-------
700C <------>
701C GapS
702C
703
704 xc=xp(i)+gapm*nx(i)
705 yc=yp(i)+gapm*ny(i)
706 zc=zp(i)+gapm*nz(i)
707 n1(i) =xi(i)-xc
708 n2(i) =yi(i)-yc
709 n3(i) =zi(i)-zc
710 dc = sqrt(n1(i)*n1(i)+ n2(i)*n2(i)+ n3(i)*n3(i))
711
712 IF(dc > em04) THEN
713 nn = one/dc
714 n1(i) = n1(i)*nn
715 n2(i) = n2(i)*nn
716 n3(i) = n3(i)*nn
717 pene(i)=max(zero,gaps(i)-dc)
718 ELSE
719C
720 nne(i)=nne(i)-gaps(i)
721 IF(-bb(i) < gapv(i)+nne(i))THEN ! Test 45 degres !
722C
723C Normale horizonale
724 n1(i) = nn1(i)
725 n2(i) = nn2(i)
726 n3(i) = nn3(i)
727C
728 pene(i)=max(zero,-nne(i))
729C
730 ELSE
731C
732C Normale verticale !
733 n1(i) = nx(i)
734 n2(i) = ny(i)
735 n3(i) = nz(i)
736C
737 pene(i)=max(zero,gapv(i)+bb(i))
738 END IF
739 END IF
740
741 ELSE !IF(NNE(I) > ZERO .AND. BB(I)+GAPM < ZERO)THEN
742
743 nne(i)=nne(i)-gaps(i)
744 IF(nne(i) >= zero)THEN
745C
746C Must remain a roundoff error <=> NNE ~ 0
747c IF(NSVG(I) > 0)THEN
748c print *,' ** possible internal error in i25dst3-3, nne should be negative',
749c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I),pene(i)
750c ELSE
751c print *,' ** possible internal error in i25dst3-3, nne should be negative',
752c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I),pene(i)
753c END IF
754
755 pene(i)=zero
756 cycle
757
758 END IF
759C
760 IF(gapv(i)+nne(i) > zero)THEN
761
762 IF(-bb(i) < gapv(i)+nne(i))THEN ! Test 45 degres !
763C
764C Normale horizonale
765 n1(i) = nn1(i)
766 n2(i) = nn2(i)
767 n3(i) = nn3(i)
768C
769 pene(i)=-nne(i)
770C
771 ELSE
772C
773C Normale verticale !
774 n1(i) = nx(i)
775 n2(i) = ny(i)
776 n3(i) = nz(i)
777C
778 pene(i)=max(zero,gapv(i)+bb(i))
779 END IF
780
781 END IF
782
783 END IF
784C-------------------------------------------
785c if(itab(nsvg(i))==10105970)
786c . print *,'dst3-apres',itab(nsvg(i)),kbord,pene(i),pene_old(5,cand_n(i)),n1(i),n2(i),n3(i),
787c . nne(I),bb(I),gapv(i)
788C-------------------------------------------
789 END DO
790
791 ELSEIF(isharp==2)THEN
792
793 DO k=1,nbord
794 i=kbord(k)
795
796 nne(i)=nne(i)-gaps(i)
797 IF(nne(i) >= zero)THEN
798C
799C Must remain a roundoff error <=> NNE ~ 0
800c IF(NSVG(I) > 0)THEN
801c print *,' ** possible internal error in i25dst3-3, nne should be negative',
802c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I)
803c ELSE
804c print *,' ** possible internal error in i25dst3-3, nne should be negative',
805c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I)
806c END IF
807
808 cycle
809
810 END IF
811C
812 IF(bb(i) > zero)THEN
813
814 IF(-bb(i) < gapv(i)+nne(i))THEN ! Test 45 degres !
815C
816C Normale horizonale
817 n1(i) = nn1(i)
818 n2(i) = nn2(i)
819 n3(i) = nn3(i)
820C
821 pene(i)=-nne(i)
822C
823 ELSE
824C
825C Normale verticale !
826 n1(i) = nx(i)
827 n2(i) = ny(i)
828 n3(i) = nz(i)
829C
830 pene(i)=max(zero,gapv(i)+bb(i))
831 END IF
832
833 ELSEIF(gapv(i)+nne(i) > zero)THEN
834C
835C Normale radiale shiftee
836 xc =xh(i)-(gapv(i)+nne(i))*nn1(i)
837 yc =yh(i)-(gapv(i)+nne(i))*nn2(i)
838 zc =zh(i)-(gapv(i)+nne(i))*nn3(i)
839 nn1(i) =xi(i)-xc
840 nn2(i) =yi(i)-yc
841 nn3(i) =zi(i)-zc
842 dc = sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i))
843 IF(dc > em04) THEN
844 nn = one/dc
845 n1(i) = nn1(i)*nn
846 n2(i) = nn2(i)*nn
847 n3(i) = nn3(i)*nn
848 pene(i)=max(zero,gapv(i)-dc)
849 ELSE
850C
851C Keep what was done in the general case : Normale ~ verticale.
852 END IF
853
854 END IF
855C-------------------------------------------
856c if(itab(nsvg(i))==10105970)
857c . print *,'dst3-apres',itab(nsvg(i)),kbord,pene(i),pene_old(5,cand_n(i)),n1(i),n2(i),n3(i),
858c . nne(I),bb(i),gapv(i)
859C-------------------------------------------
860 ENDDO
861 END IF
862C--------------------------------------------------------
863C Calculer PENE et Hi
864C--------------------------------------------------------
865 DO i=1,jlt
866C
867 ce_loc(i) = cand_e(i)
868 cn_loc(i) = cand_n(i)
869C
870 it = subtria(i)
871C
872 IF(ix3(i)/=ix4(i))THEN
873C
874 h0 = fourth * la(i)
875 IF(it==1)THEN
876 h1(i) = h0 + lb(i)
877 h2(i) = h0 + lc(i)
878 h3(i) = h0
879 h4(i) = h0
880 ELSEIF(it==2)THEN
881 h1(i) = h0
882 h2(i) = h0 + lb(i)
883 h3(i) = h0 + lc(i)
884 h4(i) = h0
885 ELSEIF(it==3)THEN
886 h1(i) = h0
887 h2(i) = h0
888 h3(i) = h0 + lb(i)
889 h4(i) = h0 + lc(i)
890 ELSEIF(it==4)THEN
891 h1(i) = h0 + lc(i)
892 h2(i) = h0
893 h3(i) = h0
894 h4(i) = h0 + lb(i)
895 END IF
896C
897 ELSE
898C
899 h1(i) = lb(i)
900 h2(i) = lc(i)
901 h3(i) = la(i)
902 h4(i) = zero
903C
904 END IF
905C
906 END DO
907C
908C-----------------------------------------------------------------------
909C PENE_OLD(5) <=> Initial Penetration
910C-----------------------------------------------------------------------
911 IF (time==zero) THEN
912 IF (inacti==5.AND.ivis2/=-1) THEN
913 DO i=1,jlt
914 IF (pene(i) == zero) cycle
915 jg = nsvg(i)
916 n = cand_n(i)
917 IF(jg > 0)THEN
918 IF(irtlm(1,n) > 0)THEN ! initial pen
919 pene_old(5,n)=max(pene(i),pene_old(5,n))
920 END IF
921 ELSE
922 jg = -jg
923 IF(irtlm_fi(nin)%P(1,jg) > 0)THEN
924 pene_oldfi(nin)%P(5,jg)=max(pene(i),pene_oldfi(nin)%P(5,jg))
925 END IF
926 END IF
927 END DO
928 END if!(INACTI==5.AND.IVIS2/=-1) THEN
929 END IF
930 IF(ivis2/=-1) THEN
931C IF(INACTI==5)THEN
932 DO i=1,jlt
933C
934 IF(pene(i) == zero)cycle
935C
936 jg = nsvg(i)
937 n = cand_n(i)
938 IF(jg > 0)THEN
939 IF(irtlm(1,n) < 0)THEN ! new impact
940 pene_old(5,n)=pene(i)
941 END IF
942 pene_shft = max(zero,pene(i)-pene_old(5,n))
943 ELSE
944 jg = -jg
945 IF(irtlm_fi(nin)%P(1,jg) < 0)THEN ! new impact
946 pene_oldfi(nin)%P(5,jg)=pene(i)
947 END IF
948 pene_shft= max(zero,pene(i)-pene_oldfi(nin)%P(5,jg))
949
950 ENDIF
951 pene(i) = pene_shft
952 ENDDO
953C END IF !IF (INACTI==5)
954C pene_old(5,N) is defined from the real_gap, gap = 2*real_gap, pene is measured from real_gap
955 ELSE ! IVIS2==-1 (Adhesion case)
956 DO i=1,jlt
957 IF(pene(i) == zero)cycle
958 jg = nsvg(i)
959 n = cand_n(i)
960 IF(jg > 0)THEN
961 IF(irtlm(1,n) < 0)THEN ! new impact
962 pene_old(5,n)= max(zero,pene(i)-half*gapv(i))
963 ENDIF
964 IF(pene(i)>=half*gapv(i)) if_adh(n) = 1
965 base_adh(i) = pene_old(5,n) + half*gapv(i)
966 ifadhi(i) = if_adh(n)
967 ELSE
968 jg = -jg
969 IF(irtlm_fi(nin)%P(1,jg) < 0)THEN
970 pene_oldfi(nin)%P(5,jg)= max(zero,pene(i)-half*gapv(i))
971 ENDIF
972 IF(pene(i)>=half*gapv(i)) if_adhfi(nin)%P(jg) = 1
973 base_adh(i) = pene_oldfi(nin)%P(5,jg) + half*gapv(i)
974 ifadhi(i) = if_adhfi(nin)%P(jg)
975 ENDIF
976 ENDDO
977 ENDIF
978C-----------------------------------------------------------------------
979 RETURN
980 END
if
if(complex_arithmetic) id
i25dst3_3
subroutine i25dst3_3(jlt, cand_n, cand_e, cn_loc, ce_loc, irtlm, xx, yy, zz, gap_nm, xi, yi, zi, gaps, gapmxl, isharp, nnx, nny, nnz, n1, n2, n3, h1, h2, h3, h4, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, inacti, kini, itab, lb, lc, penmin, eps, pene, pene_old, subtria, gapv, ivis2, if_adh, ifadhi, base_adh, mvoisn, ibound, vtx_bisector, dist, time)
Definition i25dst3_3.F:42
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
tri7box::irtlm_fi
type(int_pointer2), dimension(:), allocatable irtlm_fi
Definition tri7box.F:533
tri7box::pene_oldfi
type(real_pointer2), dimension(:), allocatable pene_oldfi
Definition tri7box.F:544