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

Go to the source code of this file.

Functions/Subroutines

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)

Function/Subroutine Documentation

◆ i25dst3_3()

subroutine i25dst3_3 ( integer jlt,
integer, dimension(*) cand_n,
integer, dimension(*) cand_e,
integer, dimension(mvsiz) cn_loc,
integer, dimension(mvsiz) ce_loc,
integer, dimension(4,*) irtlm,
xx,
yy,
zz,
gap_nm,
xi,
yi,
zi,
gaps,
gapmxl,
integer isharp,
nnx,
nny,
nnz,
n1,
n2,
n3,
h1,
h2,
h3,
h4,
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(mvsiz) kini,
integer, dimension(*) itab,
lb,
lc,
penmin,
eps,
pene,
pene_old,
integer, dimension(*) subtria,
gapv,
integer ivis2,
integer, dimension(*) if_adh,
integer, dimension(mvsiz) ifadhi,
base_adh,
integer, dimension(mvsiz,4) mvoisn,
integer, dimension(4,*) ibound,
real*4, dimension(3,2,*) vtx_bisector,
dimension(mvsiz), intent(inout) dist,
intent(in) time )

Definition at line 30 of file i25dst3_3.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, 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
my_real
#define my_real
Definition cppsort.cpp:32
if
if(complex_arithmetic) id
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
tri7box::if_adhfi
type(int_pointer), dimension(:), allocatable if_adhfi
Definition tri7box.F:440