OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i20dst3.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!|| i20dst3 ../starter/source/interfaces/inter3d1/i20dst3.F
25!||--- called by ------------------------------------------------------
26!|| i20ini3 ../starter/source/interfaces/inter3d1/i20ini3.F
27!||--- calls -----------------------------------------------------
28!|| bitget ../starter/source/interfaces/inter3d1/bitget.F
29!|| bitset ../starter/source/interfaces/inter3d1/bitget.F
30!|| bitunset ../starter/source/interfaces/inter3d1/bitget.F
31!|| i20cgap0 ../starter/source/interfaces/inter3d1/i20dst3.F
32!|| i20cgap1 ../starter/source/interfaces/inter3d1/i20dst3.F
33!||====================================================================
34 SUBROUTINE i20dst3(IGAP,GAP_SH,CAND_E,CAND_N,GAPV,
35 2 GAP ,GAP_S ,GAP_M ,GAPMAX,GAPMIN,
36 3 IRECT ,NLN ,NLG ,SOLIDN_NORMAL,NSV,
37 4 NBINFLG,TAG ,IX3 ,IX4 ,X1 ,
38 5 X2, X3, X4 ,Y1 ,Y2 ,
39 6 Y3, Y4, Z1 ,Z2 ,Z3 ,
40 7 Z4, XI, YI ,ZI ,X0 ,
41 8 Y0, Z0, NX1,NY1,NZ1,
42 9 NX2,NY2, NZ2,NX3,NY3,
43 1 NZ3,NX4, NY4,NZ4,P1 ,
44 2 P2 ,P3 ,P4 ,LB1,LB2,
45 3 LB3,LB4,LC1 ,LC2,LC3,
46 4 LC4)
47C-----------------------------------------------
48C I m p l i c i t T y p e s
49C-----------------------------------------------
50#include "implicit_f.inc"
51C-----------------------------------------------
52C G l o b a l P a r a m e t e r s
53C-----------------------------------------------
54#include "mvsiz_p.inc"
55C-----------------------------------------------
56C D u m m y A r g u m e n t s
57C-----------------------------------------------
58 INTEGER IGAP,CAND_N(*),CAND_E(*),IRECT(4,*),NLN,NLG(NLN),NSV(*),
59 . NBINFLG(*),TAG(*)
60 my_real
61 . GAP,GAP_SH(*),GAPV(MVSIZ),GAP_S(*),GAP_M(*),GAPMAX,GAPMIN,
62 . SOLIDN_NORMAL(3,*)
63 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IX3,IX4
64 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: X1,X2,X3,X4
65 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: Y1,Y2,Y3,Y4
66 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: Z1,Z2,Z3,Z4
67 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: XI,YI,ZI
68 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: X0,Y0,Z0
69 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: NX1,NY1,NZ1
70 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: nx2,ny2,nz2
71 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: nx3,ny3,nz3
72 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: nx4,ny4,nz4
73 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: p1,p2,p3,p4
74 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: lb1,lb2,lb3,lb4
75 my_real, DIMENSION(MVSIZ), INTENT(INOUT) :: lc1,lc2,lc3,lc4
76C-----------------------------------------------
77C C o m m o n B l o c k s
78C-----------------------------------------------
79#include "vect07_c.inc"
80C-----------------------------------------------
81C L o c a l V a r i a b l e s
82C-----------------------------------------------
83 INTEGER I,I1,I2,I3,I4,IG,IM,IS,IB1,IB2,IB3,IB4
84C REAL
85 my_real
86 . al1(mvsiz), al2(mvsiz), al3(mvsiz), al4(mvsiz),
87 . x01(mvsiz), x02(mvsiz), x03(mvsiz), x04(mvsiz),
88 . y01(mvsiz), y02(mvsiz), y03(mvsiz), y04(mvsiz),
89 . z01(mvsiz), z02(mvsiz), z03(mvsiz), z04(mvsiz),
90 . xi1(mvsiz), xi2(mvsiz), xi3(mvsiz), xi4(mvsiz),
91 . yi1(mvsiz), yi2(mvsiz), yi3(mvsiz), yi4(mvsiz),
92 . zi1(mvsiz), zi2(mvsiz), zi3(mvsiz), zi4(mvsiz),
93 . hlb1(mvsiz), hlc1(mvsiz), hlb2(mvsiz),hlc2(mvsiz),
94 . hlb3(mvsiz), hlc3(mvsiz), hlb4(mvsiz),hlc4(mvsiz)
95 my_real
96 . s2,a1,a2,a3,a4,d1,d2,d3,d4,
97 . x12,x23,x34,x41,xi0,sx1,sx2,sx3,sx4,sx0,sxi,
98 . y12,y23,y34,y41,yi0,sy1,sy2,sy3,sy4,sy0,syi,
99 . z12,z23,z34,z41,zi0,sz1,sz2,sz3,sz4,sz0,szi,
100 . x10,y10,z10,x20,y20,z20,x30,y30,z30,x40,y40,z40,
101 . la, hla, aaa, unssqr3
102C-----------------------------------------------
103 INTEGER BITUNSET,BITGET,BITSET
104 EXTERNAL BITUNSET,BITGET,BITSET
105C
106C--------------------------------------------------------
107C SHIFT DU GAP POUR SOLIDES AVEC GAP NUL
108C--------------------------------------------------------
109c
110c --------------------------
111c ^
112c |
113c |
114c |
115c |
116c | gap_search(search algorithm)
117c |
118c |
119c |
120c 0 Secnd |
121c | | ^
122c | gap_s | |
123c | | |
124c v | |
125c -------------------- | | gapv(penetration
126c ^ | | computation)
127c | gap_m | |
128c v v | gapv = min(gapmax,
129c Main1 0------------------0 M2 | gap_s+gap_m)
130c | ^ | + gap_sh
131c | gap_sh | |
132c | | |
133c v | v
134cM1-shift o------------------o |
135c ^ | gap_search(search algorithm)
136c | |
137c | | gap_search_engine = gapv + gap_sh
138c | gapv(penetration | = min(gapmax,
139c | computation) | gap_s+gap_m)
140c | | + 2*gap_sh
141c | gapv = min(gapmax, |
142c | gap_s+gap_m) | gap_search_starter =
143c | + gap_sh | gap_s + gap_m + 2*gap_sh
144c | | (overestimate)
145c ----------------------------
146c
147c
148C=======================================================================
149 unssqr3 = one/sqr3
150 IF(igap /= 0)THEN
151 DO i=lft,llt
152 im = cand_e(i)
153 is = cand_n(i)
154c IG = NLG(NSV(IS))
155 ig = nsv(is)
156 gapv(i) = gap_s(is) + gap_m(im)
157 gapv(i) = min(gapv(i),gapmax)
158 gapv(i) = gapv(i) + gap_sh(im)*(one-em5)
159 gapv(i) = max(gapmin,gapv(i))
160c I1 = NLG(IRECT(1,IM))
161c I2 = NLG(IRECT(2,IM))
162c I3 = NLG(IRECT(3,IM))
163c I4 = NLG(IRECT(4,IM))
164 i1 = irect(1,im)
165 i2 = irect(2,im)
166 i3 = irect(3,im)
167 i4 = irect(4,im)
168
169 ib1 = bitget(nbinflg(tag(i1)),7)
170 ib2 = bitget(nbinflg(tag(i2)),7)
171 ib3 = bitget(nbinflg(tag(i3)),7)
172 ib4 = bitget(nbinflg(tag(i4)),7)
173
174 IF(ib1+ib2+ib3+ib4 == 0)THEN
175 IF(gap_sh(im)/= zero)THEN
176C
177 sx0=(y3(i)-y1(i))*(z4(i)-z2(i))-(z3(i)-z1(i))*(y4(i)-y2(i))
178 sy0=(z3(i)-z1(i))*(x4(i)-x2(i))-(x3(i)-x1(i))*(z4(i)-z2(i))
179 sz0=(x3(i)-x1(i))*(y4(i)-y2(i))-(y3(i)-y1(i))*(x4(i)-x2(i))
180 aaa = one / sqrt(max(em20,sx0*sx0+sy0*sy0+sz0*sz0))
181 sx0=sx0*aaa
182 sy0=sy0*aaa
183 sz0=sz0*aaa
184 sxi = solidn_normal(1,ig)
185 syi = solidn_normal(2,ig)
186 szi = solidn_normal(3,ig)
187 aaa = sxi*sx0+syi*sy0+szi*sz0
188 IF(aaa > zero)THEN
189 sxi = zero
190 syi = zero
191 szi = zero
192 ENDIF
193
194 sx1 = solidn_normal(1,i1)-sxi
195 sy1 = solidn_normal(2,i1)-syi
196 sz1 = solidn_normal(3,i1)-szi
197 aaa = one / sqrt(max(em20,sx1*sx1+sy1*sy1+sz1*sz1))
198 sx1 = sx1*aaa
199 sy1 = sy1*aaa
200 sz1 = sz1*aaa
201 aaa = sx0*sx1 + sy0*sy1 + sz0*sz1
202 IF(aaa > unssqr3)THEN
203 aaa = gap_sh(im)/aaa
204 ELSE
205 aaa = unssqr3 - aaa
206 sx1 = sx1 + aaa*sx0
207 sy1 = sy1 + aaa*sy0
208 sz1 = sz1 + aaa*sz0
209 aaa = gap_sh(im)*sqr3
210 ENDIF
211 sx1 = sx1*aaa
212 sy1 = sy1*aaa
213 sz1 = sz1*aaa
214 sx2 = solidn_normal(1,i2)-sxi
215 sy2 = solidn_normal(2,i2)-syi
216 sz2 = solidn_normal(3,i2)-szi
217 aaa = one / sqrt(max(em20,sx2*sx2+sy2*sy2+sz2*sz2))
218 sx2 = sx2*aaa
219 sy2 = sy2*aaa
220 sz2 = sz2*aaa
221 aaa = sx0*sx2 + sy0*sy2 + sz0*sz2
222 IF(aaa > unssqr3)THEN
223 aaa = gap_sh(im)/aaa
224 ELSE
225 aaa = unssqr3 - aaa
226 sx2 = sx2 + aaa*sx0
227 sy2 = sy2 + aaa*sy0
228 sz2 = sz2 + aaa*sz0
229 aaa = gap_sh(im)*sqr3
230 ENDIF
231 sx2 = sx2*aaa
232 sy2 = sy2*aaa
233 sz2 = sz2*aaa
234 sx3 = solidn_normal(1,i3)-sxi
235 sy3 = solidn_normal(2,i3)-syi
236 sz3 = solidn_normal(3,i3)-szi
237 aaa = one / sqrt(max(em20,sx3*sx3+sy3*sy3+sz3*sz3))
238 sx3 = sx3*aaa
239 sy3 = sy3*aaa
240 sz3 = sz3*aaa
241 aaa = sx0*sx3 + sy0*sy3 + sz0*sz3
242 IF(aaa > unssqr3)THEN
243 aaa = gap_sh(im)/aaa
244 ELSE
245 aaa = unssqr3 - aaa
246 sx3 = sx3 + aaa*sx0
247 sy3 = sy3 + aaa*sy0
248 sz3 = sz3 + aaa*sz0
249 aaa = gap_sh(im)*sqr3
250 ENDIF
251 sx3 = sx3*aaa
252 sy3 = sy3*aaa
253 sz3 = sz3*aaa
254 sx4 = solidn_normal(1,i4)-sxi
255 sy4 = solidn_normal(2,i4)-syi
256 sz4 = solidn_normal(3,i4)-szi
257 aaa = one / sqrt(max(em20,sx4*sx4+sy4*sy4+sz4*sz4))
258 sx4 = sx4*aaa
259 sy4 = sy4*aaa
260 sz4 = sz4*aaa
261 aaa = sx0*sx4 + sy0*sy4 + sz0*sz4
262 IF(aaa > unssqr3)THEN
263 aaa = gap_sh(im)/aaa
264 ELSE
265 aaa = unssqr3 - aaa
266 sx4 = sx4 + aaa*sx0
267 sy4 = sy4 + aaa*sy0
268 sz4 = sz4 + aaa*sz0
269 aaa = gap_sh(im)*sqr3
270 ENDIF
271 sx4 = sx4*aaa
272 sy4 = sy4*aaa
273 sz4 = sz4*aaa
274 x1(i) = x1(i) - sx1
275 y1(i) = y1(i) - sy1
276 z1(i) = z1(i) - sz1
277 x2(i) = x2(i) - sx2
278 y2(i) = y2(i) - sy2
279 z2(i) = z2(i) - sz2
280 x3(i) = x3(i) - sx3
281 y3(i) = y3(i) - sy3
282 z3(i) = z3(i) - sz3
283 x4(i) = x4(i) - sx4
284 y4(i) = y4(i) - sy4
285 z4(i) = z4(i) - sz4
286 ENDIF
287 ELSEIF(i3 == i4)THEN
288c bord de coque3N
289 x10 = x1(i)
290 y10 = y1(i)
291 z10 = z1(i)
292 x20 = x2(i)
293 y20 = y2(i)
294 z20 = z2(i)
295 x30 = x3(i)
296 y30 = y3(i)
297 z30 = z3(i)
298 IF(ib1 == 1 .and. ib2 == 1)THEN
299 CALL i20cgap0(x10,y10,z10,x20,y20,z20,
300 . x30,y30,z30,x30,y30,z30,
301 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
302 . x3(i),y3(i),z3(i),x3(i),y3(i),z3(i),
303 . gap_m(im))
304 ENDIF
305 IF(ib2 == 1 .and. ib3 == 1)THEN
306 CALL i20cgap0(x20,y20,z20,x30,y30,z30,
307 . x10,y10,z10,x10,y10,z10,
308 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
309 . x1(i),y1(i),z1(i),x1(i),y1(i),z1(i),
310 . gap_m(im))
311 ENDIF
312 IF(ib3 == 1 .and. ib1 == 1)THEN
313 CALL i20cgap0(x30,y30,z30,x10,y10,z10,
314 . x20,y20,z20,x20,y20,z20,
315 . x3(i),y3(i),z3(i),x1(i),y1(i),z1(i),
316 . x2(i),y2(i),z2(i),x2(i),y2(i),z2(i),
317 . gap_m(im))
318 ENDIF
319 IF(ib1 == 1 .and. ib2+ib3 == 0)THEN
320 CALL i20cgap1(x10,y10,z10,x20,y20,z20,
321 . x30,y30,z30,
322 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
323 . x3(i),y3(i),z3(i),gap_m(im))
324 ENDIF
325 IF(ib2 == 1 .and. ib3+ib1 == 0)THEN
326 CALL i20cgap1(x20,y20,z20,x30,y30,z30,
327 . x10,y10,z10,
328 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
329 . x1(i),y1(i),z1(i),gap_m(im))
330 ENDIF
331 IF(ib3 == 1 .and. ib1+ib2 == 0)THEN
332 CALL i20cgap1(x30,y30,z30,x10,y10,z10,
333 . x20,y20,z20,
334 . x3(i),y3(i),z3(i),x1(i),y1(i),z1(i),
335 . x2(i),y2(i),z2(i),gap_m(im))
336 ENDIF
337 x4(i)=x3(i)
338 y4(i)=y3(i)
339 z4(i)=z3(i)
340 ELSE
341c bord de coque
342 x10 = x1(i)
343 y10 = y1(i)
344 z10 = z1(i)
345 x20 = x2(i)
346 y20 = y2(i)
347 z20 = z2(i)
348 x30 = x3(i)
349 y30 = y3(i)
350 z30 = z3(i)
351 x40 = x4(i)
352 y40 = y4(i)
353 z40 = z4(i)
354 IF(ib1 == 1 .and. ib2 == 1)THEN
355 CALL i20cgap0(x10,y10,z10,x20,y20,z20,
356 . x30,y30,z30,x40,y40,z40,
357 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
358 . x3(i),y3(i),z3(i),x4(i),y4(i),z4(i),
359 . gap_m(im))
360 ENDIF
361 IF(ib2 == 1 .and. ib3 == 1)THEN
362 CALL i20cgap0(x20,y20,z20,x30,y30,z30,
363 . x40,y40,z40,x10,y10,z10,
364 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
365 . x4(i),y4(i),z4(i),x1(i),y1(i),z1(i),
366 . gap_m(im))
367 ENDIF
368 IF(ib3 == 1 .and. ib4 == 1)THEN
369 CALL i20cgap0(x30,y30,z30,x40,y40,z40,
370 . x10,y10,z10,x20,y20,z20,
371 . x3(i),y3(i),z3(i),x4(i),y4(i),z4(i),
372 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
373 . gap_m(im))
374 ENDIF
375 IF(ib4 == 1 .and. ib1 == 1)THEN
376 CALL i20cgap0(x40,y40,z40,x10,y10,z10,
377 . x20,y20,z20,x30,y30,z30,
378 . x4(i),y4(i),z4(i),x1(i),y1(i),z1(i),
379 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
380 . gap_m(im))
381 ENDIF
382 IF(ib1 == 1 .and. ib2+ib4 == 0)THEN
383 CALL i20cgap1(x10,y10,z10,x20,y20,z20,
384 . x40,y40,z40,
385 . x1(i),y1(i),z1(i),x2(i),y2(i),z2(i),
386 . x4(i),y4(i),z4(i),gap_m(im))
387 ENDIF
388 IF(ib2 == 1 .and. ib3+ib1 == 0)THEN
389 CALL i20cgap1(x20,y20,z20,x30,y30,z30,
390 . x10,y10,z10,
391 . x2(i),y2(i),z2(i),x3(i),y3(i),z3(i),
392 . x1(i),y1(i),z1(i),gap_m(im))
393 ENDIF
394 IF(ib3 == 1 .and. ib4+ib2 == 0)THEN
395 CALL i20cgap1(x30,y30,z30,x40,y40,z40,
396 . x10,y10,z10,
397 . x3(i),y3(i),z3(i),x4(i),y4(i),z4(i),
398 . x2(i),y2(i),z2(i),gap_m(im))
399 ENDIF
400 IF(ib4 == 1 .and. ib1+ib3 == 0)THEN
401 CALL i20cgap1(x40,y40,z40,x10,y10,z10,
402 . x30,y30,z30,
403 . x4(i),y4(i),z4(i),x1(i),y1(i),z1(i),
404 . x3(i),y3(i),z3(i),gap_m(im))
405 ENDIF
406 ENDIF
407 ENDDO
408 ENDIF
409
410 DO i=lft,llt
411 x0(i) = fourth*(x1(i)+x2(i)+x3(i)+x4(i))
412 y0(i) = fourth*(y1(i)+y2(i)+y3(i)+y4(i))
413 z0(i) = fourth*(z1(i)+z2(i)+z3(i)+z4(i))
414 ENDDO
415C
416 DO i=lft,llt
417 IF(ix3(i)==ix4(i))THEN
418 x0(i) = x3(i)
419 y0(i) = y3(i)
420 z0(i) = z3(i)
421 ENDIF
422 ENDDO
423C
424 DO i=lft,llt
425C
426 x01(i) = x1(i) - x0(i)
427 y01(i) = y1(i) - y0(i)
428 z01(i) = z1(i) - z0(i)
429C
430 x02(i) = x2(i) - x0(i)
431 y02(i) = y2(i) - y0(i)
432 z02(i) = z2(i) - z0(i)
433C
434 x03(i) = x3(i) - x0(i)
435 y03(i) = y3(i) - y0(i)
436 z03(i) = z3(i) - z0(i)
437C
438 x04(i) = x4(i) - x0(i)
439 y04(i) = y4(i) - y0(i)
440 z04(i) = z4(i) - z0(i)
441C
442 xi0 = x0(i) - xi(i)
443 yi0 = y0(i) - yi(i)
444 zi0 = z0(i) - zi(i)
445C
446 xi1(i) = x1(i) - xi(i)
447 yi1(i) = y1(i) - yi(i)
448 zi1(i) = z1(i) - zi(i)
449C
450 xi2(i) = x2(i) - xi(i)
451 yi2(i) = y2(i) - yi(i)
452 zi2(i) = z2(i) - zi(i)
453C
454 xi3(i) = x3(i) - xi(i)
455 yi3(i) = y3(i) - yi(i)
456 zi3(i) = z3(i) - zi(i)
457C
458 xi4(i) = x4(i) - xi(i)
459 yi4(i) = y4(i) - yi(i)
460 zi4(i) = z4(i) - zi(i)
461C
462 sx1 = yi0*zi1(i) - zi0*yi1(i)
463 sy1 = zi0*xi1(i) - xi0*zi1(i)
464 sz1 = xi0*yi1(i) - yi0*xi1(i)
465C
466 sx2 = yi0*zi2(i) - zi0*yi2(i)
467 sy2 = zi0*xi2(i) - xi0*zi2(i)
468 sz2 = xi0*yi2(i) - yi0*xi2(i)
469C
470 sx0 = y01(i)*z02(i) - z01(i)*y02(i)
471 sy0 = z01(i)*x02(i) - x01(i)*z02(i)
472 sz0 = x01(i)*y02(i) - y01(i)*x02(i)
473 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
474C
475 lb1(i) = -(sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
476 lc1(i) = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
477C
478 sx3 = yi0*zi3(i) - zi0*yi3(i)
479 sy3 = zi0*xi3(i) - xi0*zi3(i)
480 sz3 = xi0*yi3(i) - yi0*xi3(i)
481C
482 sx0 = y02(i)*z03(i) - z02(i)*y03(i)
483 sy0 = z02(i)*x03(i) - x02(i)*z03(i)
484 sz0 = x02(i)*y03(i) - y02(i)*x03(i)
485 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
486C
487 lb2(i) = -(sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
488 lc2(i) = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
489C
490 sx4 = yi0*zi4(i) - zi0*yi4(i)
491 sy4 = zi0*xi4(i) - xi0*zi4(i)
492 sz4 = xi0*yi4(i) - yi0*xi4(i)
493C
494 sx0 = y03(i)*z04(i) - z03(i)*y04(i)
495 sy0 = z03(i)*x04(i) - x03(i)*z04(i)
496 sz0 = x03(i)*y04(i) - y03(i)*x04(i)
497 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
498C
499 lb3(i) = -(sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
500 lc3(i) = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
501C
502 sx0 = y04(i)*z01(i) - z04(i)*y01(i)
503 sy0 = z04(i)*x01(i) - x04(i)*z01(i)
504 sz0 = x04(i)*y01(i) - y04(i)*x01(i)
505 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
506C
507 lb4(i) = -(sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
508 lc4(i) = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
509C
510 aaa = one/max(em30,x01(i)*x01(i)+y01(i)*y01(i)+z01(i)*z01(i))
511 hlc1(i)= lc1(i)*abs(lc1(i))*aaa
512 hlb4(i)= lb4(i)*abs(lb4(i))*aaa
513 al1(i) = -(xi0*x01(i)+yi0*y01(i)+zi0*z01(i))*aaa
514 al1(i) = max(zero,min(one,al1(i)))
515 aaa = one/max(em30,x02(i)*x02(i)+y02(i)*y02(i)+z02(i)*z02(i))
516 hlc2(i)= lc2(i)*abs(lc2(i))*aaa
517 hlb1(i)= lb1(i)*abs(lb1(i))*aaa
518 al2(i) = -(xi0*x02(i)+yi0*y02(i)+zi0*z02(i))*aaa
519 al2(i) = max(zero,min(one,al2(i)))
520 aaa = one/max(em30,x03(i)*x03(i)+y03(i)*y03(i)+z03(i)*z03(i))
521 hlc3(i)= lc3(i)*abs(lc3(i))*aaa
522 hlb2(i)= lb2(i)*abs(lb2(i))*aaa
523 al3(i) = -(xi0*x03(i)+yi0*y03(i)+zi0*z03(i))*aaa
524 al3(i) = max(zero,min(one,al3(i)))
525 aaa = one/max(em30,x04(i)*x04(i)+y04(i)*y04(i)+z04(i)*z04(i))
526 hlc4(i)= lc4(i)*abs(lc4(i))*aaa
527 hlb3(i)= lb3(i)*abs(lb3(i))*aaa
528 al4(i) = -(xi0*x04(i)+yi0*y04(i)+zi0*z04(i))*aaa
529 al4(i) = max(zero,min(one,al4(i)))
530C
531 ENDDO
532C
533 DO i=lft,llt
534 x12 = x2(i) - x1(i)
535 y12 = y2(i) - y1(i)
536 z12 = z2(i) - z1(i)
537 la = one - lb1(i) - lc1(i)
538C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
539 aaa = one / max(em20,x12*x12+y12*y12+z12*z12)
540 hla= la*abs(la) * aaa
541 IF(la<zero.AND.
542 + hla<=hlb1(i).AND.hla<=hlc1(i))THEN
543 lb1(i) = (xi2(i)*x12+yi2(i)*y12+zi2(i)*z12) * aaa
544 lb1(i) = max(zero,min(one,lb1(i)))
545 lc1(i) = one - lb1(i)
546 ELSEIF(lb1(i)<zero.AND.
547 + hlb1(i)<=hlc1(i).AND.hlb1(i)<=hla)THEN
548 lb1(i) = zero
549 lc1(i) = al2(i)
550 ELSEIF(lc1(i)<zero.AND.
551 + hlc1(i)<=hla.AND.hlc1(i)<=hlb1(i))THEN
552 lc1(i) = zero
553 lb1(i) = al1(i)
554 ENDIF
555 ENDDO
556C
557 DO i=lft,llt
558 x23 = x3(i) - x2(i)
559 y23 = y3(i) - y2(i)
560 z23 = z3(i) - z2(i)
561 la = one - lb2(i) - lc2(i)
562C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
563 aaa = one / max(em20,x23*x23+y23*y23+z23*z23)
564 hla= la*abs(la) * aaa
565 IF(la<zero.AND.
566 + hla<=hlb2(i).AND.hla<=hlc2(i))THEN
567 lb2(i) = (xi3(i)*x23+yi3(i)*y23+zi3(i)*z23)*aaa
568 lb2(i) = max(zero,min(one,lb2(i)))
569 lc2(i) = one - lb2(i)
570 ELSEIF(lb2(i)<zero.AND.
571 + hlb2(i)<=hlc2(i).AND.hlb2(i)<=hla)THEN
572 lb2(i) = zero
573 lc2(i) = al3(i)
574 ELSEIF(lc2(i)<zero.AND.
575 + hlc2(i)<=hla.AND.hlc2(i)<=hlb2(i))THEN
576 lc2(i) = zero
577 lb2(i) = al2(i)
578 ENDIF
579 ENDDO
580C
581 DO i=lft,llt
582 x34 = x4(i) - x3(i)
583 y34 = y4(i) - y3(i)
584 z34 = z4(i) - z3(i)
585 la = one - lb3(i) - lc3(i)
586C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
587 aaa = one / max(em20,x34*x34+y34*y34+z34*z34)
588 hla= la*abs(la) * aaa
589 IF(la<zero.AND.
590 + hla<=hlb3(i).AND.hla<=hlc3(i))THEN
591 lb3(i) = (xi4(i)*x34+yi4(i)*y34+zi4(i)*z34)*aaa
592 lb3(i) = max(zero,min(one,lb3(i)))
593 lc3(i) = one - lb3(i)
594 ELSEIF(lb3(i)<zero.AND.
595 + hlb3(i)<=hlc3(i).AND.hlb3(i)<=hla)THEN
596 lb3(i) = zero
597 lc3(i) = al4(i)
598 ELSEIF(lc3(i)<zero.AND.
599 + hlc3(i)<=hla.AND.hlc3(i)<=hlb3(i))THEN
600 lc3(i) = zero
601 lb3(i) = al3(i)
602 ENDIF
603 ENDDO
604C
605 DO i=lft,llt
606 x41 = x1(i) - x4(i)
607 y41 = y1(i) - y4(i)
608 z41 = z1(i) - z4(i)
609 la = one - lb4(i) - lc4(i)
610C HLA, HLB1, HLC1 necessaires pour triangle angle obtu
611 aaa = one / max(em20,x41*x41+y41*y41+z41*z41)
612 hla= la*abs(la) * aaa
613 IF(la<zero.AND.
614 + hla<=hlb4(i).AND.hla<=hlc4(i))THEN
615 lb4(i) = (xi1(i)*x41+yi1(i)*y41+zi1(i)*z41)*aaa
616 lb4(i) = max(zero,min(one,lb4(i)))
617 lc4(i) = one - lb4(i)
618 ELSEIF(lb4(i)<zero.AND.
619 + hlb4(i)<=hlc4(i).AND.hlb4(i)<=hla)THEN
620 lb4(i) = zero
621 lc4(i) = al1(i)
622 ELSEIF(lc4(i)<zero.AND.
623 + hlc4(i)<=hla.AND.hlc4(i)<=hlb4(i))THEN
624 lc4(i) = zero
625 lb4(i) = al4(i)
626 ENDIF
627 ENDDO
628C
629 DO i=lft,llt
630C
631 nx1(i) = xi(i)-(x0(i) + lb1(i)*x01(i) + lc1(i)*x02(i))
632 ny1(i) = yi(i)-(y0(i) + lb1(i)*y01(i) + lc1(i)*y02(i))
633 nz1(i) = zi(i)-(z0(i) + lb1(i)*z01(i) + lc1(i)*z02(i))
634 p1(i) = nx1(i)*nx1(i) + ny1(i)*ny1(i) +nz1(i)*nz1(i)
635C
636 nx2(i) = xi(i)-(x0(i) + lb2(i)*x02(i) + lc2(i)*x03(i))
637 ny2(i) = yi(i)-(y0(i) + lb2(i)*y02(i) + lc2(i)*y03(i))
638 nz2(i) = zi(i)-(z0(i) + lb2(i)*z02(i) + lc2(i)*z03(i))
639 p2(i) = nx2(i)*nx2(i) + ny2(i)*ny2(i) +nz2(i)*nz2(i)
640C
641 nx3(i) = xi(i)-(x0(i) + lb3(i)*x03(i) + lc3(i)*x04(i))
642 ny3(i) = yi(i)-(y0(i) + lb3(i)*y03(i) + lc3(i)*y04(i))
643 nz3(i) = zi(i)-(z0(i) + lb3(i)*z03(i) + lc3(i)*z04(i))
644 p3(i) = nx3(i)*nx3(i) + ny3(i)*ny3(i) +nz3(i)*nz3(i)
645C
646 nx4(i) = xi(i)-(x0(i) + lb4(i)*x04(i) + lc4(i)*x01(i))
647 ny4(i) = yi(i)-(y0(i) + lb4(i)*y04(i) + lc4(i)*y01(i))
648 nz4(i) = zi(i)-(z0(i) + lb4(i)*z04(i) + lc4(i)*z01(i))
649 p4(i) = nx4(i)*nx4(i) + ny4(i)*ny4(i) +nz4(i)*nz4(i)
650C
651 ENDDO
652C
653 RETURN
654 END
655!||====================================================================
656!|| i20cgap0 ../starter/source/interfaces/inter3d1/i20dst3.F
657!||--- called by ------------------------------------------------------
658!|| i20dst3 ../starter/source/interfaces/inter3d1/i20dst3.F
659!||====================================================================
660 SUBROUTINE i20cgap0(X10,Y10,Z10,X20,Y20,Z20,
661 . X30,Y30,Z30,X40,Y40,Z40,
662 . X1,Y1,Z1,X2,Y2,Z2,
663 . X3,Y3,Z3,X4,Y4,Z4,GAP_M)
664C-----------------------------------------------
665C I m p l i c i t T y p e s
666C-----------------------------------------------
667#include "implicit_f.inc"
668C-----------------------------------------------
669C D u m m y A r g u m e n t s
670C-----------------------------------------------
671 my_real
672 . x10,y10,z10,x20,y20,z20,x30,y30,z30,x40,y40,z40,
673 . x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4,gap_m
674C-----------------------------------------------
675C L o c a l V a r i a b l e s
676C-----------------------------------------------
677 my_real
678 . l12,l12x,l12y,l12z,
679 . l23,l23x,l23y,l23z,l14,l14x,l14y,l14z,aaa,bbb,l1214,l1223
680C=======================================================================
681 l12x = x20-x10
682 l12y = y20-y10
683 l12z = z20-z10
684 l12 = sqrt(l12x*l12x+l12y*l12y+l12z*l12z)
685 aaa = one/max(l12,em30)
686 l12x = l12x*aaa
687 l12y = l12y*aaa
688 l12z = l12z*aaa
689
690 l14x = x40-x10
691 l14y = y40-y10
692 l14z = z40-z10
693 l14 = sqrt(l14x*l14x+l14y*l14y+l14z*l14z)
694 aaa = one/max(l14,em30)
695 l14x = l14x*aaa
696 l14y = l14y*aaa
697 l14z = l14z*aaa
698
699 l23x = x30-x20
700 l23y = y30-y20
701 l23z = z30-z20
702 l23 = sqrt(l23x*l23x+l23y*l23y+l23z*l23z)
703 aaa = one/max(l23,em30)
704 l23x = l23x*aaa
705 l23y = l23y*aaa
706 l23z = l23z*aaa
707
708 l1214 = l12x*l14x+l12y*l14y+l12z*l14z
709
710 bbb = gap_m*(one+em5)/max(sqrt(one-l1214*l1214),em30)
711 aaa = min(l12/three,bbb)
712
713 x1 = x1 + l14x*aaa
714 y1 = y1 + l14y*aaa
715 z1 = z1 + l14z*aaa
716
717 l1223 = l12x*l23x+l12y*l23y+l12z*l23z
718
719 bbb = gap_m*(one+em5)/max(sqrt(one-l1223*l1223),em30)
720 aaa = min(l14/three,bbb)
721
722 x2 = x2 + l23x*aaa
723 y2 = y2 + l23y*aaa
724 z2 = z2 + l23z*aaa
725
726 RETURN
727 END
728!||====================================================================
729!|| i20cgap1 ../starter/source/interfaces/inter3d1/i20dst3.F
730!||--- called by ------------------------------------------------------
731!|| i20dst3 ../starter/source/interfaces/inter3d1/i20dst3.F
732!||====================================================================
733 SUBROUTINE i20cgap1(X10,Y10,Z10,X20,Y20,Z20,
734 . X40,Y40,Z40,
735 . X1,Y1,Z1,X2,Y2,Z2,
736 . X4,Y4,Z4,GAP_M)
737C-----------------------------------------------
738C I m p l i c i t T y p e s
739C-----------------------------------------------
740#include "implicit_f.inc"
741C-----------------------------------------------
742C D u m m y A r g u m e n t s
743C-----------------------------------------------
744 my_real
745 . x10,y10,z10,x20,y20,z20,x40,y40,z40,
746 . x1,y1,z1,x2,y2,z2,x4,y4,z4,gap_m
747C-----------------------------------------------
748C L o c a l V a r i a b l e s
749C-----------------------------------------------
750 my_real
751 . l12,l12x,l12y,l12z,l14,l14x,l14y,l14z,aaa,bbb,l1214
752C=======================================================================
753 l12x = x20-x10
754 l12y = y20-y10
755 l12z = z20-z10
756 l12 = sqrt(l12x*l12x+l12y*l12y+l12z*l12z)
757 aaa = one/max(l12,em30)
758 l12x = l12x*aaa
759 l12y = l12y*aaa
760 l12z = l12z*aaa
761
762 l14x = x40-x10
763 l14y = y40-y10
764 l14z = z40-z10
765 l14 = sqrt(l14x*l14x+l14y*l14y+l14z*l14z)
766 aaa = one/max(l14,em30)
767 l14x = l14x*aaa
768 l14y = l14y*aaa
769 l14z = l14z*aaa
770
771 l1214 = l12x*l14x+l12y*l14y+l12z*l14z
772
773 bbb = gap_m*(one+em5)/max(sqrt(one-l1214*l1214),em30)
774 aaa = min(l12/three,bbb)
775
776 x1 = x1 + l12x*aaa
777 y1 = y1 + l12y*aaa
778 z1 = z1 + l12z*aaa
779
780 aaa = min(l14/three,bbb)
781
782 x1 = x1 + l14x*aaa
783 y1 = y1 + l14y*aaa
784 z1 = z1 + l14z*aaa
785
786 RETURN
787 END
788!||====================================================================
789!|| i20gap1 ../starter/source/interfaces/inter3d1/i20dst3.F
790!||--- called by ------------------------------------------------------
791!|| i20ini3 ../starter/source/interfaces/inter3d1/i20ini3.F
792!||--- calls -----------------------------------------------------
793!|| bitget ../starter/source/interfaces/inter3d1/bitget.F
794!|| bitset ../starter/source/interfaces/inter3d1/bitget.F
795!|| bitunset ../starter/source/interfaces/inter3d1/bitget.F
796!||====================================================================
797 SUBROUTINE i20gap1(NRTM,NSN,NLN,GAP_M,GAP_SH,
798 2 GAP_S,NBINFLG,NSV,NLG,TAG)
799C-----------------------------------------------
800C I m p l i c i t T y p e s
801C-----------------------------------------------
802#include "implicit_f.inc"
803C-----------------------------------------------
804C D u m m y A r g u m e n t s
805C-----------------------------------------------
806 INTEGER NRTM,NSN,NLN,NBINFLG(NLN),NSV(NSN),NLG(NLN),TAG(NUMNOD)
807 my_real
808 . GAP_SH(NRTM),GAP_M(NRTM),GAP_S(NSN)
809C-----------------------------------------------
810C C o m m o n B l o c k s
811C-----------------------------------------------
812#include "com04_c.inc"
813C-----------------------------------------------
814C L o c a l V a r i a b l e s
815C-----------------------------------------------
816 INTEGER I,J
817
818 INTEGER BITUNSET,BITGET,BITSET
819 EXTERNAL BITUNSET,BITGET,BITSET
820C=======================================================================
821
822 DO I=1,nrtm
823 gap_m(i) = gap_m(i)-two*gap_sh(i)
824 ENDDO
825c mise a zero du gap secnd sur le bord des coques
826c ici NVS donne encore le noeud global, NSV sera modifie dans I20NLG
827 DO i=1,nln
828 j = nlg(i)
829 tag(j)=i
830 ENDDO
831 DO i=1,nsn
832 IF(bitget(nbinflg(tag(nsv(i))),7)/=0)THEN
833 gap_s(i) = zero
834 ENDIF
835 ENDDO
836C
837 RETURN
838 END
839!||====================================================================
840!|| i20norm ../starter/source/interfaces/inter3d1/i20dst3.F
841!||--- called by ------------------------------------------------------
842!|| i20ini3 ../starter/source/interfaces/inter3d1/i20ini3.F
843!||====================================================================
844 SUBROUTINE i20norm(NRTM,IRECT,NUMNOD,X,SOLIDN_NORMAL,
845 . NMN ,MSR ,NLN,NLG,GAP_SH)
846C-----------------------------------------------
847C I m p l i c i t T y p e s
848C-----------------------------------------------
849#include "implicit_f.inc"
850C-----------------------------------------------
851C D u m m y A r g u m e n t s
852C-----------------------------------------------
853 INTEGER NRTM,NUMNOD,IRECT(4,NRTM),NMN,MSR(*),NLN,NLG(NLN)
854C REAL
855 my_real
856 . x(3,numnod), solidn_normal(3,numnod),gap_sh(*)
857C-----------------------------------------------
858C L o c a l V a r i a b l e s
859C-----------------------------------------------
860 INTEGER I ,J ,N1,N2,N3,N4
861 my_real
862 . SURFX,SURFY,SURFZ,X12,Y12,Z12,X23,Y23,Z23,
863 . X34,Y34,Z34,X41,Y41,Z41,AAA,SI,CO
864C-----------------------------------------------
865
866 DO i=1,nmn
867c N1 = NLG(MSR(I))
868 n1 = msr(i)
869 solidn_normal(1,n1) = zero
870 solidn_normal(2,n1) = zero
871 solidn_normal(3,n1) = zero
872 END DO
873
874 DO i=1,nrtm
875c test pour ne prendre en compte que les solides
876 IF(gap_sh(i)/=zero)THEN
877c N1 = NLG(IRECT(1,I))
878c N2 = NLG(IRECT(2,I))
879c N3 = NLG(IRECT(3,I))
880c N4 = NLG(IRECT(4,I))
881 n1 = irect(1,i)
882 n2 = irect(2,i)
883 n3 = irect(3,i)
884 n4 = irect(4,i)
885
886 x41 = x(1,n1) - x(1,n4)
887 y41 = x(2,n1) - x(2,n4)
888 z41 = x(3,n1) - x(3,n4)
889c AAA=ONE/MAX(EM30,SQRT(X41*X41+Y41*Y41+Z41*Z41))
890c X41 = X41 * AAA
891c Y41 = Y41 * AAA
892c Z41 = Z41 * AAA
893
894 x12 = x(1,n2) - x(1,n1)
895 y12 = x(2,n2) - x(2,n1)
896 z12 = x(3,n2) - x(3,n1)
897c AAA=ONE/MAX(EM30,SQRT(X12*X12+Y12*Y12+Z12*Z12))
898c X12 = X12 * AAA
899c Y12 = Y12 * AAA
900c Z12 = Z12 * AAA
901
902 x23 = x(1,n3) - x(1,n2)
903 y23 = x(2,n3) - x(2,n2)
904 z23 = x(3,n3) - x(3,n2)
905c AAA=ONE/MAX(EM30,SQRT(X23*X23+Y23*Y23+Z23*Z23))
906c X23 = X23 * AAA
907c Y23 = Y23 * AAA
908c Z23 = Z23 * AAA
909
910 surfx = y41*z12 - z41*y12
911 surfy = z41*x12 - x41*z12
912 surfz = x41*y12 - y41*x12
913
914 co = x41*x12 + y41*y12 + z41*z12
915 si = sqrt(surfx*surfx+surfy*surfy+surfz*surfz)
916 aaa = atan2(si,-co)*two/pi/max(em30,si)
917 solidn_normal(1,n1) = solidn_normal(1,n1) + surfx*aaa
918 solidn_normal(2,n1) = solidn_normal(2,n1) + surfy*aaa
919 solidn_normal(3,n1) = solidn_normal(3,n1) + surfz*aaa
920
921 surfx = y12*z23 - z12*y23
922 surfy = z12*x23 - x12*z23
923 surfz = x12*y23 - y12*x23
924 co = x12*x23 + y12*y23 + z12*z23
925 si = sqrt(surfx*surfx+surfy*surfy+surfz*surfz)
926 aaa = atan2(si,-co)*two/pi/max(em30,si)
927 solidn_normal(1,n2) = solidn_normal(1,n2) + surfx*aaa
928 solidn_normal(2,n2) = solidn_normal(2,n2) + surfy*aaa
929 solidn_normal(3,n2) = solidn_normal(3,n2) + surfz*aaa
930
931 IF(n3 /= n4)THEN
932 x34 = x(1,n4) - x(1,n3)
933 y34 = x(2,n4) - x(2,n3)
934 z34 = x(3,n4) - x(3,n3)
935c AAA=ONE/MAX(EM30,SQRT(X34*X34+Y34*Y34+Z34*Z34))
936c X34 = X34 * AAA
937c Y34 = Y34 * AAA
938c Z34 = Z34 * AAA
939 surfx = y23*z34 - z23*y34
940 surfy = z23*x34 - x23*z34
941 surfz = x23*y34 - y23*x34
942 co = x23*x34 + y23*y34 + z23*z34
943 si = sqrt(surfx*surfx+surfy*surfy+surfz*surfz)
944 aaa = atan2(si,-co)*two/pi/max(em30,si)
945 solidn_normal(1,n3) = solidn_normal(1,n3) + surfx*aaa
946 solidn_normal(2,n3) = solidn_normal(2,n3) + surfy*aaa
947 solidn_normal(3,n3) = solidn_normal(3,n3) + surfz*aaa
948
949 surfx = y34*z41 - z34*y41
950 surfy = z34*x41 - x34*z41
951 surfz = x34*y41 - y34*x41
952 co = x34*x41 + y34*y41 + z34*z41
953 si = sqrt(surfx*surfx+surfy*surfy+surfz*surfz)
954 aaa = atan2(si,-co)*two/pi/max(em30,si)
955 solidn_normal(1,n4) = solidn_normal(1,n4) + surfx*aaa
956 solidn_normal(2,n4) = solidn_normal(2,n4) + surfy*aaa
957 solidn_normal(3,n4) = solidn_normal(3,n4) + surfz*aaa
958 ELSE
959 surfx = y23*z41 - z23*y41
960 surfy = z23*x41 - x23*z41
961 surfz = x23*y41 - y23*x41
962 co = x41*x23 + y41*y23 + z41*z23
963 si = sqrt(surfx*surfx+surfy*surfy+surfz*surfz)
964 aaa = atan2(si,-co)*two/pi/max(em30,si)
965 solidn_normal(1,n3) = solidn_normal(1,n3) + surfx*aaa
966 solidn_normal(2,n3) = solidn_normal(2,n3) + surfy*aaa
967 solidn_normal(3,n3) = solidn_normal(3,n3) + surfz*aaa
968 ENDIF
969 ENDIF
970 ENDDO
971
972c DO I=1,NMN
973cc N1 = NLG(MSR(I))
974c N1 = MSR(I)
975c AAA = SOLIDN_NORMAL(1,N1)*SOLIDN_NORMAL(1,N1)
976c + + SOLIDN_NORMAL(2,N1)*SOLIDN_NORMAL(2,N1)
977c + + SOLIDN_NORMAL(3,N1)*SOLIDN_NORMAL(3,N1)
978c AAA = ONE/MAX(EM20,SQRT(AAA))
979c SOLIDN_NORMAL(1,N1) = SOLIDN_NORMAL(1,N1)*AAA
980c SOLIDN_NORMAL(2,N1) = SOLIDN_NORMAL(2,N1)*AAA
981c SOLIDN_NORMAL(3,N1) = SOLIDN_NORMAL(3,N1)*AAA
982c END DO
983
984 RETURN
985 END
986!||====================================================================
987!|| i20dst3e ../starter/source/interfaces/inter3d1/i20dst3.F
988!||--- called by ------------------------------------------------------
989!|| i20ini3 ../starter/source/interfaces/inter3d1/i20ini3.F
990!||====================================================================
991 SUBROUTINE i20dst3e(JLT ,GAP ,CAND_S,CAND_M,IRECTS ,
992 . IRECTM ,NX,NY,NZ,
993 . N1,N2,M1,M2,JLT_NEW,
994 . X ,IGAP,GAP_S,GAP_M, GAPV2,
995 . NLN ,NLG ,SOLIDN_NORMAL )
996C-----------------------------------------------
997C I m p l i c i t T y p e s
998C-----------------------------------------------
999#include "implicit_f.inc"
1000C-----------------------------------------------
1001C G l o b a l P a r a m e t e r s
1002C-----------------------------------------------
1003#include "mvsiz_p.inc"
1004C-----------------------------------------------
1005C D u m m y A r g u m e n t s
1006C-----------------------------------------------
1007 INTEGER JLT,JLT_NEW,IGAP
1008 INTEGER IRECTS(2,*), IRECTM(2,*),CAND_S(*),CAND_M(*),
1009 . n1(*),n2(*),m1(*),m2(*)
1010 my_real
1011 . gap,
1012 . nx(*),ny(*),nz(*),x(3,*),gap_s(*),gap_m(*), gapv2(*)
1013C-----------------------------------------------
1014C L o c a l V a r i a b l e s
1015C-----------------------------------------------
1016 INTEGER I,NLN,NLG(NLN)
1017 my_real
1018 . XS12,YS12,ZS12,XM12,YM12,ZM12,XA,XB,
1019 . XS2,XM2,XSM,XS2M2,YS2,YM2,YSM,YS2M2,ZS2,ZM2,ZSM,ZS2M2,
1020 . XX,YY,ZZ,ALS,ALM,DET,H1S,H2S,H1M,H2M,
1021 . GAP2,XXS1,YYS1,ZZS1,XXS2,YYS2,ZZS2,
1022 . XXM1,YYM1,ZZM1,XXM2,YYM2,ZZM2,
1023 . PENE2(MVSIZ),SOLIDN_NORMAL(3,*),AAA,
1024 . SX1,SX2,SX3,SX4,SY1,SY2,SY3,SY4,SZ1,SZ2,SZ3,SZ4
1025
1026C-----------------------------------------------
1027 GAP2=gap*gap
1028C
1029 IF(igap==0)THEN
1030 DO i=1,jlt
1031 gapv2(i)=gap2
1032 ENDDO
1033 ELSE
1034 DO i=1,jlt
1035 gapv2(i)=gap_s(cand_s(i))+gap_m(cand_m(i))
1036 gapv2(i)=max(gapv2(i)*gapv2(i),gap2)
1037 ENDDO
1038 ENDIF
1039C--------------------------------------------------------
1040C
1041C--------------------------------------------------------
1042C F = [A*X1+(1-A)*X2-B*X3-(1-B)*X4]^2 + [..Y..]^2 + [..Z..]^2
1043C DF/DA = 0 = (X1-X2)(A(X1-X2)+X2-X4 +B(X4-X3))+...
1044C DF/DA = 0 = A(X1-X2)^2 +X2-X4 + B(X1-X2)(X4-X3))+...
1045C DF/DA = 0 = A[(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
1046C + B[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
1047C + (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
1048C DF/DB = 0 = (X4-X3)(A(X1-X2)+X2-X4 +B(X4-X3))+...
1049C DF/DB = 0 = B[(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
1050C + A[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
1051C + (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
1052C XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
1053C XM2 = [(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
1054C XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
1055C XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
1056C XB = (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
1057C A XS2 + B XSM + XA = 0
1058C A XSM + B XM2 + XB = 0
1059C
1060C A = -(XA + B XSM)/XS2
1061C -(XA + B XSM)*XSM + B XM2*XS2 + XB*XS2 = 0
1062C -B XSM*XSM + B XM2*XS2 + XB*XS2-XA*XSM = 0
1063C B*(XM2*XS2 - XSM*XSM) = -XB*XS2+XA*XSM
1064C B = (XA*XSM-XB*XS2) / (XM2*XS2 - XSM*XSM)
1065C A = (XB*XSM-XA*XM2) / (XM2*XS2 - XSM*XSM)
1066 DO i=1,jlt
1067 n1(i)=irects(1,cand_s(i))
1068 n2(i)=irects(2,cand_s(i))
1069 m1(i)=irectm(1,cand_m(i))
1070 m2(i)=irectm(2,cand_m(i))
1071 xxs1 = x(1,n1(i))
1072 yys1 = x(2,n1(i))
1073 zzs1 = x(3,n1(i))
1074 xxs2 = x(1,n2(i))
1075 yys2 = x(2,n2(i))
1076 zzs2 = x(3,n2(i))
1077 xxm1 = x(1,m1(i))
1078 yym1 = x(2,m1(i))
1079 zzm1 = x(3,m1(i))
1080 xxm2 = x(1,m2(i))
1081 yym2 = x(2,m2(i))
1082 zzm2 = x(3,m2(i))
1083
1084 IF(igap/=0)THEN
1085 aaa = gap_s(cand_s(i))
1086 sx1 = solidn_normal(1,n1(i))*aaa
1087 sy1 = solidn_normal(2,n1(i))*aaa
1088 sz1 = solidn_normal(3,n1(i))*aaa
1089 sx2 = solidn_normal(1,n2(i))*aaa
1090 sy2 = solidn_normal(2,n2(i))*aaa
1091 sz2 = solidn_normal(3,n2(i))*aaa
1092 aaa = gap_m(cand_m(i))
1093 sx3 = solidn_normal(1,m1(i))*aaa
1094 sy3 = solidn_normal(2,m1(i))*aaa
1095 sz3 = solidn_normal(3,m1(i))*aaa
1096 sx4 = solidn_normal(1,m2(i))*aaa
1097 sy4 = solidn_normal(2,m2(i))*aaa
1098 sz4 = solidn_normal(3,m2(i))*aaa
1099 xxs1 = xxs1 - sx1
1100 yys1 = yys1 - sy1
1101 zzs1 = zzs1 - sz1
1102 xxs2 = xxs2 - sx2
1103 yys2 = yys2 - sy2
1104 zzs2 = zzs2 - sz2
1105 xxm1 = xxm1 - sx3
1106 yym1 = yym1 - sy3
1107 zzm1 = zzm1 - sz3
1108 xxm2 = xxm2 - sx4
1109 yym2 = yym2 - sy4
1110 zzm2 = zzm2 - sz4
1111 ENDIF
1112
1113 xs12 = xxs2-xxs1
1114 ys12 = yys2-yys1
1115 zs12 = zzs2-zzs1
1116 xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
1117 xm12 = xxm2-xxm1
1118 ym12 = yym2-yym1
1119 zm12 = zzm2-zzm1
1120 xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
1121 xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
1122
1123 xs2m2 = xxm2-xxs2
1124 ys2m2 = yym2-yys2
1125 zs2m2 = zzm2-zzs2
1126
1127 xa = xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
1128 xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2
1129 det = xm2*xs2 - xsm*xsm
1130 det = max(em20,det)
1131C
1132 h1m = (xa*xsm-xb*xs2) / det
1133C
1134 xs2 = max(xs2,em20)
1135 xm2 = max(xm2,em20)
1136 h1m=min(one,max(zero,h1m))
1137 h1s = -(xa + h1m*xsm) / xs2
1138 h1s=min(one,max(zero,h1s))
1139 h1m = -(xb + h1s*xsm) / xm2
1140 h1m=min(one,max(zero,h1m))
1141C
1142 h2s = one - h1s
1143 h2m = one - h1m
1144C !!!!!!!!!!!!!!!!!!!!!!!
1145C PENE = GAP^2 - DIST^2 UTILISE POUR TESTER SI NON NUL
1146C!!!!!!!!!!!!!!!!!!!!!!!!
1147 nx(i) = h1s*xxs1 + h2s*xxs2
1148 . - h1m*xxm1 - h2m*xxm2
1149 ny(i) = h1s*yys1 + h2s*yys2
1150 . - h1m*yym1 - h2m*yym2
1151 nz(i) = h1s*zzs1 + h2s*zzs2
1152 . - h1m*zzm1 - h2m*zzm2
1153 pene2(i) = gapv2(i) - nx(i)*nx(i) - ny(i)*ny(i) - nz(i)*nz(i)
1154 pene2(i) = max(zero,pene2(i))
1155C
1156 ENDDO
1157 DO i=1,jlt
1158 IF(pene2(i)/=zero)THEN
1159 jlt_new = jlt_new + 1
1160 cand_s(jlt_new) = cand_s(i)
1161 cand_m(jlt_new) = cand_m(i)
1162 n1(jlt_new) = n1(i)
1163 n2(jlt_new) = n2(i)
1164 m1(jlt_new) = m1(i)
1165 m2(jlt_new) = m2(i)
1166 nx(jlt_new) = nx(i)
1167 ny(jlt_new) = ny(i)
1168 nz(jlt_new) = nz(i)
1169 gapv2(jlt_new) = gapv2(i)
1170 ENDIF
1171 ENDDO
1172C-----------
1173 RETURN
1174 END
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
i20dst3
subroutine i20dst3(igap, gap_sh, cand_e, cand_n, gapv, gap, gap_s, gap_m, gapmax, gapmin, irect, nln, nlg, solidn_normal, nsv, nbinflg, tag, ix3, ix4, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, x0, y0, z0, nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3, nx4, ny4, nz4, p1, p2, p3, p4, lb1, lb2, lb3, lb4, lc1, lc2, lc3, lc4)
Definition i20dst3.F:47
i20dst3e
subroutine i20dst3e(jlt, gap, cand_s, cand_m, irects, irectm, nx, ny, nz, n1, n2, m1, m2, jlt_new, x, igap, gap_s, gap_m, gapv2, nln, nlg, solidn_normal)
Definition i20dst3.F:996
i20gap1
subroutine i20gap1(nrtm, nsn, nln, gap_m, gap_sh, gap_s, nbinflg, nsv, nlg, tag)
Definition i20dst3.F:799
i20cgap1
subroutine i20cgap1(x10, y10, z10, x20, y20, z20, x40, y40, z40, x1, y1, z1, x2, y2, z2, x4, y4, z4, gap_m)
Definition i20dst3.F:737
i20norm
subroutine i20norm(nrtm, irect, numnod, x, solidn_normal, nmn, msr, nln, nlg, gap_sh)
Definition i20dst3.F:846
i20cgap0
subroutine i20cgap0(x10, y10, z10, x20, y20, z20, x30, y30, z30, x40, y40, z40, x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, gap_m)
Definition i20dst3.F:664