OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i25pen3.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/.
23C
24!||====================================================================
25!|| i25pen3 ../starter/source/interfaces/inter3d1/i25pen3.F
26!||--- called by ------------------------------------------------------
27!|| inint3 ../starter/source/interfaces/inter3d1/inint3.F
28!||--- uses -----------------------------------------------------
29!|| message_mod ../starter/share/message_module/message_mod.F
30!||====================================================================
31 SUBROUTINE i25pen3(
32 1 JLT ,CAND_N ,CAND_E ,PENMIN ,PENMAX,
33 2 X1 ,X2 ,X3 ,X4 ,X0 ,
34 3 Y1 ,Y2 ,Y3 ,Y4 ,Y0 ,
35 4 Z1 ,Z2 ,Z3 ,Z4 ,Z0 ,
36 5 XI ,YI ,ZI ,NSN ,IX1 ,
37 6 IX2 ,IX3 ,IX4 ,NSVG ,NRTM ,
38 7 MSEGLO ,GAPS ,IRECT ,IRTLM ,
39 9 TIME_S ,PENE_OLD,ITAB ,MSEGTYP,ISHARP,
40 A NNX ,NNY ,NNZ ,GAP_NM ,MVOISN,
41 B GAPMXL ,IVIS2 ,IBOUND,VTX_BISECTOR,ILEV,
42 C INACTI)
43C
44 USE message_mod
45C
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"
54C-----------------------------------------------
55C C o m m o n B l o c k s
56C-----------------------------------------------
57#include "scr05_c.inc"
58C-----------------------------------------------
59 INTEGER JLT, NSVG(*), NSN, MSEGLO(*), MSEGTYP(*), NRTM, IVIS2, ISHARP
60 INTEGER , INTENT(IN) :: INACTI, ILEV
61 my_real
62 . GAPS(*), GAP_NM(4,*), GAPMXL(*)
63 INTEGER IRECT(4,*), ITAB(*),CAND_E(*),CAND_N(*),IRTLM(4,*)
64 INTEGER IX1(MVSIZ),IX2(MVSIZ),IX3(MVSIZ),IX4(MVSIZ),
65 . MVOISN(MVSIZ,4), IBOUND(4,MVSIZ)
66 my_real
67 . PENMAX,PENMIN,PENE_OLD(5,NSN),TIME_S(NSN)
68 my_real
69 . X1(MVSIZ), X2(MVSIZ), X3(MVSIZ), X4(MVSIZ), X0(MVSIZ),
70 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz), y0(mvsiz),
71 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz), z0(mvsiz),
72 . xi(mvsiz), yi(mvsiz), zi(mvsiz),
73 . nnx(mvsiz,5), nny(mvsiz,5), nnz(mvsiz,5)
74C-----------------------------------------------
75C L o c a l V a r i a b l e s
76C-----------------------------------------------
77 INTEGER I, J, K, L, N, IT, N4N, MGLOB, ETYP(MVSIZ),
78 . NINDX, INDX(MVSIZ), I4N(MVSIZ),
79 . far(mvsiz,4), subtria(mvsiz), i1, i2, itria(2,4), jj
80 INTEGER IB1, IB2, IB3, IBX, IX, IY, IZ, NBORD, KBORD(MVSIZ)
81 my_real
82 . AAA,BB(MVSIZ,4)
83 my_real
84 . XIJ(MVSIZ,4), XI0V(MVSIZ),
85 . YIJ(MVSIZ,4), YI0V(MVSIZ),
86 . ZIJ(MVSIZ,4), ZI0V(MVSIZ),
87 . DX, DY, DZ, DD(MVSIZ,4),
88 . xp, yp, zp,
89 . gap_mm(mvsiz), gapv(mvsiz,4), gap2,norm
90 my_real
91 . pene,
92 . sx1,sx2,sx3,sx4,
93 . sy1,sy2,sy3,sy4,
94 . sz1,sz2,sz3,sz4,
95 . x0n(mvsiz,4), y0n(mvsiz,4), z0n(mvsiz,4),
96 . xn(mvsiz,4),yn(mvsiz,4),zn(mvsiz,4),
97 . side(mvsiz,4),
98 . x12(mvsiz), y12(mvsiz), z12(mvsiz),
99 . px, py, pz, pp, ll, nn, pn, vx, vy, vz,
100 . xx(mvsiz,5),yy(mvsiz,5),zz(mvsiz,5),
101 . pent(mvsiz,4), dist(mvsiz), lb(mvsiz,4), lc(mvsiz,4),
102 . lap, lbp(mvsiz,4), lcp(mvsiz,4), hla, hlb(mvsiz,4), hlc(mvsiz,4),
103 . al(mvsiz,4),
104 . nn1(mvsiz), nn2(mvsiz), nn3(mvsiz),
105 . lx, lax, lbx, lcx, ld(mvsiz), dmin,
106 . nne(mvsiz), xh(mvsiz), yh(mvsiz), zh(mvsiz), xc, yc, zc, dc, p1, p2, d1, d2, d3, gapm
107 real*4 vtx_bisector(3,2,*)
108 DATA itria/1,2,2,3,3,4,4,1/
109C-----------------------------------------------
110 nindx=0
111 DO i=1,jlt
112 l = cand_e(i)
113 etyp(i) =msegtyp(l)
114 IF(etyp(i)==0)THEN
115C
116C Computing initial penetrations vs solids
117 nindx=nindx+1
118 indx(nindx)=i
119 ELSEIF((etyp(i) /= 0 .AND. iabs(etyp(i)) <= nrtm) .OR. etyp(i) > nrtm)THEN
120C
121C Computing initial penetrations vs shells except opposite segments to the coating shells
122 nindx=nindx+1
123 indx(nindx)=i
124 END IF
125 END DO
126C-----------------------------------------------
127 IF (iresp==1.AND.penmin<=em06) penmin = two*em06
128C--------
129 far(1:jlt,1:4) = 0
130 dd(1:jlt,1:4) = ep20
131 dist(1:jlt) = ep20
132 ld(1:jlt) = ep20
133C--------
134 DO k=1,nindx
135 i=indx(k)
136
137 pent(i,1:4)=zero
138C
139 xx(i,1) = x1(i)
140 yy(i,1) = y1(i)
141 zz(i,1) = z1(i)
142C
143 xx(i,2) = x2(i)
144 yy(i,2) = y2(i)
145 zz(i,2) = z2(i)
146C
147 xx(i,3) = x3(i)
148 yy(i,3) = y3(i)
149 zz(i,3) = z3(i)
150C
151 xx(i,4) = x4(i)
152 yy(i,4) = y4(i)
153 zz(i,4) = z4(i)
154C
155 xx(i,5) = x0(i)
156 yy(i,5) = y0(i)
157 zz(i,5) = z0(i)
158
159 ENDDO
160C--------------------------------------------------------
161 DO k=1,nindx
162 i=indx(k)
163C
164 x0n(i,1) = xx(i,1) - xx(i,5)
165 y0n(i,1) = yy(i,1) - yy(i,5)
166 z0n(i,1) = zz(i,1) - zz(i,5)
167C
168 x0n(i,2) = xx(i,2) - xx(i,5)
169 y0n(i,2) = yy(i,2) - yy(i,5)
170 z0n(i,2) = zz(i,2) - zz(i,5)
171C
172 x0n(i,3) = xx(i,3) - xx(i,5)
173 y0n(i,3) = yy(i,3) - yy(i,5)
174 z0n(i,3) = zz(i,3) - zz(i,5)
175C
176 x0n(i,4) = xx(i,4) - xx(i,5)
177 y0n(i,4) = yy(i,4) - yy(i,5)
178 z0n(i,4) = zz(i,4) - zz(i,5)
179C
180 IF(ix3(i)/=ix4(i))THEN
181 gap_mm(i)=fourth*(gap_nm(1,i)+gap_nm(2,i)+gap_nm(3,i)+gap_nm(4,i))
182 ELSE
183 gap_mm(i)=gap_nm(3,i)
184 END IF
185C
186 ENDDO
187C--------------------------------------------------------
188 DO k=1,nindx
189 i=indx(k)
190C
191 xi0v(i) = xx(i,5) - xi(i)
192 yi0v(i) = yy(i,5) - yi(i)
193 zi0v(i) = zz(i,5) - zi(i)
194C
195 xij(i,1) = xx(i,1) - xi(i)
196 yij(i,1) = yy(i,1) - yi(i)
197 zij(i,1) = zz(i,1) - zi(i)
198C
199 xij(i,2) = xx(i,2) - xi(i)
200 yij(i,2) = yy(i,2) - yi(i)
201 zij(i,2) = zz(i,2) - zi(i)
202C
203 xij(i,3) = xx(i,3) - xi(i)
204 yij(i,3) = yy(i,3) - yi(i)
205 zij(i,3) = zz(i,3) - zi(i)
206C
207 xij(i,4) = xx(i,4) - xi(i)
208 yij(i,4) = yy(i,4) - yi(i)
209 zij(i,4) = zz(i,4) - zi(i)
210C
211 sx1 = yi0v(i)*zij(i,1) - zi0v(i)*yij(i,1)
212 sy1 = zi0v(i)*xij(i,1) - xi0v(i)*zij(i,1)
213 sz1 = xi0v(i)*yij(i,1) - yi0v(i)*xij(i,1)
214C
215 sx2 = yi0v(i)*zij(i,2) - zi0v(i)*yij(i,2)
216 sy2 = zi0v(i)*xij(i,2) - xi0v(i)*zij(i,2)
217 sz2 = xi0v(i)*yij(i,2) - yi0v(i)*xij(i,2)
218C
219 xn(i,1) = y0n(i,1)*z0n(i,2) - z0n(i,1)*y0n(i,2)
220 yn(i,1) = z0n(i,1)*x0n(i,2) - x0n(i,1)*z0n(i,2)
221 zn(i,1) = x0n(i,1)*y0n(i,2) - y0n(i,1)*x0n(i,2)
222 nn = one/
223 . max(em30,sqrt(xn(i,1)*xn(i,1)+ yn(i,1)*yn(i,1)+ zn(i,1)*zn(i,1)))
224 xn(i,1)=xn(i,1)*nn
225 yn(i,1)=yn(i,1)*nn
226 zn(i,1)=zn(i,1)*nn
227 lb(i,1) = -(xn(i,1)*sx2 + yn(i,1)*sy2 + zn(i,1)*sz2)*nn
228 lc(i,1) = (xn(i,1)*sx1 + yn(i,1)*sy1 + zn(i,1)*sz1)*nn
229C
230 sx3 = yi0v(i)*zij(i,3) - zi0v(i)*yij(i,3)
231 sy3 = zi0v(i)*xij(i,3) - xi0v(i)*zij(i,3)
232 sz3 = xi0v(i)*yij(i,3) - yi0v(i)*xij(i,3)
233C
234 xn(i,2) = y0n(i,2)*z0n(i,3) - z0n(i,2)*y0n(i,3)
235 yn(i,2) = z0n(i,2)*x0n(i,3) - x0n(i,2)*z0n(i,3)
236 zn(i,2) = x0n(i,2)*y0n(i,3) - y0n(i,2)*x0n(i,3)
237 nn = one/
238 . max(em30,sqrt(xn(i,2)*xn(i,2)+ yn(i,2)*yn(i,2)+ zn(i,2)*zn(i,2)))
239 xn(i,2)=xn(i,2)*nn
240 yn(i,2)=yn(i,2)*nn
241 zn(i,2)=zn(i,2)*nn
242 lb(i,2) = -(xn(i,2)*sx3 + yn(i,2)*sy3 + zn(i,2)*sz3)*nn
243 lc(i,2) = (xn(i,2)*sx2 + yn(i,2)*sy2 + zn(i,2)*sz2)*nn
244C
245 sx4 = yi0v(i)*zij(i,4) - zi0v(i)*yij(i,4)
246 sy4 = zi0v(i)*xij(i,4) - xi0v(i)*zij(i,4)
247 sz4 = xi0v(i)*yij(i,4) - yi0v(i)*xij(i,4)
248C
249 xn(i,3) = y0n(i,3)*z0n(i,4) - z0n(i,3)*y0n(i,4)
250 yn(i,3) = z0n(i,3)*x0n(i,4) - x0n(i,3)*z0n(i,4)
251 zn(i,3) = x0n(i,3)*y0n(i,4) - y0n(i,3)*x0n(i,4)
252 nn = one/
253 . max(em30,sqrt(xn(i,3)*xn(i,3)+ yn(i,3)*yn(i,3)+ zn(i,3)*zn(i,3)))
254 xn(i,3)=xn(i,3)*nn
255 yn(i,3)=yn(i,3)*nn
256 zn(i,3)=zn(i,3)*nn
257 lb(i,3) = -(xn(i,3)*sx4 + yn(i,3)*sy4 + zn(i,3)*sz4)*nn
258 lc(i,3) = (xn(i,3)*sx3 + yn(i,3)*sy3 + zn(i,3)*sz3)*nn
259C
260 xn(i,4) = y0n(i,4)*z0n(i,1) - z0n(i,4)*y0n(i,1)
261 yn(i,4) = z0n(i,4)*x0n(i,1) - x0n(i,4)*z0n(i,1)
262 zn(i,4) = x0n(i,4)*y0n(i,1) - y0n(i,4)*x0n(i,1)
263 nn = one/
264 . max(em30,sqrt(xn(i,4)*xn(i,4)+ yn(i,4)*yn(i,4)+ zn(i,4)*zn(i,4)))
265 xn(i,4)=xn(i,4)*nn
266 yn(i,4)=yn(i,4)*nn
267 zn(i,4)=zn(i,4)*nn
268 lb(i,4) = -(xn(i,4)*sx1 + yn(i,4)*sy1 + zn(i,4)*sz1)*nn
269 lc(i,4) = (xn(i,4)*sx4 + yn(i,4)*sy4 + zn(i,4)*sz4)*nn
270C
271 END DO
272C--------------------------------------------------------
273 DO k=1,nindx
274 i=indx(k)
275C
276 aaa = one/max(em30,x0n(i,1)*x0n(i,1)+y0n(i,1)*y0n(i,1)+z0n(i,1)*z0n(i,1))
277 hlc(i,1) = lc(i,1)*abs(lc(i,1))*aaa
278 hlb(i,4) = lb(i,4)*abs(lb(i,4))*aaa
279 al(i,1) = -(xi0v(i)*x0n(i,1)+yi0v(i)*y0n(i,1)+zi0v(i)*z0n(i,1))*aaa
280 al(i,1) = max(zero,min(one,al(i,1)))
281 aaa = one/max(em30,x0n(i,2)*x0n(i,2)+y0n(i,2)*y0n(i,2)+z0n(i,2)*z0n(i,2))
282 hlc(i,2) = lc(i,2)*abs(lc(i,2))*aaa
283 hlb(i,1) = lb(i,1)*abs(lb(i,1))*aaa
284 al(i,2) = -(xi0v(i)*x0n(i,2)+yi0v(i)*y0n(i,2)+zi0v(i)*z0n(i,2))*aaa
285 al(i,2) = max(zero,min(one,al(i,2)))
286 aaa = one/max(em30,x0n(i,3)*x0n(i,3)+y0n(i,3)*y0n(i,3)+z0n(i,3)*z0n(i,3))
287 hlc(i,3) = lc(i,3)*abs(lc(i,3))*aaa
288 hlb(i,2) = lb(i,2)*abs(lb(i,2))*aaa
289 al(i,3) = -(xi0v(i)*x0n(i,3)+yi0v(i)*y0n(i,3)+zi0v(i)*z0n(i,3))*aaa
290 al(i,3) = max(zero,min(one,al(i,3)))
291 aaa = one/max(em30,x0n(i,4)*x0n(i,4)+y0n(i,4)*y0n(i,4)+z0n(i,4)*z0n(i,4))
292 hlc(i,4) = lc(i,4)*abs(lc(i,4))*aaa
293 hlb(i,3) = lb(i,3)*abs(lb(i,3))*aaa
294 al(i,4) = -(xi0v(i)*x0n(i,4)+yi0v(i)*y0n(i,4)+zi0v(i)*z0n(i,4))*aaa
295 al(i,4) = max(zero,min(one,al(i,4)))
296C
297 END DO
298C--------------------------------------------------------
299 it=1
300 i1=itria(1,it)
301 i2=itria(2,it)
302 DO k=1,nindx
303 i=indx(k)
304C
305 x12(i) = xx(i,i2) - xx(i,i1)
306 y12(i) = yy(i,i2) - yy(i,i1)
307 z12(i) = zz(i,i2) - zz(i,i1)
308C
309 lbp(i,it) = lb(i,it)
310 lcp(i,it) = lc(i,it)
311 lap = one-lbp(i,it)-lcp(i,it)
312C HLA, HLB, HLC necessaires pour triangle angle obtu
313 aaa = one / max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
314 hla= lap*abs(lap) * aaa
315 IF(lap<zero.AND.
316 + hla<=hlb(i,it).AND.hla<=hlc(i,it))THEN
317 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
318 lbp(i,it) = max(zero,min(one,lbp(i,it)))
319 lcp(i,it) = one - lbp(i,it)
320 ELSEIF(lbp(i,it)<zero.AND.
321 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)THEN
322 lbp(i,it) = zero
323 lcp(i,it) = al(i,i2)
324 ELSEIF(lcp(i,it)<zero.AND.
325 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))THEN
326 lcp(i,it) = zero
327 lbp(i,it) = al(i,i1)
328 ENDIF
329
330 lap = one-lbp(i,it)-lcp(i,it)
331 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
332 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
333 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
334 dx =xi(i)-xp
335 dy =yi(i)-yp
336 dz =zi(i)-zp
337 dd(i,it)=dx*dx+dy*dy+dz*dz
338
339 END DO
340C--------------------------------------------------------
341 DO k=1,nindx
342 i=indx(k)
343C
344 lap = one-lbp(i,it)-lcp(i,it)
345C
346 gapv(i,it) = min(gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i),gapmxl(i))
347 gap2 = gapv(i,it)**2
348C
349 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))
350
351 IF(etyp(i)==0)THEN
352 IF(bb(i,it) > zero)THEN
353 pent(i,it)=gapv(i,it)+bb(i,it)
354 ELSE
355 pent(i,it)=zero
356 END IF
357 ELSEIF(etyp(i) > nrtm)THEN ! Coating shell (same orientation than solid)
358 IF(bb(i,it) > zero)THEN
359 pent(i,it)=gapv(i,it)+bb(i,it)
360 ELSE
361 pent(i,it)=max(zero,gapv(i,it)-sqrt(dd(i,it)))
362 END IF
363 ELSE
364 IF(bb(i,it) > zero)THEN
365C
366C penetration approximee (cf zone limite interpolation des normales)
367 pent(i,it)=zero
368 ELSE
369 pent(i,it)=max(zero,gapv(i,it)-sqrt(dd(i,it)))
370 END IF
371 END IF
372C
373 ENDDO
374C--------------------------------------------------------
375 n4n=0
376 DO k=1,nindx
377 i=indx(k)
378C
379 IF(ix3(i)/=ix4(i))THEN
380 n4n = n4n+1
381 i4n(n4n)=i
382 ENDIF
383 ENDDO
384C--------------------------------------------------------
385 DO it=2,4
386
387 i1=itria(1,it)
388 i2=itria(2,it)
389
390 DO k=1,n4n
391 i=i4n(k)
392C
393 x12(i) = xx(i,i2) - xx(i,i1)
394 y12(i) = yy(i,i2) - yy(i,i1)
395 z12(i) = zz(i,i2) - zz(i,i1)
396C
397 lbp(i,it) = lb(i,it)
398 lcp(i,it) = lc(i,it)
399 lap = one-lbp(i,it)-lcp(i,it)
400C HLA, HLB, HLC necessaires pour triangle angle obtu
401 aaa = one / max(em20,x12(i)*x12(i)+y12(i)*y12(i)+z12(i)*z12(i))
402 hla= lap*abs(lap) * aaa
403 IF(lap<zero.AND.
404 + hla<=hlb(i,it).AND.hla<=hlc(i,it))THEN
405 lbp(i,it) = (xij(i,i2)*x12(i)+yij(i,i2)*y12(i)+zij(i,i2)*z12(i)) * aaa
406 lbp(i,it) = max(zero,min(one,lbp(i,it)))
407 lcp(i,it) = one - lbp(i,it)
408 ELSEIF(lbp(i,it)<zero.AND.
409 + hlb(i,it)<=hlc(i,it).AND.hlb(i,it)<=hla)THEN
410 lbp(i,it) = zero
411 lcp(i,it) = al(i,i2)
412 ELSEIF(lcp(i,it)<zero.AND.
413 + hlc(i,it)<=hla.AND.hlc(i,it)<=hlb(i,it))THEN
414 lcp(i,it) = zero
415 lbp(i,it) = al(i,i1)
416 ENDIF
417
418 lap = one-lbp(i,it)-lcp(i,it)
419 xp =lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
420 yp =lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
421 zp =lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
422 dx =xi(i)-xp
423 dy =yi(i)-yp
424 dz =zi(i)-zp
425 dd(i,it)=dx*dx+dy*dy+dz*dz
426
427 END DO
428
429 DO k=1,n4n
430 i=i4n(k)
431C
432 lap = one-lbp(i,it)-lcp(i,it)
433C
434 gapv(i,it) = min(gaps(i)+lap*gap_mm(i)+lbp(i,it)*gap_nm(i1,i)+lcp(i,it)*gap_nm(i2,i),gapmxl(i))
435 gap2 = gapv(i,it)**2
436C
437 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))
438
439 l = cand_e(i)
440 IF(etyp(i)==0)THEN
441 IF(bb(i,it) > zero)THEN
442 pent(i,it)=gapv(i,it)+bb(i,it)
443 ELSE
444 pent(i,it)=zero
445 END IF
446 ELSEIF(etyp(i) > nrtm)THEN ! Coating shell (same orientation than solid)
447 IF(bb(i,it) > zero)THEN
448 pent(i,it)=gapv(i,it)+bb(i,it)
449 ELSE
450 pent(i,it)=max(zero,gapv(i,it)-sqrt(dd(i,it)))
451 END IF
452 ELSE
453 IF(bb(i,it) > zero)THEN
454C
455C penetration approximee (cf zone limite interpolation des normales)
456 pent(i,it)=zero
457 ELSE
458 pent(i,it)=max(zero,gapv(i,it)-sqrt(dd(i,it)))
459 END IF
460 END IF
461C
462 ENDDO
463 END DO ! DO IT=2,4
464C--------------------------------------------------------
465 DO k=1,nindx
466 i=indx(k)
467C
468C Erase Subtria (cf intersections)
469 subtria(i)=0
470C
471 IF(ix3(i)/=ix4(i))THEN
472 dmin=min(dd(i,1),dd(i,2),dd(i,3),dd(i,4))
473 DO jj=1,4
474 IF(dd(i,jj) <= onep03*dmin)THEN
475 lbx = lb(i,jj)
476 lcx = lc(i,jj)
477 lax = one-lb(i,jj)-lc(i,jj)
478C
479C Privilegier secteur dans lequel on se trouve.
480 IF(lbx >= zero .AND. lcx >= zero)THEN
481 lx=zero
482 ELSE
483 lx=max(zero,dd(i,jj)-bb(i,jj)*bb(i,jj))
484 END IF
485
486 IF(lx < ld(i)) THEN
487 subtria(i)= jj
488 dist(i) = dd(i,jj)
489 ld(i)=lx
490 END IF
491
492 END IF
493 END DO
494 ELSE
495 IF(dd(i,1) <= dist(i))THEN
496 subtria(i)= 1
497 dist(i) = dd(i,1)
498 END IF
499 END IF
500C
501 it = subtria(i)
502 DO j=1,4
503 IF(j /= it) pent(i,j)=zero
504 END DO
505C
506 END DO
507C--------------------------------------------------------
508C Look for nodes projecting out of the free edges <=> Zero Penetration
509C--------------------------------------------------------
510 DO k=1,nindx
511 i=indx(k)
512C
513 it=subtria(i)
514C IF(PENT(I,IT)==ZERO) CYCLE ! In Starter, needs to check for FAR=2
515
516 i1=itria(1,it)
517 i2=itria(2,it)
518C
519 xh(i)=xi(i)+bb(i,it)*xn(i,it)
520 yh(i)=yi(i)+bb(i,it)*yn(i,it)
521 zh(i)=zi(i)+bb(i,it)*zn(i,it)
522C
523C Projection en dehors du triangle
524 IF(lb(i,it) < -em03 .OR. lc(i,it) < -em03 .OR.
525 . lb(i,it)+lc(i,it) > one+em03) far(i,it)=1
526C
527 IF(ix3(i) /= ix4(i))THEN
528C
529 ib1=ibound(i1,i)
530 ib2=ibound(i2,i)
531 IF(mvoisn(i,it)==0)THEN
532C
533 IF( (xh(i)-xx(i,i1))* nnx(i,it)+
534 . (yh(i)-yy(i,i1))* nny(i,it)+
535 . (zh(i)-zz(i,i1))* nnz(i,it) >= gaps(i)) far(i,it) =2
536C
537 ELSEIF((ib1 /= 0 .AND. ib2 == 0).OR.
538 . (ib2 /= 0 .AND. ib1 == 0))THEN
539C
540 ibx=max(ib1,ib2)
541 IF(ib1/=0)THEN
542 ix =i1
543 ELSEIF(ib2/=0)THEN
544 ix =i2
545 END IF
546C
547 IF(vtx_bisector(1,1,ibx)/=zero.OR.
548 . vtx_bisector(2,1,ibx)/=zero.OR.
549 . vtx_bisector(3,1,ibx)/=zero.OR.
550 . vtx_bisector(1,2,ibx)/=zero.OR.
551 . vtx_bisector(2,2,ibx)/=zero.OR.
552 . vtx_bisector(3,2,ibx)/=zero)THEN
553 p1 = (xh(i)-xx(i,ix))* vtx_bisector(1,1,ibx)+
554 . (yh(i)-yy(i,ix))* vtx_bisector(2,1,ibx)+
555 . (zh(i)-zz(i,ix))* vtx_bisector(3,1,ibx)
556 p2 = (xh(i)-xx(i,ix))* vtx_bisector(1,2,ibx)+
557 . (yh(i)-yy(i,ix))* vtx_bisector(2,2,ibx)+
558 . (zh(i)-zz(i,ix))* vtx_bisector(3,2,ibx)
559 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
560 ELSE
561 vx = x0n(i,ix) ! fake bisector of angle at vertex IX
562 vy = y0n(i,ix)
563 vz = z0n(i,ix)
564 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
565 pn = ((xh(i)-xx(i,ix))*vx+(yh(i)-yy(i,ix))*vy+(zh(i)-zz(i,ix))*vz)*nn
566 IF(pn >= gaps(i)) far(i,it) =2
567 END IF
568C
569 END IF
570C
571 IF(far(i,it)==1 .OR. bb(i,it) <= zero)THEN
572C
573 x12(i)=xx(i,i2)-xx(i,i1)
574 y12(i)=yy(i,i2)-yy(i,i1)
575 z12(i)=zz(i,i2)-zz(i,i1)
576C
577C normal to the bisecting plane (pointing toward the inside)
578 px = z12(i)*nny(i,it)-y12(i)*nnz(i,it)
579 py = x12(i)*nnz(i,it)-z12(i)*nnx(i,it)
580 pz = y12(i)*nnx(i,it)-x12(i)*nny(i,it)
581 pp = px*px+py*py+pz*pz
582C
583 ll = xij(i,i1)*xij(i,i1)+yij(i,i1)*yij(i,i1)+zij(i,i1)*zij(i,i1)
584C
585 side(i,it)=-(xij(i,i1)*px+yij(i,i1)*py+zij(i,i1)*pz)*sqrt(one/max(em30,ll*pp))
586C
587C Projette a l'exterieur et Hors secteur angulaire (cf bissectrices) <=> Glissement
588 IF(side(i,it) < -zep01) far(i,it) =2
589C
590 END IF
591C
592 ELSE
593C
594 ib1=ibound(1,i)
595 ib2=ibound(2,i)
596 ib3=ibound(3,i)
597C
598 d1=(xh(i)-xx(i,1))* nnx(i,1)+
599 . (yh(i)-yy(i,1))* nny(i,1)+
600 . (zh(i)-zz(i,1))* nnz(i,1)
601 d2=(xh(i)-xx(i,2))* nnx(i,2)+
602 . (yh(i)-yy(i,2))* nny(i,2)+
603 . (zh(i)-zz(i,2))* nnz(i,2)
604 d3=(xh(i)-xx(i,3))* nnx(i,4)+
605 . (yh(i)-yy(i,3))* nny(i,4)+
606 . (zh(i)-zz(i,3))* nnz(i,4)
607C
608 IF( (mvoisn(i,1) == 0 .AND. d1 >= gaps(i)).OR.
609 . (mvoisn(i,2) == 0 .AND. d2 >= gaps(i)).OR.
610 . (mvoisn(i,4) == 0 .AND. d3 >= gaps(i)) )THEN
611C
612 far(i,it)=2
613C
614 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
615 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
616 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
617C
618 ibx=max(ib1,ib2,ib3)
619 IF(ib1/=0)THEN
620 ix =1
621 iy =2
622 iz =3
623 ELSEIF(ib2/=0)THEN
624 ix =2
625 iy =3
626 iz =1
627 ELSEIF(ib3/=0)THEN
628 ix =3
629 iy =1
630 iz =2
631 END IF
632C
633 IF(vtx_bisector(1,1,ibx)/=zero.OR.
634 . vtx_bisector(2,1,ibx)/=zero.OR.
635 . vtx_bisector(3,1,ibx)/=zero.OR.
636 . vtx_bisector(1,2,ibx)/=zero.OR.
637 . vtx_bisector(2,2,ibx)/=zero.OR.
638 . vtx_bisector(3,2,ibx)/=zero)THEN
639 p1 = (xh(i)-xx(i,ix))* vtx_bisector(1,1,ibx)+
640 . (yh(i)-yy(i,ix))* vtx_bisector(2,1,ibx)+
641 . (zh(i)-zz(i,ix))* vtx_bisector(3,1,ibx)
642 p2 = (xh(i)-xx(i,ix))* vtx_bisector(1,2,ibx)+
643 . (yh(i)-yy(i,ix))* vtx_bisector(2,2,ibx)+
644 . (zh(i)-zz(i,ix))* vtx_bisector(3,2,ibx)
645 IF(p1 >= gaps(i) .AND. p2 >= gaps(i)) far(i,it) =2
646 ELSE
647 vx = two*xx(i,ix)-(xx(i,iy)+xx(i,iz)) ! fake bisector of angle 2,1,3
648 vy = two*yy(i,ix)-(yy(i,iy)+yy(i,iz))
649 vz = two*zz(i,ix)-(zz(i,iy)+zz(i,iz))
650 nn = one/max(em20,sqrt(vx*vx+vy*vy+vz*vz))
651 pn = ((xh(i)-xx(i,ix))*vx+(yh(i)-yy(i,ix))*vy+(zh(i)-zz(i,ix))*vz)*nn
652 IF(pn >= gaps(i)) far(i,it) =2
653 END IF
654C
655 END IF
656C
657 IF(far(i,it)==1 .OR. bb(i,it) <= zero)THEN
658
659 IF(mvoisn(i,1) /= 0)THEN
660C
661C INGAP = 1 => tjrs regarder si l'on est dans le cone (permet de determiner dans quel cone on se trouve)
662C
663C Glissement => utiliser le cone pour choisir entre les segments voisins
664C mais garder le contact si le nd a glisse de plus d'un elt (cf supersonic)
665C
666 x12(i)=xx(i,2)-xx(i,1)
667 y12(i)=yy(i,2)-yy(i,1)
668 z12(i)=zz(i,2)-zz(i,1)
669C
670C normal to the bisecting plane (pointing toward the inside)
671 px = z12(i)*nny(i,1)-y12(i)*nnz(i,1)
672 py = x12(i)*nnz(i,1)-z12(i)*nnx(i,1)
673 pz = y12(i)*nnx(i,1)-x12(i)*nny(i,1)
674 pp = px*px+py*py+pz*pz
675C
676 ll = xij(i,1)*xij(i,1)+yij(i,1)*yij(i,1)+zij(i,1)*zij(i,1)
677C
678 side(i,1)=-(xij(i,1)*px+yij(i,1)*py+zij(i,1)*pz)*sqrt(one/max(em30,ll*pp))
679 IF(side(i,1) < -zep01) far(i,it) =2
680C
681 END IF
682
683 IF(mvoisn(i,2) /= 0)THEN
684C
685C INGAP = 1 => tjrs regarder si l'on est dans le cone (permet de determiner dans quel cone on se trouve)
686C
687C Glissement => utiliser le cone pour choisir entre les segments voisins
688C mais garder le contact si le nd a glisse de plus d'un elt (cf supersonic)
689C
690 x12(i)=xx(i,3)-xx(i,2)
691 y12(i)=yy(i,3)-yy(i,2)
692 z12(i)=zz(i,3)-zz(i,2)
693C
694C normal to the bisecting plane (pointing toward the inside)
695 px = z12(i)*nny(i,2)-y12(i)*nnz(i,2)
696 py = x12(i)*nnz(i,2)-z12(i)*nnx(i,2)
697 pz = y12(i)*nnx(i,2)-x12(i)*nny(i,2)
698 pp = px*px+py*py+pz*pz
699C
700 ll = xij(i,2)*xij(i,2)+yij(i,2)*yij(i,2)+zij(i,2)*zij(i,2)
701C
702 side(i,2)=-(xij(i,2)*px+yij(i,2)*py+zij(i,2)*pz)*sqrt(one/max(em30,ll*pp))
703 IF(side(i,2) < -zep01) far(i,it) =2
704C
705 END IF
706
707 IF(mvoisn(i,4) /= 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,1)-xx(i,3)
715 y12(i)=yy(i,1)-yy(i,3)
716 z12(i)=zz(i,1)-zz(i,3)
717C
718C normal to the bisecting plane (pointing toward the inside)
719 px = z12(i)*nny(i,4)-y12(i)*nnz(i,4)
720 py = x12(i)*nnz(i,4)-z12(i)*nnx(i,4)
721 pz = y12(i)*nnx(i,4)-x12(i)*nny(i,4)
722 pp = px*px+py*py+pz*pz
723C
724 ll = xij(i,3)*xij(i,3)+yij(i,3)*yij(i,3)+zij(i,3)*zij(i,3)
725C
726 side(i,4)=-(xij(i,3)*px+yij(i,3)*py+zij(i,3)*pz)*sqrt(one/max(em30,ll*pp))
727 IF(side(i,4) < -zep01) far(i,it) =2
728C
729 END IF
730 END IF
731C
732 END IF
733 IF(far(i,it)==2) pent(i,it)=zero
734C
735C-------------------------------------------
736c if(itab(nsvg(i))==2810875.or.
737c . itab(ix1(i))==2810875.or.itab(ix2(i))==2810875.or.itab(ix3(i))==2810875.or.itab(ix4(i))==2810875)
738c . print *,'avant',itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),mvoisn(i,1:4),
739c . pent(i,it)
740C-------------------------------------------
741
742 ENDDO
743C--------------------------------------------------------
744C Exact value of the penetration close to the edges still needs to be computed
745C--------------------------------------------------------
746 nbord=0
747 DO k=1,nindx
748 i=indx(k)
749C
750 it=subtria(i)
751 i1=itria(1,it)
752 i2=itria(2,it)
753C
754 nn1(i) = zero
755 nn2(i) = zero
756 nn3(i) = zero
757
758 IF(pent(i,it)==zero) THEN
759 cycle
760 ENDIF
761C
762 IF(ix3(i) /= ix4(i))THEN
763C
764 ib1=ibound(i1,i)
765 ib2=ibound(i2,i)
766C
767 IF(mvoisn(i,it)==0)THEN
768C
769C Distance signee a l'arete
770C
771C upper skin --------------------
772C |
773C I |
774C -BB: | |
775C | NNE |
776C neutral fiber - - - C - - H < - - >
777C <-------------->|
778C Gap |
779C |
780C --------------------
781C
782 nn1(i)=nnx(i,it)
783 nn2(i)=nny(i,it)
784 nn3(i)=nnz(i,it)
785 nne(i)= (xh(i)-xx(i,i1))*nn1(i)+ (yh(i)-yy(i,i1))*nn2(i)+ (zh(i)-zz(i,i1))*nn3(i)
786C
787 nbord=nbord+1
788 kbord(nbord)=i
789C
790 ELSEIF((ib1/=0.AND.ib2==0).OR.
791 . (ib2/=0.AND.ib1==0))THEN
792C
793 ibx=max(ib1,ib2)
794 IF(ib1/=0)THEN
795 ix =i1
796 ELSEIF(ib2/=0)THEN
797 ix =i2
798 END IF
799C
800 p1 = (xh(i)-xx(i,ix))* vtx_bisector(1,1,ibx)+
801 . (yh(i)-yy(i,ix))* vtx_bisector(2,1,ibx)+
802 . (zh(i)-zz(i,ix))* vtx_bisector(3,1,ibx)
803 p2 = (xh(i)-xx(i,ix))* vtx_bisector(1,2,ibx)+
804 . (yh(i)-yy(i,ix))* vtx_bisector(2,2,ibx)+
805 . (zh(i)-zz(i,ix))* vtx_bisector(3,2,ibx)
806C
807 IF(p1 < gaps(i) .AND. p2 < gaps(i))THEN
808C
809 nn1(i)=vtx_bisector(1,1,ibx)+vtx_bisector(1,2,ibx)
810 nn2(i)=vtx_bisector(2,1,ibx)+vtx_bisector(2,2,ibx)
811 nn3(i)=vtx_bisector(3,1,ibx)+vtx_bisector(3,2,ibx)
812C
813 nn=sqrt(nn1(i)*nn1(i)+nn2(i)*nn2(i)+nn3(i)*nn3(i))
814 nn = one/max(em30,nn)
815 nn1(i)=nn1(i)*nn
816 nn2(i)=nn2(i)*nn
817 nn3(i)=nn3(i)*nn
818C
819 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
820C
821 nbord=nbord+1
822 kbord(nbord)=i
823C
824 ELSEIF(p1 < gaps(i))THEN
825
826 nn1(i)= vtx_bisector(1,1,ibx)
827 nn2(i)= vtx_bisector(2,1,ibx)
828 nn3(i)= vtx_bisector(3,1,ibx)
829 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
830C
831 nbord=nbord+1
832 kbord(nbord)=i
833
834 ELSEIF(p2 < gaps(i))THEN
835
836 nn1(i)= vtx_bisector(1,2,ibx)
837 nn2(i)= vtx_bisector(2,2,ibx)
838 nn3(i)= vtx_bisector(3,2,ibx)
839 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
840C
841 nbord=nbord+1
842 kbord(nbord)=i
843
844 ELSE
845
846C print *,' ** possible internal error wrt p1,p2 in i25pen3',ib1,ib2,
847C . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),p1,p2
848
849 END IF
850
851 END IF
852C
853 ELSE ! IF(IX3(I) /= IX4(I))THEN
854
855 ib1=ibound(1,i)
856 ib2=ibound(2,i)
857 ib3=ibound(3,i)
858C
859C Projection sur le plan et Distance signee au plan
860 xh(i)=xi(i)+bb(i,it)*xn(i,it)
861 yh(i)=yi(i)+bb(i,it)*yn(i,it)
862 zh(i)=zi(i)+bb(i,it)*zn(i,it)
863C
864 IF(mvoisn(i,1)==0.OR.
865 . mvoisn(i,2)==0.OR.
866 . mvoisn(i,4)==0)THEN
867
868 nne(i)=gaps(i)
869 IF(mvoisn(i,1)==0)THEN
870C
871 nn = (xh(i)-xx(i,i1))*nnx(i,i1)+(yh(i)-yy(i,i1))*nny(i,i1)+(zh(i)-zz(i,i1))*nnz(i,i1)
872
873 IF(nn < nne(i)) THEN
874 nne(i)=nn
875 nn1(i)=nnx(i,i1)
876 nn2(i)=nny(i,i1)
877 nn3(i)=nnz(i,i1)
878C
879 nbord=nbord+1
880 kbord(nbord)=i
881
882 END IF
883C
884 END IF
885C
886 IF(mvoisn(i,2)==0)THEN
887C
888C Distance signee a l'arete
889 nn = (xh(i)-xx(i,i2))*nnx(i,i2)+(yh(i)-yy(i,i2))*nny(i,i2)+(zh(i)-zz(i,i2))*nnz(i,i2)
890
891 IF(nn < nne(i)) THEN
892 nne(i)=nn
893 nn1(i)=nnx(i,i2)
894 nn2(i)=nny(i,i2)
895 nn3(i)=nnz(i,i2)
896C
897 nbord=nbord+1
898 kbord(nbord)=i
899
900 END IF
901C
902 END IF
903C
904 IF(mvoisn(i,4)==0)THEN
905C
906C Distance signee a l'arete
907 nn = (xh(i)-xx(i,5))*nnx(i,5)+(yh(i)-yy(i,5))*nny(i,5)+(zh(i)-zz(i,5))*nnz(i,5)
908
909 IF(nn < nne(i)) THEN
910 nne(i)=nn
911 nn1(i)=nnx(i,5)
912 nn2(i)=nny(i,5)
913 nn3(i)=nnz(i,5)
914C
915 nbord=nbord+1
916 kbord(nbord)=i
917
918 END IF
919C
920 END IF
921C
922 ELSEIF((ib1/=0 .AND. ib2==0 .AND. ib3==0).OR.
923 . (ib2/=0 .AND. ib3==0 .AND. ib1==0).OR.
924 . (ib3/=0 .AND. ib1==0 .AND. ib2==0))THEN
925C
926 ibx=max(ib1,ib2,ib3)
927 IF(ib1/=0)THEN
928 ix =1
929 ELSEIF(ib2/=0)THEN
930 ix =2
931 ELSEIF(ib3/=0)THEN
932 ix =3
933 END IF
934C
935 p1 = (xh(i)-xx(i,ix))* vtx_bisector(1,1,ibx)+
936 . (yh(i)-yy(i,ix))* vtx_bisector(2,1,ibx)+
937 . (zh(i)-zz(i,ix))* vtx_bisector(3,1,ibx)
938 p2 = (xh(i)-xx(i,ix))* vtx_bisector(1,2,ibx)+
939 . (yh(i)-yy(i,ix))* vtx_bisector(2,2,ibx)+
940 . (zh(i)-zz(i,ix))* vtx_bisector(3,2,ibx)
941
942 IF(p1 < gaps(i) .AND. p2 < gaps(i))THEN
943C
944 nn1(i)=vtx_bisector(1,1,ibx)+vtx_bisector(1,2,ibx)
945 nn2(i)=vtx_bisector(2,1,ibx)+vtx_bisector(2,2,ibx)
946 nn3(i)=vtx_bisector(3,1,ibx)+vtx_bisector(3,2,ibx)
947C
948 nn=sqrt(nn1(i)*nn1(i)+nn2(i)*nn2(i)+nn3(i)*nn3(i))
949 nn = one/max(em30,nn)
950 nn1(i)=nn1(i)*nn
951 nn2(i)=nn2(i)*nn
952 nn3(i)=nn3(i)*nn
953C
954 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
955C
956 nbord=nbord+1
957 kbord(nbord)=i
958C
959 ELSEIF(p1 < gaps(i))THEN
960
961 nn1(i)= vtx_bisector(1,1,ibx)
962 nn2(i)= vtx_bisector(2,1,ibx)
963 nn3(i)= vtx_bisector(3,1,ibx)
964 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
965C
966 nbord=nbord+1
967 kbord(nbord)=i
968
969 ELSEIF(p2 < gaps(i))THEN
970
971 nn1(i)= vtx_bisector(1,2,ibx)
972 nn2(i)= vtx_bisector(2,2,ibx)
973 nn3(i)= vtx_bisector(3,2,ibx)
974 nne(i)= (xh(i)-xx(i,ix))*nn1(i)+ (yh(i)-yy(i,ix))*nn2(i)+ (zh(i)-zz(i,ix))*nn3(i)
975C
976 nbord=nbord+1
977 kbord(nbord)=i
978
979 ELSE
980
981c print *,' ** possible internal error wrt p1,p2 in i25pen3',ib1,ib2,
982c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),p1,p2
983
984 END IF
985 END IF
986
987 END IF ! IF(IX3(I) /= IX4(I))THEN
988 END DO
989C-------------------------------------------
990 IF(isharp==1)THEN
991
992 DO k=1,nbord
993 i=kbord(k)
994
995 it=subtria(i)
996 i1=itria(1,it)
997 i2=itria(2,it)
998
999 gapm = gapv(i,it)-gaps(i)
1000 IF(nne(i) > zero .AND. bb(i,it)+gapm < zero)THEN
1001C
1002C ___________________________ xxx
1003C |xxxxxx
1004C GapS |xxxxxxxx <=> Upper Round Corner
1005C ___________________________|xxxxxxxx
1006C | |
1007C GapM | | <=> Straight shape of the edge
1008C ---------------------------M-------
1009C <------>
1010C GapS
1011C
1012
1013C
1014 lap = one-lbp(i,it)-lcp(i,it)
1015 xp=lap*xx(i,5)+lbp(i,it)*xx(i,i1) + lcp(i,it)*xx(i,i2)
1016 yp=lap*yy(i,5)+lbp(i,it)*yy(i,i1) + lcp(i,it)*yy(i,i2)
1017 zp=lap*zz(i,5)+lbp(i,it)*zz(i,i1) + lcp(i,it)*zz(i,i2)
1018 xc=xp+gapm*xn(i,it)
1019 yc=yp+gapm*yn(i,it)
1020 zc=zp+gapm*zn(i,it)
1021 nn1(i) =xi(i)-xc
1022 nn2(i) =yi(i)-yc
1023 nn3(i) =zi(i)-zc
1024 dc = sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i))
1025
1026 IF(dc > em04) THEN
1027C NN = ONE/DC
1028C N1(I) = N1(I)*NN
1029C N2(I) = N2(I)*NN
1030C N3(I) = N3(I)*NN
1031 pent(i,it)=max(zero,gaps(i)-dc)
1032 ELSE
1033C
1034 nne(i)=nne(i)-gaps(i)
1035 IF(-bb(i,it) < gapv(i,it)+nne(i))THEN ! Test 45 degrees !
1036C
1037C Normale horizonale
1038C N1(I) = NN1(I)
1039C N2(I) = NN2(I)
1040C N3(I) = NN3(I)
1041C
1042 pent(i,it)=max(zero,-nne(i))
1043C
1044 ELSE
1045C
1046C Normale verticale !
1047C N1(I) = XN(I,IT)
1048C N2(I) = YN(I,IT)
1049C N3(I) = ZN(I,IT)
1050C
1051 pent(i,it)=max(zero,gapv(i,it)+bb(i,it))
1052 END IF
1053 END IF
1054
1055 ELSE !IF(NNE(I) > ZERO .AND. BB(I)+GAPM < ZERO)THEN
1056C
1057 nne(i)=nne(i)-gaps(i)
1058 IF(nne(i) >= zero)THEN
1059C
1060C Must remain a roundoff error <=> NNE ~ 0
1061c IF(NSVG(I) > 0)THEN
1062c print *,' ** possible internal error in i25dst3-3, nne should be negative',
1063c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I),pent(i,it)
1064c ELSE
1065c print *,' ** possible internal error in i25dst3-3, nne should be negative',
1066c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I),pent(i,it)
1067c END IF
1068
1069 pent(i,it)=zero
1070 cycle
1071
1072 END IF
1073
1074 IF(gapv(i,it)+nne(i) > zero)THEN
1075C
1076 IF(-bb(i,it) < gapv(i,it)+nne(i))THEN
1077C
1078C Normale horizonale
1079C N1(I) = NN1(I)
1080C N2(I) = NN2(I)
1081C N3(I) = NN3(I)
1082C
1083 pent(i,it)=-nne(i)
1084C
1085 ELSE
1086C
1087C Normale verticale !
1088C N1(I) = XN(I,IT)
1089C N2(I) = YN(I,IT)
1090C N3(I) = ZN(I,IT)
1091C
1092 pent(i,it)=max(zero,gapv(i,it)+bb(i,it))
1093 END IF
1094
1095 END IF
1096
1097 END IF
1098
1099 END DO
1100
1101 ELSEIF(isharp==2)THEN
1102
1103 DO k=1,nbord
1104 i=kbord(k)
1105
1106 it=subtria(i)
1107 i1=itria(1,it)
1108 i2=itria(2,it)
1109
1110 nne(i)=nne(i)-gaps(i)
1111 IF(nne(i) >= zero)THEN
1112C
1113C Must remain a roundoff error <=> NNE ~ 0
1114c IF(NSVG(I) > 0)THEN
1115c print *,' ** possible internal error in i25dst3-3, nne should be negative',
1116c . itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I)
1117c ELSE
1118c print *,' ** possible internal error in i25dst3-3, nne should be negative',
1119c . itafi(nin)%p(-nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),'nne =',nne(I)
1120c END IF
1121
1122 cycle
1123
1124 END IF
1125C
1126 IF(gapv(i,it)+nne(i) > zero)THEN
1127C
1128C Normale radiale shiftee
1129 xc =xh(i)-(gapv(i,it)+nne(i))*nn1(i)
1130 yc =yh(i)-(gapv(i,it)+nne(i))*nn2(i)
1131 zc =zh(i)-(gapv(i,it)+nne(i))*nn3(i)
1132 nn1(i) =xi(i)-xc
1133 nn2(i) =yi(i)-yc
1134 nn3(i) =zi(i)-zc
1135 dc = sqrt(nn1(i)*nn1(i)+ nn2(i)*nn2(i)+ nn3(i)*nn3(i))
1136 IF(dc > em04) THEN
1137C NN = ONE/DC
1138C N1(I) = NN1(I)*NN
1139C N2(I) = NN2(I)*NN
1140C N3(I) = NN3(I)*NN
1141 pent(i,it)=max(zero,gapv(i,it)-dc)
1142 ELSE
1143C
1144C Keep what was done in the general case : Normale ~ verticale.
1145 END IF
1146
1147 END IF
1148C-------------------------------------------
1149c if(itab(nsvg(i))==2810875.or.
1150c . itab(ix1(i))==2810875.or.itab(ix2(i))==2810875.or.itab(ix3(i))==2810875.or.itab(ix4(i))==2810875)
1151c . print *,'apres',itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),mvoisn(i,1:4),
1152c . it,kbord,gapv(i,it),nne(i),bb(i,it),pent(i,it)
1153C-------------------------------------------
1154 END DO
1155
1156 END IF ! IF(ISHARP==1)THEN
1157C--------------------------------------------------------
1158 IF(ivis2/=-1) THEN
1159
1160 DO k=1,nindx
1161 i=indx(k)
1162
1163 it = subtria(i)
1164
1165 IF(etyp(i) > nrtm)THEN
1166 IF(far(i,it)==2) cycle ! retain closest segment whose cone includes the secondary node
1167 ! => same impact (no sliding) at time=0 in RD Engine
1168 ! Case of coating shell, penetration may be larger than the gap
1169 ELSEIF(etyp(i)/=0)THEN
1170 IF(far(i,it)==2) cycle ! retain closest segment whose cone includes the SECONDARY node
1171 ! => same impact (no sliding) at time=0 in RD Engine
1172 IF(bb(i,it) > zero) cycle
1173 END IF
1174
1175 n=cand_n(i)
1176 l=cand_e(i)
1177
1178 mglob = mseglo(l)
1179
1180 pene=pent(i,it) ! PENE /= ZERO .OR. DIST(I) < PENMIN**2
1181
1182C------- exclude wrong normal dir PEN_OLD(1-3): normal of nn
1183 IF (msegtyp(l)==0.OR.msegtyp(l)>nrtm) THEN
1184 IF(inacti/=0.AND.ilev/=3)THEN
1185 norm = nn1(i)*pene_old(1,n)+nn2(i)*pene_old(2,n)
1186 + +nn3(i)*pene_old(3,n)
1187 IF(norm > zero) pene = zero
1188 END IF
1189 ENDIF
1190C if(itab(nsvg(i))==299344.or.itab(ix1(i))==299344.or.itab(ix2(i))==299344
1191C . .or.itab(ix3(i))==299344.or.itab(ix4(i))==299344)
1192C . print *,itab(nsvg(i)),itab(ix1(i)),itab(ix2(i)),itab(ix3(i)),itab(ix4(i)),
1193C . it,etyp(i),nrtm,pent(i,it),far(i,it),dist(i)
1194C
1195C Looking for penetration vs closest segment (cf initial penetrations vs solids)
1196 IF(dist(i) < time_s(n) .OR.
1197 * (dist(i) == time_s(n) .AND. irtlm(1,n) < mglob))THEN
1198 IF(far(i,it)==0 .AND. abs(dd(i,it)-gapv(i,it)*gapv(i,it)) < penmin**2)THEN
1199 irtlm(1,n)=mglob
1200 irtlm(2,n)=it
1201 irtlm(3,n)=l
1202 time_s(n) =dist(i)
1203 pene_old(5,n)=pene
1204 ELSE
1205 irtlm(1,n)=mglob
1206 irtlm(2,n)=it
1207 irtlm(3,n)=l
1208 time_s(n) =dist(i)
1209 pene_old(5,n)=pene
1210 END IF
1211 END IF
1212 END DO
1213 ELSE ! adhesion case, pene(5,n) datum is at 0.5*gapv, if pene < 0.5*gapv then pene_old(5,n) = 0
1214
1215 DO k=1,nindx
1216 i=indx(k)
1217
1218 it = subtria(i)
1219 IF(far(i,it)==2)THEN ! retain closest segment whose cone includes the SECONDARY node
1220 ! => same impact (no sliding) at time=0 in RD Engine
1221 it = 0
1222 subtria(i) = it
1223 END IF
1224 IF(it == 0)cycle
1225
1226 n=cand_n(i)
1227 l=cand_e(i)
1228
1229 mglob = mseglo(l)
1230
1231 IF(.NOT.(etyp(i)==0.OR.msegtyp(l) > nrtm))THEN
1232 IF(pent(i,it)==zero) cycle
1233 END IF
1234
1235
1236 pene=max(zero,pent(i,it)-half*gapv(i,it)) ! PENE /= ZERO .OR. DIST(I) < PENMIN**2
1237
1238C------- exclude wrong normal dir PEN_OLD(1-3): normal of nn
1239 IF (msegtyp(l)==0.OR.msegtyp(l)>nrtm) THEN
1240 IF(inacti/=0.AND.ilev/=3)THEN
1241 norm = nn1(i)*pene_old(1,n)+nn2(i)*pene_old(2,n)
1242 + +nn3(i)*pene_old(3,n)
1243 IF(norm > zero) pene = zero
1244 END IF
1245 ENDIF
1246C
1247C Looking for penetration vs closest segment (cf initial penetrations vs solids)
1248 IF(dist(i) < time_s(n) .OR.
1249 * (dist(i) == time_s(n) .AND. irtlm(1,n) < mglob))THEN
1250 IF(far(i,it)==0 .AND. abs(dd(i,it)-gapv(i,it)*gapv(i,it)) < penmin**2)THEN
1251 irtlm(1,n)=mglob
1252 irtlm(2,n)=it
1253 irtlm(3,n)=l
1254 time_s(n) =pene
1255 pene_old(5,n)=pene
1256 ELSE
1257 irtlm(1,n)=mglob
1258 irtlm(2,n)=it
1259 irtlm(3,n)=l
1260 time_s(n) =pene
1261 pene_old(5,n)=pene
1262 END IF
1263 END IF
1264 END DO
1265 ENDIF ! ivis2 if
1266C
1267 RETURN
1268 END
the
end diagonal values have been computed in the(sparse) matrix id.SOL
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
i25pen3
subroutine i25pen3(jlt, cand_n, cand_e, penmin, penmax, x1, x2, x3, x4, x0, y1, y2, y3, y4, y0, z1, z2, z3, z4, z0, xi, yi, zi, nsn, ix1, ix2, ix3, ix4, nsvg, nrtm, mseglo, gaps, irect, irtlm, time_s, pene_old, itab, msegtyp, isharp, nnx, nny, nnz, gap_nm, mvoisn, gapmxl, ivis2, ibound, vtx_bisector, ilev, inacti)
Definition i25pen3.F:43