OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i24dst3.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!|| i24dst3 ../engine/source/interfaces/int24/i24dst3.F
25!||--- called by ------------------------------------------------------
26!|| i24mainf ../engine/source/interfaces/int24/i24main.F
27!||--- calls -----------------------------------------------------
28!|| i24nexttria ../engine/source/interfaces/int24/i24dst3.F
29!|| intersecsh ../engine/source/interfaces/int24/i24dst3.F
30!|| intersecv ../engine/source/interfaces/int24/i24dst3.F
31!|| intersecv0 ../engine/source/interfaces/int24/i24dst3.F
32!|| s_in_m4 ../engine/source/interfaces/int24/i24dst3.F
33!||--- uses -----------------------------------------------------
34!|| tri7box ../engine/share/modules/tri7box.F
35!||====================================================================
36 SUBROUTINE i24dst3(
37 1 JLT ,CAND_N ,CAND_E ,CN_LOC ,CE_LOC ,
38 2 X1 ,X2 ,X3 ,X4 ,Y1 ,
39 3 Y2 ,Y3 ,Y4 ,Z1 ,Z2 ,
40 4 Z3 ,Z4 ,XI ,YI ,ZI ,
41 5 VX1 ,VX2 ,VX3 ,VX4 ,VXI ,
42 6 VY1 ,VY2 ,VY3 ,VY4 ,VYI ,
43 7 VZ1 ,VZ2 ,VZ3 ,VZ4 ,VZI ,
44 8 N1 ,N2 ,N3 ,H1 ,H2 ,
45 9 H3 ,H4 ,NIN ,NSN ,IX1 ,
46 A IX2 ,IX3 ,IX4 ,NSVG ,STIF ,
47 B JLT_NEW,INACTI ,MSEGLO ,GAPS ,GAP_NM ,
48 C KINI ,IRECT ,IRTLM ,TIME_S ,SUBTRIA,
49 D INTTH ,NSMS ,PENE ,XX0 ,YY0 ,
50 E ZZ0 ,VX ,VY ,VZ ,IXX ,
51 F MVOISIN,PMAX_GAP,SECND_FR,GAP_M ,PENE_OLD,
52 G STIF_OLD,ITRIV ,ITAB ,CAND_T ,IEDGE ,
53 H NFT ,PENMIN,EPS0 ,NM1 ,NM2 ,
54 I NM3 ,INTPLY ,DGAP_NM,ICONT_I,MARGE ,
55 J NSNE ,ISPT2 ,IZERO ,IKNON ,PENREF )
56C-----------------------------------------------
57C M o d u l e s
58C-----------------------------------------------
59 USE tri7box
60C-----------------------------------------------
61C I m p l i c i t T y p e s
62C-----------------------------------------------
63#include "implicit_f.inc"
64C-----------------------------------------------
65C G l o b a l P a r a m e t e r s
66C-----------------------------------------------
67#include "mvsiz_p.inc"
68#include "impl1_c.inc"
69#include "comlock.inc"
70C-----------------------------------------------
71C D u m m y A r g u m e n t s
72C-----------------------------------------------
73 INTEGER JLT, JLT_NEW,NIN,NSN,INACTI,IEDGE, NFT,
74 . CAND_N(*),CN_LOC(MVSIZ),
75 . CAND_E(*),CE_LOC(MVSIZ), NSVG(MVSIZ), KINI(*), CAND_T(*)
76 INTEGER IX1(MVSIZ), IX2(MVSIZ), IX3(MVSIZ), IX4(MVSIZ),
77 . INTTH,MSEGLO(*),IRTLM(2,NSN) ,SUBTRIA(MVSIZ),
78 . nsms(*),ixx(mvsiz,13),mvoisin(4,*),itriv(4,mvsiz),
79 . intply ,nsne,ispt2(*),izero
80 my_real
81 . n1(mvsiz), n2(mvsiz), n3(mvsiz), pmax_gap,
82 . h1(mvsiz), h2(mvsiz), h3(mvsiz), h4(mvsiz),
83 . kb1(mvsiz), kb2(mvsiz), kb3(mvsiz), kb4(mvsiz),
84 . kc1(mvsiz), kc2(mvsiz), kc3(mvsiz), kc4(mvsiz),
85 . x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz),
86 . y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
87 . z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),
88 . xi(mvsiz), yi(mvsiz), zi(mvsiz), stif(mvsiz),
89 . vx1(mvsiz),vx2(mvsiz),vx3(mvsiz),vx4(mvsiz),vy1(mvsiz),
90 . vy2(mvsiz),vy3(mvsiz),vy4(mvsiz),vz1(mvsiz),vz2(mvsiz),
91 . vz3(mvsiz),vz4(mvsiz),vxi(mvsiz),vyi(mvsiz),vzi(mvsiz),
92 . time_s(nsn),pene(mvsiz),secnd_fr(6,nsn),
93 . xx0(mvsiz,17),yy0(mvsiz,17),zz0(mvsiz,17),gap_m(*),
94 . vx(mvsiz,17),vy(mvsiz,17),vz(mvsiz,17),gaps(*),gap_nm(12,*),
95 . pene_old(5,nsn),stif_old(2,nsn),penmin,eps0,
96 . nm1(mvsiz), nm2(mvsiz), nm3(mvsiz),dgap_nm(4,*),marge
97 INTEGER IRECT(4,*),ITAB(*),ICONT_I(*)
98 INTEGER, DIMENSION(MVSIZ), INTENT(IN) :: IKNON
99 my_real, DIMENSION(MVSIZ), INTENT(OUT) :: penref
100C-----------------------------------------------
101C C o m m o n B l o c k s
102C-----------------------------------------------
103#include "scr05_c.inc"
104#include "com08_c.inc"
105C-----------------------------------------------
106C L o c a l V a r i a b l e s
107C-----------------------------------------------
108 INTEGER I, IG, J, K, L, INTERSECT,MG,N,II,NSLID,ITR,ITS,ITT,ITQ,NS
109 my_real
110 . X5(MVSIZ), Y5(MVSIZ), Z5(MVSIZ),
111 . VX5(MVSIZ), VY5(MVSIZ), VZ5(MVSIZ),
112 . AL1(MVSIZ), AL2(MVSIZ), AL3(MVSIZ), AL4(MVSIZ),
113 . NXX(MVSIZ,17),NYY(MVSIZ,17),NZZ(MVSIZ,17),
114 . X51, X52, X53, X54,
115 . Y51, Y52, Y53, Y54,
116 . z51, z52, z53, z54,
117 . xo16, xo17, xo15, xo14, x163, x153, x152, x142, x141,
118 . yo16, yo17, yo15, yo14, y163, y153, y152, y142, y141,
119 . zo16, zo17, zo15, zo14, z163, z153, z152, z142, z141,
120 . x164, x174, x171,
121 . y164, y174, y171,
122 . z164, z174, z171,
123 . xi1, xi2, xi3, xi4, xi5,
124 . yi1, yi2, yi3, yi4, yi5,
125 . zi1, zi2, zi3, zi4, zi5,
126 . xia, xib, xic,
127 . yia, yib, yic,
128 . zia, zib, zic,
129 . xo1, xo2, xo3, xo4, xo5, xoi,
130 . yo1, yo2, yo3, yo4, yo5, yoi,
131 . zo1, zo2, zo3, zo4, zo5, zoi,xs,ys,zs,
132 . xm1, xm2, xm3, xm4, xm5,
133 . ym1, ym2, ym3, ym4, ym5,
134 . zm1, zm2, zm3, zm4, zm5
135 my_real
136 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
137 . sx125(mvsiz),sx235(mvsiz),sx345(mvsiz),sx415(mvsiz),
138 . sy125(mvsiz),sy235(mvsiz),sy345(mvsiz),sy415(mvsiz),
139 . sz125(mvsiz),sz235(mvsiz),sz345(mvsiz),sz415(mvsiz),
140 . sx1250(mvsiz),sx2350(mvsiz),sx3450(mvsiz),sx4150(mvsiz),
141 . sy1250(mvsiz),sy2350(mvsiz),sy3450(mvsiz),sy4150(mvsiz),
142 . sz1250(mvsiz),sz2350(mvsiz),sz3450(mvsiz),sz4150(mvsiz),
143 . sx2114(mvsiz),sx3215(mvsiz),sx4316(mvsiz),sx1417(mvsiz),
144 . sy2114(mvsiz),sy3215(mvsiz),sy4316(mvsiz),sy1417(mvsiz),
145 . sz2114(mvsiz),sz3215(mvsiz),sz4316(mvsiz),sz1417(mvsiz),
146 . sx21140(mvsiz),sx32150(mvsiz),sx43160(mvsiz),sx14170(mvsiz),
147 . sy21140(mvsiz),sy32150(mvsiz),sy43160(mvsiz),sy14170(mvsiz),
148 . sz21140(mvsiz),sz32150(mvsiz),sz43160(mvsiz),sz14170(mvsiz),
149 . la(mvsiz),lb(mvsiz), lc(mvsiz),
150 . xa(mvsiz),ya(mvsiz),za(mvsiz),
151 . xb(mvsiz),yb(mvsiz),zb(mvsiz),
152 . xc(mvsiz),yc(mvsiz),zc(mvsiz),
153 . xd(mvsiz),yd(mvsiz),zd(mvsiz),
154 . xe(mvsiz),ye(mvsiz),ze(mvsiz),
155 . xf(mvsiz),yf(mvsiz),zf(mvsiz),adt0(mvsiz),
156 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),
157 . xxl(17),yyl(17),zzl(17)
158 my_real
159 . s2,d1,d2,d3,d4,unp,zerom,h0,nqx,nqy,nqz,la2,lb2,lc2,
160 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
161 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
162 . sz,sz1,sz2,sz3,sz4,sz5,sz0,
163 . v12,v23,v34,v41,
164 . nx1, nx2, nx3, nx4,
165 . ny1, ny2, ny3, ny4,
166 . nz1, nz2, nz3, nz4,
167 . nx14, nx15, nx16, nx17,
168 . ny14, ny15, ny16, ny17,
169 . nz14, nz15, nz16, nz17,
170 . sox125,sox235,sox345,sox415,
171 . soy125,soy235,soy345,soy415,
172 . soz125,soz235,soz345,soz415,
173 . sox2114,sox3215,sox4316,sox1417,
174 . soy2114,soy3215,soy4316,soy1417,
175 . soz2114,soz3215,soz4316,soz1417,
176 . xo12,xo23,xo34,xo41,sxo1,sxo2,sxo3,sxo4,sxo5,
177 . yo12,yo23,yo34,yo41,syo1,syo2,syo3,syo4,syo5,
178 . zo12,zo23,zo34,zo41,szo1,szo2,szo3,szo4,szo5,
179 . sxo12,sxo23,sxo34,sxo41,vo12,vo23,vo34,vo41,
180 . syo12,syo23,syo34,syo41,
181 . szo12,szo23,szo34,szo41,
182 . x13_2,y13_2,z13_2,x7_3 ,y7_3 ,z7_3 ,
183 . x9_4 ,y9_4 ,z9_4 ,x11_1,y11_1,z11_1,
184 . x6_4 ,y6_4 ,z6_4 ,x8_1 ,y8_1 ,z8_1 ,
185 . x10_2,y10_2,z10_2,x12_3,y12_3,z12_3,
186 . prx,pry,prz,psx,psy,psz,ptx,pty,ptz,
187 . nqrx,nqry,nqrz,nqsx,nqsy,nqsz,nqtx,nqty,nqtz,
188 . nrx,nry,nrz,nsx,nsy,nsz,ntx,nty,ntz,
189 . nax,nay,naz,nbx,nby,nbz,ncx,ncy,ncz,
190 . sax,say,saz,sbx,sby,sbz,scx,scy,scz,
191 . trx,try,trz,tsx,tsy,tsz,ttx,tty,ttz,
192 . tr2,ts2,tt2,aaa,bbb,vr,vs,vt,nnx,nny,nnz,
193 . xab,xbc,xca,yab,ybc,yca,zab,zbc,zca,
194 . xce,yce,zce,xaf,yaf,zaf,
195 . xpa,ypa,zpa,xpb,ypb,zpb,xpc,ypc,zpc,xpi,ypi,zpi,
196 . abx,bcx,cax,aby,bcy,cay,abz,bcz,caz,xbd,ybd,zbd
197 my_real
198 . gap2, ds2,t1,t2,t3,xxx,yyy,zzz,nni,ni2,
199 . al1num,al2num,al3num,al4num,al1den,al2den,al3den,al4den,
200 . x23d,y23d,z23d,x34d,y34d,z34d,x41d,y41d,z41d,
201 . x12d,y12d,z12d,gap2d,xi5d,yi5d,zi5d,s2d, la0,lb0,lc0,
202 . hla, adt,overw,eps,pene_sh,overw0,vol,eps2,dgap,rx,ry,rz,
203 . pene_tm(mvsiz),la_tm(mvsiz),lb_tm(mvsiz),n_tm(3,mvsiz),
204 . area(mvsiz),f_pene,fac_k,fac_p
205 INTEGER ICONTACT(MVSIZ),IT0(3,20),IT(3,20,MVSIZ),IC(10,20),
206 . ISTEP(MVSIZ),IRTLM_OLD(MVSIZ),SUBTRIA_N(MVSIZ),
207 . NSS,MGLOB,ibug,JG,IZLIM,ip,IKEEP,LARGEP,IPEN0(MVSIZ),IER,
208 . NOD1,NOD2,NOR_OLD,ISHEL(MVSIZ),ICONT_R(MVSIZ),
209 . SUBTRIA_TM(MVSIZ),IFIRST
210 my_real
211 . UNPT,ZEROMT,EPS02,TOLE,EPSEG,MARGE025,TOL_INTS,TOL_B,EPSEXT
212! +----+----------+----------+-------------+
213! | IC |
214! +----+----------+----------+-------------+
215! | ST | noeuds T | noeuds TV| noeuds Quad |
216! +----+----------+----------+-------------+ 11-------10
217! | 1 | 5 1 2 | 14 3 4 | 3 4 1 2 | |\ 19 /|
218! | 2 | 5 2 3 | 15 4 1 | 4 1 2 3 | | \ / |
219! | 3 | 5 3 4 | 16 1 2 | 1 2 3 4 | | \ / |
220! | 4 | 5 4 1 | 17 2 3 | 2 3 4 1 | | 16 |
221! +----+----------+----------+-------------+ |15/ \11|
222! | 5 | 14 2 1 | 5 6 7 | 6 7 2 1 | | / \ |
223! | 6 | 15 3 2 | 5 8 9 | 8 9 3 2 | |/ 7 \|
224! | 7 | 16 4 3 | 5 10 11 | 10 11 4 3 |12-------4-------3-------9
225! | 8 | 17 1 4 | 5 12 13 | 12 13 1 4 | |\ 12 /|\ /|\ 14 /|
226! +----+----------+----------+-------------+ | \ / | \ 3 / | \ / |
227! | 9 | 14 1 6 |(13) 7 2 | 7 2 1 6 | | \ / | \ /2 |6 \ /18|
228! | 10 | 15 2 8 |( 7) 9 3 | 9 3 2 8 | | 17 | 5 | 15 |
229! | 11 | 16 3 10 |( 9)11 4 | 11 4 3 10 | |20/ \ 8| 4/ \ | / \ |
230! | 12 | 17 4 12 |(11)13 1 | 13 1 4 12 | | / \ | / 1 \ | / \ |
231! +----+----------+----------+-------------+ |/ 16 \|/ \|/ 10 \|
232! | 13 | 14 7 2 |( 8) 1 6 | 1 6 7 2 |13-------1-------2-------8
233! | 14 | 15 9 3 |(10) 2 8 | 2 8 9 3 | |\ 5 /|
234! | 15 | 16 11 4 |(12) 3 10 | 3 10 11 4 | | \ / |
235! | 16 | 17 13 1 |( 6) 4 12 | 4 12 13 1 | |9 \ /13|
236! +----+----------+----------+-------------+ | 14 |
237! | 17 | 14 6 7 | - 2 1 | 2 1 6 7 | | / \ |
238! | 18 | 15 8 9 | - 3 2 | 3 2 8 9 | | / \ |
239! | 19 | 16 10 11 | - 4 3 | 4 3 10 11 | |/ 17 \|
240! | 20 | 17 12 13 | - 1 4 | 1 4 12 13 | 6-------7
241! +----+----------+----------+-------------+
242 DATA ic /
243 1 5, 1, 2 , 14, 3, 4 , 3, 4, 1, 2,
244 2 5, 2, 3 , 15, 4, 1 , 4, 1, 2, 3,
245 3 5, 3, 4 , 16, 1, 2 , 1, 2, 3, 4,
246 4 5, 4, 1 , 17, 2, 3 , 2, 3, 4, 1,
247 5 14, 2, 1 , 5, 6, 7 , 6, 7, 2, 1,
248 6 15, 3, 2 , 5, 8, 9 , 8, 9, 3, 2,
249 7 16, 4, 3 , 5,10,11 , 10,11, 4, 3,
250 8 17, 1, 4 , 5,12,13 , 12,13, 1, 4,
251 9 14, 1, 6 , 0, 7, 2 , 7, 2, 1, 6,
252 . 15, 2, 8 , 0, 9, 3 , 9, 3, 2, 8,
253 1 16, 3,10 , 0,11, 4 , 11, 4, 3,10,
254 2 17, 4,12 , 0,13, 1 , 13, 1, 4,12,
255 3 14, 7, 2 , 0, 1, 6 , 1, 6, 7, 2,
256 4 15, 9, 3 , 0, 2, 8 , 2, 8, 9, 3,
257 5 16,11, 4 , 0, 3,10 , 3,10,11, 4,
258 6 17,13, 1 , 0, 4,12 , 4,12,13, 1,
259 7 14, 6, 7 , 0, 2, 1 , 2, 1, 6, 7,
260 8 15, 8, 9 , 0, 3, 2 , 3, 2, 8, 9,
261 9 16,10,11 , 0, 4, 3 , 4, 3,10,11,
262 . 17,12,13 , 0, 1, 4 , 1, 4,12,13/
263
264! +------+-------------------------------------------+
265! | | sous-triangle voisins |
266! |sous +----------+----------+---------------------+
267! |trian | n3/=n4 | n3=n4 |
268! | +----------+----------+----------+----------+
269! | | IT0 | 2,4,6,8=0| | 2,4,6,8=0|
270! +------+----------+----------+----------+----------+
271! | 1 | 5 2 4 | 5 2 4 | 5 6 8 | 5 6 8 |
272! | 2 | 6 3 1 | 6 3 1 | - - - | - - - |
273! | 3 | 7 4 2 | 7 4 2 | - - - | - - - |
274! | 4 | 8 1 3 | 8 1 3 | - - - | - - - |
275! +------+----------+----------+----------+----------+
276! | 5 | 1 9 13 | 1 -1 -1 | 1 9 13 | 1 -1 -1 |
277! | 6 | 2 10 14 | 2 -1 -1 | 1 10 14 | 1 -1 -1 |
278! | 7 | 3 11 15 | 3 -1 -1 | - - - | - - - |
279! | 8 | 4 12 16 | 4 -1 -1 | 1 12 16 | 1 -1 -1 |
280! +------+----------+----------+----------+----------+
281! | 9 | -1 17 5 | - - - | -1 17 5 | - - - |
282! | 10 | -1 18 6 | - - - | -1 18 6 | - - - |
283! | 11 | -1 19 7 | - - - | - - - | - - - |
284! | 12 | -1 20 8 | - - - | -1 20 8 | - - - |
285! +------+----------+----------+----------+----------+
286! | 13 | -1 5 17 | - - - | -1 5 17 | - - - |
287! | 14 | -1 6 18 | - - - | -1 6 18 | - - - |
288! | 15 | -1 7 19 | - - - | - - - | - - - |
289! | 16 | -1 8 20 | - - - | -1 8 20 | - - - |
290! +------+----------+----------+----------+----------+
291! | 17 | -1 9 13 | - - - | -1 9 13 | - - - |
292! | 18 | -1 10 14 | - - - | -1 10 14 | - - - |
293! | 19 | -1 11 15 | - - - | - - - | - - - |
294! | 20 | -1 12 16 | - - - | -1 12 16 | - - - |
295! +------+----------+----------+----------+----------+
296 DATA it0 /
297 1 5, 2, 4,
298 2 6, 3, 1,
299 3 7, 4, 2,
300 4 8, 1, 3,
301 5 1, 9,13,
302 6 2,10,14,
303 7 3,11,15,
304 8 4,12,16,
305 9 -1,17, 5,
306 . -1,18, 6,
307 1 -1,19, 7,
308 2 -1,20, 8,
309 3 -1, 5,17,
310 4 -1, 6,18,
311 5 -1, 7,19,
312 6 -1, 8,20,
313 7 -1, 9,13,
314 8 -1,10,14,
315 9 -1,11,15,
316 . -1,12,16/
317C-----------------------------------------------
318c if(nft==45952)then
319c ibug=1
320c endif
321 nor_old = 0
322 unp = one + em4
323 zerom = zero - em4
324 eps = (two+half)/hundred
325 unpt = one + eps
326 zeromt = zero - eps
327C------used for high penetration case (extension criterion)
328 eps2=eps*eps
329C------0.1%----
330 eps = em3
331 eps02=eps*eps
332 eps = eps0
333 overw0=one
334 marge025=fourth*marge
335C---
336C--- TOL_INTS,TOL_B : tols for intersection w/ tiny shell, for returning old seg w/ higher pene
337 IF (iresp==1) THEN
338 tol_ints = two*em3
339 tol_b = two*em6
340 ELSE
341 tol_ints = em08
342 tol_b = em08
343 END IF
344 nod1=0
345 nod2=0
346c write(iout,*)'t=',tt,'i24dst3 1'
347c if (ncycle>=0) ip=1
348c if (ncycle<100) write(iout,*)'JLT,ncycle=',JLT,ncycle
349 DO i=1,jlt
350 n=cand_n(i)
351c IRTLM(1,N)>0 node previously impacted on global MAIN IRTLM(1,N)
352c treat one contact but no new intersection
353c IRTLM(1,N)<0 node currently impacted on global MAIN -IRTLM(1,N)
354c treat multiple contact during one cycle
355c IRTLM(1,N)==0 node not impacted
356 ipen0(i)=0
357 icont_r(i) = 0
358 IF(n <= nsn)THEN
359 icontact(i) = irtlm(1,n)
360 subtria(i)=irtlm(2,n)
361 ELSE
362 icontact(i) = irtlm_fi(nin)%P(1,n-nsn)
363 subtria(i) = irtlm_fi(nin)%P(2,n-nsn)
364 ENDIF
365 irtlm_old(i) = icontact(i)
366
367 IF(icontact(i) <= zero)THEN
368 istep(i) = 2
369 icontact(i) = 0
370 ELSEIF(icontact(i) == mseglo(cand_e(i)))THEN
371 istep(i) = 1
372 ELSE
373 icontact(i) = 0
374 istep(i) = 0
375 ENDIF
376 IF(stif(i) == zero)THEN
377 icontact(i) = 0
378 istep(i) = 0
379 ENDIF
380 IF(n <= nsn)THEN
381 IF (pene_old(5,n)> zero.AND.istep(i) == 1) ipen0(i)=1
382 ELSE
383 IF (pene_oldfi(nin)%P(5,n-nsn)> zero.AND.istep(i) == 1) ipen0(i)=1
384 ENDIF
385C-----special Case plan/plan
386 IF(istep(i) == 1.AND.subtria(i) > 20) THEN
387 subtria(i) = subtria(i) - 20
388 istep(i) = 6
389 ENDIF
390 IF(iedge/=0)THEN
391 IF(cand_t(i)==2)THEN
392c CAND_T(I)==2 => only edge candidate
393 icontact(i) = 0
394 istep(i) = 0
395 ENDIF
396 ENDIF
397 nm1(i) = zero
398 nm2(i) = zero
399 nm3(i) = zero
400 ENDDO
401C--------------------------------------------------------
402C ISTEP
403c 0 : rien faire (>STIF=0...)
404c 1 : calcul Hi,pene... si ancien contact
405c si change de segment ou soustrianle st => ISTEP=3
406c sinon (ISTEP=0) plus rien faire
407c 2 : algo d'intersection si pas de contact prcdent
408c si intersection => ISTEP=3
409c sinon plus rien faire
410c 3 : algo de glissement
411c si ne glisse pas => ISTEP=5
412c si glisse mais pas hors surface => ISTEP=4
413c si glisse hors surface => ISTEP=0 plus rien faire
414c 4 : algo de glissement !!!! not used
415c si ne glisse pas hors surface => ISTEP=5
416c si glisse mais pas hors surface => ISTEP=4
417c si glisse hors surface => ISTEP=0 plus rien faire
418c 5 : calcul Hi,pene ... pour nouveau contact ou ancien avec glissement
419c
420c 6 : special shell plane/shell plance case: treated by node/edge(line) with Iedge=1
421C--------------------------------------------------------
422
423C--------------------------------------------------------
424C 0 common initialization node 1 to 5
425C--------------------------------------------------------
426
427 DO i=1,jlt
428 ishel(i)=0
429 IF(istep(i) == 0)cycle
430
431c write(iout,*)'i=',i,'i24dst3 21'
432 xo14 = xx0(i,14)
433 yo14 = yy0(i,14)
434 zo14 = zz0(i,14)
435
436 xo15 = xx0(i,15)
437 yo15 = yy0(i,15)
438 zo15 = zz0(i,15)
439
440 xo16 = xx0(i,16)
441 yo16 = yy0(i,16)
442 zo16 = zz0(i,16)
443
444 xo17 = xx0(i,17)
445 yo17 = yy0(i,17)
446 zo17 = zz0(i,17)
447
448 x51 = xx0(i,1) - xx0(i,5)
449 y51 = yy0(i,1) - yy0(i,5)
450 z51 = zz0(i,1) - zz0(i,5)
451
452 x52 = xx0(i,2) - xx0(i,5)
453 y52 = yy0(i,2) - yy0(i,5)
454 z52 = zz0(i,2) - zz0(i,5)
455
456 x53 = xx0(i,3) - xx0(i,5)
457 y53 = yy0(i,3) - yy0(i,5)
458 z53 = zz0(i,3) - zz0(i,5)
459
460 x54 = xx0(i,4) - xx0(i,5)
461 y54 = yy0(i,4) - yy0(i,5)
462 z54 = zz0(i,4) - zz0(i,5)
463
464
465 sx1250(i) = y51*z52 - z51*y52
466 sy1250(i) = z51*x52 - x51*z52
467 sz1250(i) = x51*y52 - y51*x52
468
469 sx2350(i) = y52*z53 - z52*y53
470 sy2350(i) = z52*x53 - x52*z53
471 sz2350(i) = x52*y53 - y52*x53
472
473 sx3450(i) = y53*z54 - z53*y54
474 sy3450(i) = z53*x54 - x53*z54
475 sz3450(i) = x53*y54 - y53*x54
476
477 sx4150(i) = y54*z51 - z54*y51
478 sy4150(i) = z54*x51 - x54*z51
479 sz4150(i) = x54*y51 - y54*x51
480
481
482 x141 = xx0(i,1) - xx0(i,14)
483 y141 = yy0(i,1) - yy0(i,14)
484 z141 = zz0(i,1) - zz0(i,14)
485
486 x142 = xx0(i,2) - xx0(i,14)
487 y142 = yy0(i,2) - yy0(i,14)
488 z142 = zz0(i,2) - zz0(i,14)
489
490 x152 = xx0(i,2) - xx0(i,15)
491 y152 = yy0(i,2) - yy0(i,15)
492 z152 = zz0(i,2) - zz0(i,15)
493
494 x153 = xx0(i,3) - xx0(i,15)
495 y153 = yy0(i,3) - yy0(i,15)
496 z153 = zz0(i,3) - zz0(i,15)
497
498 x163 = xx0(i,3) - xx0(i,16)
499 y163 = yy0(i,3) - yy0(i,16)
500 z163 = zz0(i,3) - zz0(i,16)
501
502 x164 = xx0(i,4) - xx0(i,16)
503 y164 = yy0(i,4) - yy0(i,16)
504 z164 = zz0(i,4) - zz0(i,16)
505
506 x174 = xx0(i,4) - xx0(i,17)
507 y174 = yy0(i,4) - yy0(i,17)
508 z174 = zz0(i,4) - zz0(i,17)
509
510 x171 = xx0(i,1) - xx0(i,17)
511 y171 = yy0(i,1) - yy0(i,17)
512 z171 = zz0(i,1) - zz0(i,17)
513
514
515 x13_2 = xx0(i,2) - xx0(i,13)
516 y13_2 = yy0(i,2) - yy0(i,13)
517 z13_2 = zz0(i,2) - zz0(i,13)
518
519 x7_3 = xx0(i,3) - xx0(i,7)
520 y7_3 = yy0(i,3) - yy0(i,7)
521 z7_3 = zz0(i,3) - zz0(i,7)
522
523 x9_4 = xx0(i,4) - xx0(i,9)
524 y9_4 = yy0(i,4) - yy0(i,9)
525 z9_4 = zz0(i,4) - zz0(i,9)
526
527 x11_1 = xx0(i,1) - xx0(i,11)
528 y11_1 = yy0(i,1) - yy0(i,11)
529 z11_1 = zz0(i,1) - zz0(i,11)
530
531 x6_4 = xx0(i,4) - xx0(i,6)
532 y6_4 = yy0(i,4) - yy0(i,6)
533 z6_4 = zz0(i,4) - zz0(i,6)
534
535 x8_1 = xx0(i,1) - xx0(i,8)
536 y8_1 = yy0(i,1) - yy0(i,8)
537 z8_1 = zz0(i,1) - zz0(i,8)
538
539 x10_2 = xx0(i,2) - xx0(i,10)
540 y10_2 = yy0(i,2) - yy0(i,10)
541 z10_2 = zz0(i,2) - zz0(i,10)
542
543 x12_3 = xx0(i,3) - xx0(i,12)
544 y12_3 = yy0(i,3) - yy0(i,12)
545 z12_3 = zz0(i,3) - zz0(i,12)
546
547c write(iout,*)'i=',i,'i24dst3 22'
548
549 IF(mvoisin(1,cand_e(i))/=0)THEN
550 sx21140(i) = y142*z141 - z142*y141
551 sy21140(i) = z142*x141 - x142*z141
552 sz21140(i) = x142*y141 - y142*x141
553 ELSE
554 sx21140(i) = sx1250(i)
555 sy21140(i) = sy1250(i)
556 sz21140(i) = sz1250(i)
557 ENDIF
558
559c write(iout,*)'i=',i,'i24dst3 23'
560 IF(ixx(i,3) /= ixx(i,4))THEN
561 IF(mvoisin(2,cand_e(i))/=0)THEN
562 sx32150(i) = y153*z152 - z153*y152
563 sy32150(i) = z153*x152 - x153*z152
564 sz32150(i) = x153*y152 - y153*x152
565 ELSE
566 sx32150(i) = sx2350(i)
567 sy32150(i) = sy2350(i)
568 sz32150(i) = sz2350(i)
569 ENDIF
570
571 IF(mvoisin(3,cand_e(i))/=0)THEN
572 sx43160(i) = y164*z163 - z164*y163
573 sy43160(i) = z164*x163 - x164*z163
574 sz43160(i) = x164*y163 - y164*x163
575 ELSE
576 sx43160(i) = sx3450(i)
577 sy43160(i) = sy3450(i)
578 sz43160(i) = sz3450(i)
579 ENDIF
580
581 IF(mvoisin(4,cand_e(i))/=0)THEN
582 sx14170(i) = y171*z174 - z171*y174
583 sy14170(i) = z171*x174 - x171*z174
584 sz14170(i) = x171*y174 - y171*x174
585 ELSE
586 sx14170(i) = sx4150(i)
587 sy14170(i) = sy4150(i)
588 sz14170(i) = sz4150(i)
589 ENDIF
590
591
592 nxx(i,1) = y13_2 * z6_4 - z13_2 * y6_4
593 nyy(i,1) = z13_2 * x6_4 - x13_2 * z6_4
594 nzz(i,1) = x13_2 * y6_4 - y13_2 * x6_4
595
596 nxx(i,2) = y7_3 * z8_1 - z7_3 * y8_1
597 nyy(i,2) = z7_3 * x8_1 - x7_3 * z8_1
598 nzz(i,2) = x7_3 * y8_1 - y7_3 * x8_1
599
600 nxx(i,3) = y9_4 * z10_2 - z9_4 * y10_2
601 nyy(i,3) = z9_4 * x10_2 - x9_4 * z10_2
602 nzz(i,3) = x9_4 * y10_2 - y9_4 * x10_2
603
604 nxx(i,4) = y11_1 * z12_3 - z11_1 * y12_3
605 nyy(i,4) = z11_1 * x12_3 - x11_1 * z12_3
606 nzz(i,4) = x11_1 * y12_3 - y11_1 * x12_3
607
608 nxx(i,5) = fourth*(nxx(i,1)+nxx(i,2)+nxx(i,3)+nxx(i,4))
609 nyy(i,5) = fourth*(nyy(i,1)+nyy(i,2)+nyy(i,3)+nyy(i,4))
610 nzz(i,5) = fourth*(nzz(i,1)+nzz(i,2)+nzz(i,3)+nzz(i,4))
611
612 ELSE
613
614 IF(mvoisin(2,cand_e(i))/=0)THEN
615 sx32150(i) = y153*z152 - z153*y152
616 sy32150(i) = z153*x152 - x153*z152
617 sz32150(i) = x153*y152 - y153*x152
618 ELSE
619 sx32150(i) = sx1250(i)
620 sy32150(i) = sy1250(i)
621 sz32150(i) = sz1250(i)
622 ENDIF
623
624 IF(mvoisin(3,cand_e(i))/=0)THEN
625 sx43160(i) = y164*z163 - z164*y163
626 sy43160(i) = z164*x163 - x164*z163
627 sz43160(i) = x164*y163 - y164*x163
628 ELSE
629 sx43160(i) = sx1250(i)
630 sy43160(i) = sy1250(i)
631 sz43160(i) = sz1250(i)
632 ENDIF
633
634 IF(mvoisin(4,cand_e(i))/=0)THEN
635 sx14170(i) = y171*z174 - z171*y174
636 sy14170(i) = z171*x174 - x171*z174
637 sz14170(i) = x171*y174 - y171*x174
638 ELSE
639 sx14170(i) = sx1250(i)
640 sy14170(i) = sy1250(i)
641 sz14170(i) = sz1250(i)
642 ENDIF
643
644 nxx(i,1) =sx1250(i)+sx14170(i)+sx21140(i)
645 nyy(i,1) =sy1250(i)+sy14170(i)+sy21140(i)
646 nzz(i,1) =sz1250(i)+sz14170(i)+sz21140(i)
647
648 nxx(i,2) =sx1250(i)+sx21140(i)+sx32150(i)
649 nyy(i,2) =sy1250(i)+sy21140(i)+sy32150(i)
650 nzz(i,2) =sz1250(i)+sz21140(i)+sz32150(i)
651
652 nxx(i,3) =sx1250(i)+sx32150(i)+sx14170(i)
653 nyy(i,3) =sy1250(i)+sy32150(i)+sy14170(i)
654 nzz(i,3) =sz1250(i)+sz32150(i)+sz14170(i)
655
656 nxx(i,4) =nxx(i,3)
657 nyy(i,4) =nyy(i,3)
658 nzz(i,4) =nzz(i,3)
659
660 nxx(i,5) =nxx(i,3)
661 nyy(i,5) =nyy(i,3)
662 nzz(i,5) =nzz(i,3)
663
664 ENDIF
665
666c GAP : coordinates correction
667
668 n=cand_n(i)
669 l=cand_e(i)
670c IF(GAPS(I)+GAP_M(CAND_E(I))
671c . +GAP_NM(5,CAND_E(I))+GAP_NM(6,CAND_E(I))
672c . +GAP_NM(7,CAND_E(I))+GAP_NM(8,CAND_E(I))
673c . +GAP_NM(9,CAND_E(I))+GAP_NM(10,CAND_E(I))
674c . +GAP_NM(11,CAND_E(I))+GAP_NM(12,CAND_E(I)) /= ZERO)THEN
675 IF(gaps(i)>zero.OR.gap_m(l)>zero) THEN
676 ishel(i)=1
677 ELSEIF(gap_nm(5,l)>zero.OR.gap_nm(6,l)>zero.OR.gap_nm(7,l)>zero) THEN
678 ishel(i)=1
679 ELSEIF(gap_nm(8,l)>zero.OR.gap_nm(9,l)>zero.OR.gap_nm(10,l)>zero) THEN
680 ishel(i)=1
681 ELSEIF(gap_nm(11,l)>zero.OR.gap_nm(12,l)>zero) THEN
682 ishel(i)=1
683 END IF
684C-now ISHEL(I)=2 means shell-to-shell since we set GAPS=EM20 instead of 0 when Igap0=1
685 IF(gaps(i)>zero.AND.gap_m(l)>zero) ishel(i)=2
686 IF(ishel(i)>0) THEN
687 nx1 = nxx(i,1)
688 ny1 = nyy(i,1)
689 nz1 = nzz(i,1)
690 aaa = gaps(i)+gap_nm(1,l)
691 aaa = aaa/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
692 nx1 = nx1 * aaa
693 ny1 = ny1 * aaa
694 nz1 = nz1 * aaa
695
696 nx2 = nxx(i,2)
697 ny2 = nyy(i,2)
698 nz2 = nzz(i,2)
699 aaa = gaps(i)+gap_nm(2,l)
700 aaa = aaa/sqrt(nx2*nx2+ny2*ny2+nz2*nz2)
701 nx2 = nx2 * aaa
702 ny2 = ny2 * aaa
703 nz2 = nz2 * aaa
704
705 IF(ixx(i,3) /= ixx(i,4))THEN
706 nx3 = nxx(i,3)
707 ny3 = nyy(i,3)
708 nz3 = nzz(i,3)
709 aaa = gaps(i)+gap_nm(3,l)
710 aaa = aaa/sqrt(nx3*nx3+ny3*ny3+nz3*nz3)
711 nx3 = nx3 * aaa
712 ny3 = ny3 * aaa
713 nz3 = nz3 * aaa
714
715 nx4 = nxx(i,4)
716 ny4 = nyy(i,4)
717 nz4 = nzz(i,4)
718 aaa = gaps(i)+gap_nm(4,l)
719 aaa = aaa/sqrt(nx4*nx4+ny4*ny4+nz4*nz4)
720 nx4 = nx4 * aaa
721 ny4 = ny4 * aaa
722 nz4 = nz4 * aaa
723
724 xx(i,1) = xx0(i,1) + nx1
725 yy(i,1) = yy0(i,1) + ny1
726 zz(i,1) = zz0(i,1) + nz1
727
728 xx(i,2) = xx0(i,2) + nx2
729 yy(i,2) = yy0(i,2) + ny2
730 zz(i,2) = zz0(i,2) + nz2
731
732 xx(i,3) = xx0(i,3) + nx3
733 yy(i,3) = yy0(i,3) + ny3
734 zz(i,3) = zz0(i,3) + nz3
735
736 xx(i,4) = xx0(i,4) + nx4
737 yy(i,4) = yy0(i,4) + ny4
738 zz(i,4) = zz0(i,4) + nz4
739
740 xx(i,5) = fourth*(xx(i,1)+xx(i,2)+xx(i,3)+xx(i,4))
741 yy(i,5) = fourth*(yy(i,1)+yy(i,2)+yy(i,3)+yy(i,4))
742 zz(i,5) = fourth*(zz(i,1)+zz(i,2)+zz(i,3)+zz(i,4))
743 IF(ixx(i,10)==ixx(i,11))THEN
744 xx(i,16) = xx0(i,10)
745 yy(i,16) = yy0(i,10)
746 zz(i,16) = zz0(i,10)
747 ELSE
748 xx(i,16) = fourth*(xx0(i,4)+xx0(i,3)+xx0(i,10)+xx0(i,11))
749 yy(i,16) = fourth*(yy0(i,4)+yy0(i,3)+yy0(i,10)+yy0(i,11))
750 zz(i,16) = fourth*(zz0(i,4)+zz0(i,3)+zz0(i,10)+zz0(i,11))
751 ENDIF
752 nx16 = sx43160(i)
753 ny16 = sy43160(i)
754 nz16 = sz43160(i)
755 aaa = gaps(i)+fourth*(
756 . gap_nm(3,l)+gap_nm(4,l)
757 . + gap_nm(9,l)+gap_nm(10,l) )
758 aaa = aaa/sqrt(nx16*nx16+ny16*ny16+nz16*nz16)
759 xx(i,16) = xx(i,16) + nx16 * aaa
760 yy(i,16) = yy(i,16) + ny16 * aaa
761 zz(i,16) = zz(i,16) + nz16 * aaa
762 ELSE
763 nx3 = nxx(i,3)
764 ny3 = nyy(i,3)
765 nz3 = nzz(i,3)
766 aaa = gaps(i)+gap_nm(3,l)
767 aaa = aaa/sqrt(nx3*nx3+ny3*ny3+nz3*nz3)
768 nx3 = nx3 * aaa
769 ny3 = ny3 * aaa
770 nz3 = nz3 * aaa
771
772 nx4 = nxx(i,4)
773 ny4 = nyy(i,4)
774 nz4 = nzz(i,4)
775
776 xx(i,1) = xx0(i,1) + nx1
777 yy(i,1) = yy0(i,1) + ny1
778 zz(i,1) = zz0(i,1) + nz1
779
780 xx(i,2) = xx0(i,2) + nx2
781 yy(i,2) = yy0(i,2) + ny2
782 zz(i,2) = zz0(i,2) + nz2
783
784 xx(i,3) = xx0(i,3) + nx3
785 yy(i,3) = yy0(i,3) + ny3
786 zz(i,3) = zz0(i,3) + nz3
787
788 xx(i,4) = xx(i,3)
789 yy(i,4) = yy(i,3)
790 zz(i,4) = zz(i,3)
791
792 xx(i,5) = xx(i,3)
793 yy(i,5) = yy(i,3)
794 zz(i,5) = zz(i,3)
795 ENDIF
796 IF(ixx(i,6) == ixx(i,7))THEN
797 xx(i,14) = xx0(i,6)
798 yy(i,14) = yy0(i,6)
799 zz(i,14) = zz0(i,6)
800 ELSE
801 xx(i,14) = fourth*(xx0(i,2)+xx0(i,1)+xx0(i,6)+xx0(i,7))
802 yy(i,14) = fourth*(yy0(i,2)+yy0(i,1)+yy0(i,6)+yy0(i,7))
803 zz(i,14) = fourth*(zz0(i,2)+zz0(i,1)+zz0(i,6)+zz0(i,7))
804 ENDIF
805 nx14 = sx21140(i)
806 ny14 = sy21140(i)
807 nz14 = sz21140(i)
808 aaa = gaps(i)+fourth*(
809 . gap_nm(1,l)+gap_nm(2,l)
810 . + gap_nm(5,l)+gap_nm(6,l) )
811 aaa = aaa/sqrt(nx14*nx14+ny14*ny14+nz14*nz14)
812 xx(i,14) = xx(i,14) + nx14 * aaa
813 yy(i,14) = yy(i,14) + ny14 * aaa
814 zz(i,14) = zz(i,14) + nz14 * aaa
815 IF(ixx(i,8)==ixx(i,9))THEN
816 xx(i,15) = xx0(i,8)
817 yy(i,15) = yy0(i,8)
818 zz(i,15) = zz0(i,8)
819 ELSE
820 xx(i,15) = fourth*(xx0(i,3)+xx0(i,2)+xx0(i,8)+xx0(i,9))
821 yy(i,15) = fourth*(yy0(i,3)+yy0(i,2)+yy0(i,8)+yy0(i,9))
822 zz(i,15) = fourth*(zz0(i,3)+zz0(i,2)+zz0(i,8)+zz0(i,9))
823 ENDIF
824 nx15 = sx32150(i)
825 ny15 = sy32150(i)
826 nz15 = sz32150(i)
827 aaa = gaps(i)+fourth*(
828 . gap_nm(2,l)+gap_nm(3,l)
829 . + gap_nm(7,l)+gap_nm(8,l) )
830 aaa = aaa/sqrt(nx15*nx15+ny15*ny15+nz15*nz15)
831 xx(i,15) = xx(i,15) + nx15 * aaa
832 yy(i,15) = yy(i,15) + ny15 * aaa
833 zz(i,15) = zz(i,15) + nz15 * aaa
834 IF(ixx(i,12)==ixx(i,13))THEN
835 xx(i,17) = xx0(i,12)
836 yy(i,17) = yy0(i,12)
837 zz(i,17) = zz0(i,12)
838 ELSE
839 xx(i,17) = fourth*(xx0(i,1)+xx0(i,4)+xx0(i,12)+xx0(i,13))
840 yy(i,17) = fourth*(yy0(i,1)+yy0(i,4)+yy0(i,12)+yy0(i,13))
841 zz(i,17) = fourth*(zz0(i,1)+zz0(i,4)+zz0(i,12)+zz0(i,13))
842 ENDIF
843 nx17 = sx14170(i)
844 ny17 = sy14170(i)
845 nz17 = sz14170(i)
846 aaa = gaps(i)+fourth*(
847 . gap_nm(4,l)+gap_nm(1,l)
848 . + gap_nm(11,l)+gap_nm(12,l) )
849 aaa = aaa/sqrt(nx17*nx17+ny17*ny17+nz17*nz17)
850 xx(i,17) = xx(i,17) + nx17 * aaa
851 yy(i,17) = yy(i,17) + ny17 * aaa
852 zz(i,17) = zz(i,17) + nz17 * aaa
853
854 x51 = xx(i,1) - xx(i,5)
855 y51 = yy(i,1) - yy(i,5)
856 z51 = zz(i,1) - zz(i,5)
857
858 x52 = xx(i,2) - xx(i,5)
859 y52 = yy(i,2) - yy(i,5)
860 z52 = zz(i,2) - zz(i,5)
861
862 x53 = xx(i,3) - xx(i,5)
863 y53 = yy(i,3) - yy(i,5)
864 z53 = zz(i,3) - zz(i,5)
865
866 x54 = xx(i,4) - xx(i,5)
867 y54 = yy(i,4) - yy(i,5)
868 z54 = zz(i,4) - zz(i,5)
869
870
871 sx125(i) = y51*z52 - z51*y52
872 sy125(i) = z51*x52 - x51*z52
873 sz125(i) = x51*y52 - y51*x52
874
875 sx235(i) = y52*z53 - z52*y53
876 sy235(i) = z52*x53 - x52*z53
877 sz235(i) = x52*y53 - y52*x53
878
879 sx345(i) = y53*z54 - z53*y54
880 sy345(i) = z53*x54 - x53*z54
881 sz345(i) = x53*y54 - y53*x54
882
883 sx415(i) = y54*z51 - z54*y51
884 sy415(i) = z54*x51 - x54*z51
885 sz415(i) = x54*y51 - y54*x51
886
887 ELSE
888
889 xx(i,1:17) = xx0(i,1:17)
890 yy(i,1:17) = yy0(i,1:17)
891 zz(i,1:17) = zz0(i,1:17)
892
893 sx125(i)=sx1250(i)
894 sy125(i)=sy1250(i)
895 sz125(i)=sz1250(i)
896
897 sx235(i) = sx2350(i)
898 sy235(i) = sy2350(i)
899 sz235(i) = sz2350(i)
900
901 sx345(i) = sx3450(i)
902 sy345(i) = sy3450(i)
903 sz345(i) = sz3450(i)
904
905 sx415(i) = sx4150(i)
906 sy415(i) = sy4150(i)
907 sz415(i) = sz4150(i)
908 ENDIF
909c write(iout,*)'i=',i,'i24dst3 212'
910
911 ENDDO
912C--------------------------------------------------------
913C 1 Previous contact compute Hi pene ...
914C--------------------------------------------------------
915c write(iout,*)'jlt=',jlt,'i24dst3 3'
916
917 DO i=1,jlt
918 pene(i) = zero
919 pene_tm(i) = zero
920 IF(istep(i) == 0)cycle
921 IF(istep(i) /= 1.AND.tt /= zero)cycle
922
923C
924 IF (tt==zero) THEN
925 itq = 1
926 ELSE
927 itq = subtria(i)
928 END IF
929 k = ic(1,itq)
930 xa(i) = xx(i,k)
931 ya(i) = yy(i,k)
932 za(i) = zz(i,k)
933
934 k = ic(2,itq)
935 xb(i) = xx(i,k)
936 yb(i) = yy(i,k)
937 zb(i) = zz(i,k)
938
939 k = ic(3,itq)
940 xc(i) = xx(i,k)
941 yc(i) = yy(i,k)
942 zc(i) = zz(i,k)
943
944 k = ic(4,itq)
945 xd(i) = xx(i,k)
946 yd(i) = yy(i,k)
947 zd(i) = zz(i,k)
948
949 k = 15 ! only used for triangles
950 xe(i) = xx(i,k)
951 ye(i) = yy(i,k)
952 ze(i) = zz(i,k)
953
954 k = 17 ! only used for triangles
955 xf(i) = xx(i,k)
956 yf(i) = yy(i,k)
957 zf(i) = zz(i,k)
958
959 IF(itq ==1)THEN
960 nqx = sx125(i)
961 nqy = sy125(i)
962 nqz = sz125(i)
963 ELSEIF(itq ==2)THEN
964 nqx = sx235(i)
965 nqy = sy235(i)
966 nqz = sz235(i)
967 ELSEIF(itq ==3)THEN
968 nqx = sx345(i)
969 nqy = sy345(i)
970 nqz = sz345(i)
971 ELSEIF(itq ==4)THEN
972 nqx = sx415(i)
973 nqy = sy415(i)
974 nqz = sz415(i)
975 ENDIF
976
977 s2 = max(em30,sqrt(nqx*nqx + nqy*nqy + nqz*nqz))
978 area(i) = s2
979 aaa = one/s2
980 nqx = nqx*aaa
981 nqy = nqy*aaa
982 nqz = nqz*aaa
983
984 bbb = (xi(i)-xa(i))*nqx+(yi(i)-ya(i))*nqy+(zi(i)-za(i))*nqz
985C------Add Tiny contact here (replace the starter part)
986 IF (tt == zero .AND. istep(i)/= 1) THEN
987 IF (abs(bbb) < penmin.AND.bbb/=zero) THEN
988C--------------if SECONDARY (projected) is inside segment
989 CALL s_in_m4(xx(i,1),yy(i,1),zz(i,1),xx(i,2),yy(i,2),zz(i,2),
990 1 xx(i,3),yy(i,3),zz(i,3),xx(i,4),yy(i,4),zz(i,4),
991 1 xi(i),yi(i),zi(i),ier)
992 IF (ier ==0) THEN
993 pene(i) = abs(bbb)
994 n=cand_n(i)
995#include "lockon.inc"
996 IF(n <= nsn)THEN
997 IF(irtlm(1,n)==0)THEN
998 irtlm(1,n)=-mseglo(cand_e(i))
999 irtlm(2,n)=1
1000 time_s(n) = pene(i)
1001 ELSE
1002 IF( time_s(n)<pene(i) .OR.
1003 * (time_s(n)==pene(i).AND.-irtlm(1,n) < mseglo(cand_e(i)) ) )THEN
1004 irtlm(1,n)=-mseglo(cand_e(i))
1005 irtlm(2,n)=1
1006 time_s(n) = pene(i)
1007 ENDIF
1008 ENDIF
1009 ELSE
1010 IF(irtlm_fi(nin)%P(1,n-nsn) == 0)THEN
1011 irtlm_fi(nin)%P(1,n-nsn)=-mseglo(cand_e(i))
1012 irtlm_fi(nin)%P(2,n-nsn)=1
1013 time_sfi(nin)%P(n-nsn) = pene(i)
1014 ELSE
1015 IF( time_sfi(nin)%P(n-nsn)<pene(i) .OR.
1016 * (time_sfi(nin)%P(n-nsn)==pene(i).AND.-irtlm_fi(nin)%P(1,n-nsn) < mseglo(cand_e(i)) ) )THEN
1017 irtlm_fi(nin)%P(1,n-nsn)=-mseglo(cand_e(i))
1018 irtlm_fi(nin)%P(2,n-nsn)=1
1019 time_sfi(nin)%P(n-nsn) = pene(i)
1020 ENDIF
1021 ENDIF
1022 ENDIF
1023#include "lockoff.inc"
1024 istep(i)= 1
1025 subtria(i) =1
1026 ELSE
1027 istep(i)= 0
1028 cycle
1029 END IF
1030 ELSE
1031 cycle
1032 END IF
1033 ELSE
1034 pene(i) = -min(zero,bbb)
1035 END IF
1036 largep = 0
1037 IF(pene(i)*pene(i) > eps2*s2.AND.ipen0(i)==0) largep = 1
1038 pene_tm(i) =pene(i)
1039c PMAX_GAP = MAX(PMAX_GAP,ABS(PENE(I)))
1040
1041 n1(i) = nqx
1042 n2(i) = nqy
1043 n3(i) = nqz
1044 n_tm(1,i) = n1(i)
1045 n_tm(2,i) = n2(i)
1046 n_tm(3,i) = n3(i)
1047
1048c project triangle on SECONDARY node
1049
1050
1051 ! a
1052 ! / \
1053 ! / \ a,b,c position at current time
1054 ! / A \ A,B,C projected on SECONDARY node
1055 ! / /\ \
1056 ! / / \ \
1057 ! / / \ \ normal at point a:
1058 ! / / \ \ na = a - a
1059 ! / / \ \
1060 ! / / \ \
1061 ! / / \ \
1062 ! / B--------------C \
1063 ! / \
1064 ! b-------------------------c
1065 !
1066
1067 xab = xb(i)-xa(i)
1068 yab = yb(i)-ya(i)
1069 zab = zb(i)-za(i)
1070 xbc = xc(i)-xb(i)
1071 ybc = yc(i)-yb(i)
1072 zbc = zc(i)-zb(i)
1073 xca = xa(i)-xc(i)
1074 yca = ya(i)-yc(i)
1075 zca = za(i)-zc(i)
1076
1077 xbd = xd(i)-xb(i)
1078 ybd = yd(i)-yb(i)
1079 zbd = zd(i)-zb(i)
1080 xce = xe(i)-xc(i)
1081 yce = ye(i)-yc(i)
1082 zce = ze(i)-zc(i)
1083 xaf = xf(i)-xa(i)
1084 yaf = yf(i)-ya(i)
1085 zaf = zf(i)-za(i)
1086
1087 xia = xa(i)-xi(i)
1088 yia = ya(i)-yi(i)
1089 zia = za(i)-zi(i)
1090 xib = xb(i)-xi(i)
1091 yib = yb(i)-yi(i)
1092 zib = zb(i)-zi(i)
1093 xic = xc(i)-xi(i)
1094 yic = yc(i)-yi(i)
1095 zic = zc(i)-zi(i)
1096 sx = - yab*zca + zab*yca
1097 sy = - zab*xca + xab*zca
1098 sz = - xab*yca + yab*xca
1099 s2 = sx*sx+sy*sy+sz*sz
1100 sax = yib*zic - zib*yic
1101 say = zib*xic - xib*zic
1102 saz = xib*yic - yib*xic
1103 la(i) = (sx*sax+sy*say+sz*saz)/s2
1104 sbx = yic*zia - zic*yia
1105 sby = zic*xia - xic*zia
1106 sbz = xic*yia - yic*xia
1107 lb(i) = (sx*sbx+sy*sby+sz*sbz)/s2
1108 lc(i) = one - la(i) - lb(i)
1109 izlim=0
1110 subtria_tm(i) = subtria(i)
1111 la_tm(i) = la(i)
1112 lb_tm(i) = lb(i)
1113 epseg =eps
1114 epsext = eps
1115C------NSNE: IEDGE4, remove all SECONDARY nodes external extension if Iedge=1
1116 IF (nsne >0) epsext=em04
1117 IF (ispt2(i) >0) epseg=em04
1118 IF(ixx(i,3) /= ixx(i,4))THEN
1119c------------------------------------------------------
1120c quadrangles
1121c check if outside side CA
1122 aaa = zerom
1123 IF(lb(i)<aaa)THEN
1124 istep(i)=3
1125 subtria_n(i)=it0(2,itq)
1126 subtria(i)=it0(2,itq)
1127 xxl(1:17)=xx(i,1:17)
1128 yyl(1:17)=yy(i,1:17)
1129 zzl(1:17)=zz(i,1:17)
1130 izlim=-1
1131 gap2 = gaps(i)+gap_m(cand_e(i))
1132 CALL i24nexttria(
1133 1 izlim ,istep(i),subtria_n(i),subtria(i),
1134 2 largep ,pene(i) ,xxl ,yyl ,zzl ,
1135 3 sx125(i),sy125(i),sz125(i),sx235(i),sy235(i),
1136 4 sz235(i),sx345(i),sy345(i),sz345(i),sx415(i),
1137 5 sy415(i),sz415(i),xi(i) ,yi(i) ,zi(i) ,
1138 6 n1(i) ,n2(i) ,n3(i) ,la(i) ,lb(i) ,
1139 7 lc(i) ,gap2 ,epseg )
1140 subtria_tm(i) = subtria(i)
1141c check if outside side AB
1142 ELSEIF(lc(i)<aaa)THEN
1143 istep(i)=3
1144 subtria_n(i)=it0(3,itq)
1145 subtria(i)=it0(3,itq)
1146 xxl(1:17)=xx(i,1:17)
1147 yyl(1:17)=yy(i,1:17)
1148 zzl(1:17)=zz(i,1:17)
1149 izlim=-1
1150 gap2 = gaps(i)+gap_m(cand_e(i))
1151 CALL i24nexttria(
1152 1 izlim ,istep(i),subtria_n(i),subtria(i),
1153 2 largep ,pene(i) ,xxl ,yyl ,zzl ,
1154 3 sx125(i),sy125(i),sz125(i),sx235(i),sy235(i),
1155 4 sz235(i),sx345(i),sy345(i),sz345(i),sx415(i),
1156 5 sy415(i),sz415(i),xi(i) ,yi(i) ,zi(i) ,
1157 6 n1(i) ,n2(i) ,n3(i) ,la(i) ,lb(i) ,
1158 7 lc(i) ,gap2 ,epseg )
1159 subtria_tm(i) = subtria(i)
1160c check if outside side BC, or when there is rebound and strictly outside
1161 ELSEIF(la(i)<-epsext.OR.(pene(i)==zero.AND.la(i)<zero))THEN
1162c hors zone d'extension de surface
1163 istep(i)=3
1164 subtria_n(i)=it0(1,itq)
1165 izlim=-1
1166 ELSEIF(la(i)<epseg)THEN
1167c zone limite interpolation des normales si convex
1168 vol = n1(i)*xbd + n2(i)*ybd + n3(i)*zbd
1169 gap2 = gaps(i)+gap_m(cand_e(i))
1170 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1171 IF(largep == 1)THEN
1172c large penetration inside +- extension zone
1173 istep(i)=3
1174 subtria_n(i)=it0(1,itq)
1175 izlim = -1
1176 ELSEIF(vol < zero)THEN
1177c convex
1178 izlim = 1
1179 ELSEIF(la(i)<zerom)THEN
1180 istep(i)=3
1181 subtria_n(i)=it0(1,itq)
1182 izlim = -1
1183 ENDIF
1184 ENDIF
1185 IF(izlim == -1)cycle
1186
1187 ELSE
1188c------------------------------------------------------
1189c triangles
1190 IF(la(i)<-epsext.OR.(pene(i)==zero.AND.la(i)<zero))THEN
1191 istep(i)=3
1192 subtria_n(i)=5
1193 izlim = -1
1194 ELSEIF(la(i)<epseg)THEN
1195c zone limite interpolation des normales si convex
1196 vol = n1(i)*xbd + n2(i)*ybd + n3(i)*zbd
1197 gap2 = gaps(i)+gap_m(cand_e(i))
1198 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1199 IF(largep == 1)THEN
1200c large penetration inside +- extension zone
1201 istep(i)=3
1202 subtria_n(i)=5
1203 izlim = -1
1204 ELSEIF(vol < zero)THEN
1205c convex
1206 izlim=1
1207 ELSEIF(la(i)<zerom)THEN
1208 istep(i)=3
1209 subtria_n(i)=5
1210 izlim = -1
1211 ENDIF
1212 ENDIF
1213 IF(lb(i)<-epsext.OR.(pene(i)==zero.AND.lb(i)<zero))THEN
1214 istep(i)=3
1215 IF(lb(i) < la(i)) subtria_n(i)=6
1216 izlim = -1
1217 ELSEIF(lb(i)<epseg)THEN
1218c zone limite interpolation des normales si convex
1219 vol = n1(i)*xce + n2(i)*yce + n3(i)*zce
1220 gap2 = gaps(i)+gap_m(cand_e(i))
1221 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1222 IF(largep == 1)THEN
1223c large penetration inside +- extension zone
1224 istep(i)=3
1225 IF(lb(i) < la(i).OR.izlim==1) subtria_n(i)=6
1226 izlim = -1
1227 ELSEIF(vol < zero)THEN
1228c convex
1229 izlim=1
1230 ELSEIF(lb(i)<zerom)THEN
1231 istep(i)=3
1232 IF(lb(i) < la(i).OR.izlim==1) subtria_n(i)=6
1233 izlim = -1
1234 ENDIF
1235 ENDIF
1236 IF(lc(i)<-epsext.OR.(pene(i)==zero.AND.lc(i)<zero))THEN
1237 istep(i)=3
1238 IF(lc(i) < la(i).AND.lc(i) < lb(i)) subtria_n(i)=8
1239 izlim = -1
1240 ELSEIF(lc(i)<epseg)THEN
1241c zone limite interpolation des normales si convex
1242 vol = n1(i)*xaf + n2(i)*yaf + n3(i)*zaf
1243 gap2 = gaps(i)+gap_m(cand_e(i))
1244 IF (gap2<epseg.AND.(lb(i)<epseg.OR.lc(i)<epseg))vol=-abs(vol)
1245 IF(largep == 1)THEN
1246c large penetration inside +- extension zone
1247 istep(i)=3
1248 IF((lc(i) < la(i).AND.lc(i) < lb(i)).OR.izlim==1) subtria_n(i)=8
1249 izlim = -1
1250 ELSEIF(vol < zero)THEN
1251c convex
1252 izlim=1
1253 ELSEIF(lc(i)<zerom)THEN
1254 istep(i)=3
1255 IF((lc(i) < la(i).AND.lc(i) < lb(i)).OR.izlim==1) subtria_n(i)=8
1256 izlim = -1
1257 ENDIF
1258 ENDIF
1259 IF(izlim == -1)cycle
1260 ENDIF
1261
1262 la_tm(i) = la(i)
1263 lb_tm(i) = lb(i)
1264 IF(la(i)<zero)THEN
1265 IF(lb(i)<zero)THEN
1266 la(i) = zero
1267 lb(i) = zero
1268 lc(i) = one
1269 ELSEIF(lc(i)<zero)THEN
1270 lc(i) = zero
1271 la(i) = zero
1272 lb(i) = one
1273 ELSE
1274 la(i) = zero
1275 aaa = lb(i) + lc(i)
1276 lb(i) = lb(i)/aaa
1277 lc(i) = lc(i)/aaa
1278 ENDIF
1279 ELSEIF(lb(i)<zero)THEN
1280 IF(lc(i)<zero)THEN
1281 lb(i) = zero
1282 lc(i) = zero
1283 la(i) = one
1284 ELSE
1285 lb(i) = zero
1286 aaa = lc(i) + la(i)
1287 lc(i) = lc(i)/aaa
1288 la(i) = la(i)/aaa
1289 ENDIF
1290 ELSEIF(lc(i)<zero)THEN
1291 lc(i) = zero
1292 aaa = la(i) + lb(i)
1293 la(i) = la(i)/aaa
1294 lb(i) = lb(i)/aaa
1295 ENDIF
1296
1297c interpolation des normales si impact dans la zone
1298c +/- extension de surface
1299 IF(izlim==1 .OR. ispt2(i)>0)THEN
1300C--------------- Use for calculate consisting Dpene in i24forc3
1301 nm1(i) = n1(i)
1302 nm2(i) = n2(i)
1303 nm3(i) = n3(i)
1304 k = ic(1,itq)
1305 nax = nxx(i,k)
1306 nay = nyy(i,k)
1307 naz = nzz(i,k)
1308C
1309 k = ic(2,itq)
1310 nbx = nxx(i,k)
1311 nby = nyy(i,k)
1312 nbz = nzz(i,k)
1313C
1314 k = ic(3,itq)
1315 ncx = nxx(i,k)
1316 ncy = nyy(i,k)
1317 ncz = nzz(i,k)
1318 n1(i) = la(i)*nax + lb(i)*nbx + lc(i)*ncx
1319 n2(i) = la(i)*nay + lb(i)*nby + lc(i)*ncy
1320 n3(i) = la(i)*naz + lb(i)*nbz + lc(i)*ncz
1321
1322
1323 nni = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
1324 ni2 = n1(i)*n1(i) + n2(i)*n2(i) + n3(i)*n3(i)
1325 IF(two*nni*nni < ni2)THEN
1326c scharp angle bound nodal normal to 45 from segment normal
1327 aaa = sqrt(ni2-nni*nni) - nni
1328 n1(i) = n1(i) + aaa*nqx
1329 n2(i) = n2(i) + aaa*nqy
1330 n3(i) = n3(i) + aaa*nqz
1331 ni2 = n1(i)*n1(i) + n2(i)*n2(i) + n3(i)*n3(i)
1332 ENDIF
1333 s2 = one/max(em30,sqrt(ni2))
1334 n1(i) = n1(i)*s2
1335 n2(i) = n2(i)*s2
1336 n3(i) = n3(i)*s2
1337C-------special case, skew mesh------------
1338 nni = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
1339 IF (nni<em20) THEN
1340 n1(i) = nqx
1341 n2(i) = nqy
1342 n3(i) = nqz
1343 END IF
1344C NNI = NQX*N1(I) + NQY*N2(I) + NQZ*N3(I)
1345C AAA = MAX(ONE,ONE/MAX(NNI,EM20))
1346C PENE(I) = PENE(I)*AAA
1347
1348 ENDIF
1349
1350 IF(pene(i) == zero)THEN
1351 icontact(i) = 0
1352 istep(i) = 0
1353 cycle
1354 ENDIF
1355
1356 ix1(i) = ixx(i,ic( 7,itq))
1357 ix2(i) = ixx(i,ic( 8,itq))
1358 ix3(i) = ixx(i,ic( 9,itq))
1359 ix4(i) = ixx(i,ic(10,itq))
1360 x1(i) = xx0(i,ic( 7,itq))
1361 x2(i) = xx0(i,ic( 8,itq))
1362 x3(i) = xx0(i,ic( 9,itq))
1363 x4(i) = xx0(i,ic(10,itq))
1364 y1(i) = yy0(i,ic( 7,itq))
1365 y2(i) = yy0(i,ic( 8,itq))
1366 y3(i) = yy0(i,ic( 9,itq))
1367 y4(i) = yy0(i,ic(10,itq))
1368 z1(i) = zz0(i,ic( 7,itq))
1369 z2(i) = zz0(i,ic( 8,itq))
1370 z3(i) = zz0(i,ic( 9,itq))
1371 z4(i) = zz0(i,ic(10,itq))
1372
1373 vx1(i) = vx(i,ic( 7,itq))
1374 vx2(i) = vx(i,ic( 8,itq))
1375 vx3(i) = vx(i,ic( 9,itq))
1376 vx4(i) = vx(i,ic(10,itq))
1377 vy1(i) = vy(i,ic( 7,itq))
1378 vy2(i) = vy(i,ic( 8,itq))
1379 vy3(i) = vy(i,ic( 9,itq))
1380 vy4(i) = vy(i,ic(10,itq))
1381 vz1(i) = vz(i,ic( 7,itq))
1382 vz2(i) = vz(i,ic( 8,itq))
1383 vz3(i) = vz(i,ic( 9,itq))
1384 vz4(i) = vz(i,ic(10,itq))
1385
1386 IF (ix1(i) == ix2(i))THEN
1387 h1(i) = la(i)
1388 h2(i) = zero
1389 h3(i) = lb(i)
1390 h4(i) = lc(i)
1391 ELSEIF(ix2(i) == ix3(i))THEN
1392 h1(i) = la(i)
1393 h2(i) = lb(i)
1394 h3(i) = zero
1395 h4(i) = lc(i)
1396c ELSEIF(IX3(I) == IX4(I))THEN impossible
1397 ELSEIF(ix4(i) == ix1(i))THEN
1398 h1(i) = lc(i)
1399 h2(i) = la(i)
1400 h3(i) = lb(i)
1401 h4(i) = zero
1402 ELSE
1403 h0 = fourth * la(i)
1404 h1(i) = h0
1405 h2(i) = h0
1406 h3(i) = h0 + lb(i)
1407 h4(i) = h0 + lc(i)
1408 ENDIF
1409 ENDDO
1410C--------------------------------------------------------
1411C 2 Search intersection
1412C--------------------------------------------------------
1413c write(iout,*)'i24dst3 5'
1414 DO i=1,jlt
1415 IF(istep(i) /= 2)cycle ! previous contact
1416c
1417
1418 xi5 = xx(i,5) - xi(i)
1419 yi5 = yy(i,5) - yi(i)
1420 zi5 = zz(i,5) - zi(i)
1421
1422 v12 = sx125(i)*xi5 + sy125(i)*yi5 + sz125(i)*zi5
1423 v23 = sx235(i)*xi5 + sy235(i)*yi5 + sz235(i)*zi5
1424 v34 = sx345(i)*xi5 + sy345(i)*yi5 + sz345(i)*zi5
1425 v41 = sx415(i)*xi5 + sy415(i)*yi5 + sz415(i)*zi5
1426
1427 IF(v12 < zero .and. v23 < zero .and.
1428 . v34 < zero .and. v41 < zero ) THEN
1429c INTERSECTION IMPOSSIBLE
1430 cycle
1431 ENDIF
1432
1433 n=cand_n(i)
1434
1435c POSSIBLE INTERSECTION
1436
1437C--------------------------------------------------------
1438C OLD COORDINATES
1439C--------------------------------------------------------
1440 IF(intply > 0) THEN
1441C
1442 dgap = dgap_nm(1,cand_e(i))
1443 nx1 = nxx(i,1)
1444 ny1 = nyy(i,1)
1445 nz1 = nzz(i,1)
1446 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1447 nx1 = nx1 * aaa
1448 ny1 = ny1 * aaa
1449 nz1 = nz1 * aaa
1450C
1451 xo1 = xx(i,1) - dt1*vx(i,1) - dgap*nx1
1452 yo1 = yy(i,1) - dt1*vy(i,1) - dgap*ny1
1453 zo1 = zz(i,1) - dt1*vz(i,1) - dgap*nz1
1454 dgap = dgap_nm(2,cand_e(i))
1455 nx1 = nxx(i,2)
1456 ny1 = nyy(i,2)
1457 nz1 = nzz(i,2)
1458 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1459 nx1 = nx1 * aaa
1460 ny1 = ny1 * aaa
1461 nz1 = nz1 * aaa
1462C
1463 xo2 = xx(i,2) - dt1*vx(i,2)- dgap*nx1
1464 yo2 = yy(i,2) - dt1*vy(i,2)- dgap*ny1
1465 zo2 = zz(i,2) - dt1*vz(i,2)- dgap*nz1
1466 dgap = dgap_nm(3,cand_e(i))
1467 nx1 = nxx(i,3)
1468 ny1 = nyy(i,3)
1469 nz1 = nzz(i,3)
1470 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1471 nx1 = nx1 * aaa
1472 ny1 = ny1 * aaa
1473 nz1 = nz1 * aaa
1474C
1475 xo3 = xx(i,3) - dt1*vx(i,3)- dgap*nx1
1476 yo3 = yy(i,3) - dt1*vy(i,3)- dgap*ny1
1477 zo3 = zz(i,3) - dt1*vz(i,3)- dgap*nz1
1478 dgap = dgap_nm(4,cand_e(i))
1479 nx1 = nxx(i,4)
1480 ny1 = nyy(i,4)
1481 nz1 = nzz(i,4)
1482 aaa = one/sqrt(nx1*nx1+ny1*ny1+nz1*nz1)
1483 nx1 = nx1 * aaa
1484 ny1 = ny1 * aaa
1485 nz1 = nz1 * aaa
1486C
1487 xo4 = xx(i,4) - dt1*vx(i,4)- dgap*nx1
1488 yo4 = yy(i,4) - dt1*vy(i,4)- dgap*ny1
1489 zo4 = zz(i,4) - dt1*vz(i,4)- dgap*nz1
1490 ELSE
1491 xo1 = xx(i,1) - dt1*vx(i,1)
1492 yo1 = yy(i,1) - dt1*vy(i,1)
1493 zo1 = zz(i,1) - dt1*vz(i,1)
1494C
1495 xo2 = xx(i,2) - dt1*vx(i,2)
1496 yo2 = yy(i,2) - dt1*vy(i,2)
1497 zo2 = zz(i,2) - dt1*vz(i,2)
1498c
1499 xo3 = xx(i,3) - dt1*vx(i,3)
1500 yo3 = yy(i,3) - dt1*vy(i,3)
1501 zo3 = zz(i,3) - dt1*vz(i,3)
1502C
1503 xo4 = xx(i,4) - dt1*vx(i,4)
1504 yo4 = yy(i,4) - dt1*vy(i,4)
1505 zo4 = zz(i,4) - dt1*vz(i,4)
1506 ENDIF
1507
1508 xoi = xi(i) - dt1*vxi(i)
1509 yoi = yi(i) - dt1*vyi(i)
1510 zoi = zi(i) - dt1*vzi(i)
1511
1512 IF(ixx(i,3) /= ixx(i,4))THEN
1513 xo5 = fourth*(xo1+xo2+xo3+xo4)
1514 yo5 = fourth*(yo1+yo2+yo3+yo4)
1515 zo5 = fourth*(zo1+zo2+zo3+zo4)
1516 ELSE
1517 xo5 = xo3
1518 yo5 = yo3
1519 zo5 = zo3
1520 ENDIF
1521
1522c compute volumes at previous time step
1523
1524 x51 = xo1 - xo5
1525 y51 = yo1 - yo5
1526 z51 = zo1 - zo5
1527
1528 x52 = xo2 - xo5
1529 y52 = yo2 - yo5
1530 z52 = zo2 - zo5
1531
1532 x53 = xo3 - xo5
1533 y53 = yo3 - yo5
1534 z53 = zo3 - zo5
1535
1536 x54 = xo4 - xo5
1537 y54 = yo4 - yo5
1538 z54 = zo4 - zo5
1539
1540 xi5 = xo5 - xoi
1541 yi5 = yo5 - yoi
1542 zi5 = zo5 - zoi
1543
1544 sox125 = y51*z52 - z51*y52
1545 soy125 = z51*x52 - x51*z52
1546 soz125 = x51*y52 - y51*x52
1547 vo12 = sox125*xi5 + soy125*yi5 + soz125*zi5
1548
1549 sox235 = y52*z53 - z52*y53
1550 soy235 = z52*x53 - x52*z53
1551 soz235 = x52*y53 - y52*x53
1552 vo23 = sox235*xi5 + soy235*yi5 + soz235*zi5
1553
1554 sox345 = y53*z54 - z53*y54
1555 soy345 = z53*x54 - x53*z54
1556 soz345 = x53*y54 - y53*x54
1557 vo34 = sox345*xi5 + soy345*yi5 + soz345*zi5
1558
1559 sox415 = y54*z51 - z54*y51
1560 soy415 = z54*x51 - x54*z51
1561 soz415 = x54*y51 - y54*x51
1562 vo41 = sox415*xi5 + soy415*yi5 + soz415*zi5
1563
1564 adt0(i) = -one
1565
1566c compute intersection time (linear approximation)
1567c compute coordinates at intersection time
1568c intersection time can be different for each sub-triangle
1569
1570 subtria(i)= 0 ! subtriangle 0 => no contact
1571
1572 intersect = 0
1573 CALL intersecv(
1574 1 istep(i) ,subtria(i),ixx(i,3) ,ixx(i,4),
1575 1 intersect,v12 ,v23 ,v34 ,v41 ,
1576 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
1577 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2) ,
1578 4 xx(i,3 ) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4) ,
1579 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),adt0(i) ,
1580 6 vxi(i) ,vyi(i) ,vzi(i) ,xo1 ,xo2 ,
1581 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1582 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1583 9 zo3 ,zo4 ,zo5 ,nx(i) ,ny(i) ,
1584 a nz(i) ,sx125(i),sx235(i),sx345(i),sx415(i),
1585 b sy125(i) ,sy235(i),sy345(i),sy415(i),sz125(i),
1586 c sz235(i) ,sz345(i),sz415(i),xi(i) ,yi(i) ,
1587 d zi(i) ,zerom ,unp ,zeromt,unpt )
1588 IF(intersect == 0.AND.ishel(i)>1.AND.ispt2(i)>0)THEN
1589 xi5 = xx0(i,5) - xi(i)
1590 yi5 = yy0(i,5) - yi(i)
1591 zi5 = zz0(i,5) - zi(i)
1592 aaa = abs(sx1250(i)*xi5 + sy1250(i)*yi5 + sz1250(i)*zi5)
1593 aaa = max(aaa,(sx2350(i)*xi5 + sy2350(i)*yi5 + sz2350(i)*zi5))
1594 aaa = max(aaa,(sx3450(i)*xi5 + sy3450(i)*yi5 + sz3450(i)*zi5))
1595 aaa = max(aaa,(sx4150(i)*xi5 + sy4150(i)*yi5 + sz4150(i)*zi5))
1596 IF (aaa<five*em04)
1597 * CALL intersecsh(
1598 1 istep(i) ,subtria(i),ixx(i,3) ,ixx(i,4),
1599 1 intersect,xx(i,1),
1600 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2) ,
1601 4 xx(i,3 ) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4) ,
1602 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),adt0(i) ,
1603 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
1604 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1605 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1606 9 zo3 ,zo4 ,zo5 ,nx(i) ,ny(i) ,
1607 a nz(i) ,sx125(i),sx235(i),sx345(i),sx415(i),
1608 b sy125(i) ,sy235(i),sy345(i),sy415(i),sz125(i),
1609 c sz235(i) ,sz345(i),sz415(i),xi(i) ,yi(i) ,
1610 d zi(i) ,sox125 ,sox235 ,sox345 ,sox415 ,
1611 b soy125 ,soy235 ,soy345 ,soy415 ,soz125 ,
1612 c soz235 ,soz345 ,soz415 ,mvoisin(1,cand_e(i)) )
1613 END IF
1614 ikeep = 1
1615 IF(intersect == 0.AND.n <= nsn)THEN
1616 IF (icont_i(n) < 0 ) ikeep = 0
1617 ELSEIF(intersect == 0) THEN
1618 IF (icont_i_fi(nin)%P(n-nsn) < 0 ) ikeep = 0
1619 ENDIF
1620C----- Inacti=5 all 1er contact has 0 force therefore no risk
1621 IF(ikeep == 1.AND.ishel(i)>0.AND.ispt2(i)==0.AND.inacti/=5)THEN
1622 IF ((gaps(i)+gap_m(cand_e(i)))>tol_ints) ikeep = 0
1623 END IF
1624 IF(ikeep == 1.AND.intersect == 0)THEN
1625 CALL intersecv0(
1626 1 istep(i) ,subtria(i),ixx(i,3) ,ixx(i,4),
1627 1 intersect,v12 ,v23 ,v34 ,v41 ,
1628 2 vo12 ,vo23 ,vo34 ,vo41 ,xx(i,1),
1629 3 yy(i,1) ,zz(i,1),xx(i,2),yy(i,2),zz(i,2) ,
1630 4 xx(i,3 ) ,yy(i,3),zz(i,3),xx(i,4),yy(i,4) ,
1631 5 zz(i,4) ,xx(i,5),yy(i,5),zz(i,5),adt0(i) ,
1632 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
1633 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
1634 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
1635 9 zo3 ,zo4 ,zo5 ,nx(i) ,ny(i) ,
1636 a nz(i) ,sx125(i),sx235(i),sx345(i),sx415(i),
1637 b sy125(i) ,sy235(i),sy345(i),sy415(i),sz125(i),
1638 c sz235(i) ,sz345(i),sz415(i),xi(i) ,yi(i) ,
1639 d zi(i) ,zerom ,unp )
1640C----- due to different algorithms used in starter/engine, it's more assurant to do the correction here
1641 IF (tt==zero.AND.intersect>0.AND.inacti==0) THEN
1642 IF(n <= nsn)THEN
1643 icont_i(n) = -cand_e(i)
1644 ELSE
1645 icont_i_fi(nin)%P(n-nsn) = -cand_e(i)
1646 ENDIF
1647 intersect = 0
1648 istep(i) = 2
1649 END IF
1650 icont_r(i) = intersect
1651 END if!(IKEEP == 1)THEN
1652C
1653 ENDDO
1654! numerotation des sous-triangles
1655!
1656!
1657! 6o------o5
1658! |\ 19 /|
1659! | \ / |
1660! | \/ | 7,8=0o------3=4
1661! |15/\11| | /|\
1662! | / \ | | 8 / | \
1663! |/ 7 \| | / | \
1664! 7o------4------3------o4 | / | \
1665! |\ 12 /|\ /|\ 14 /| | / 1 | 6 \
1666! | \ / | \3 / | \ / | |/ | \
1667! | \/ | \/2 |6 \/18| 1------2------o3,4=0
1668! |20/\ 8| 4/\ | /\ | | /
1669! | / \ | / 1\ | / \ | | 5 /
1670! |/ 16 \|/ \|/ 10 \| | /
1671! 8o------1------2------o3 | /
1672! |\ 5 /| | /
1673! | \ / | |/
1674! |9 \/13| 1o2=0
1675! | /\ |
1676! | / \ |
1677! |/ 17 \|
1678! 1o------o2
1679!
1680! +------+-------------------------------------------+
1681! | | sous-triangle voisins |
1682! |sous +----------+----------+---------------------+
1683! |trian | n3/=n4 | n3=n4 |
1684! | +----------+----------+----------+----------+
1685! | | IT0 | 2,4,6,8=0| | 2,4,6,8=0|
1686! +------+----------+----------+----------+----------+
1687! | 1 | 5 2 4 | 5 2 4 | 5 6 8 | 5 6 8 |
1688! | 2 | 6 3 1 | 6 3 1 | - - - | - - - |
1689! | 3 | 7 4 2 | 7 4 2 | - - - | - - - |
1690! | 4 | 8 1 3 | 8 1 3 | - - - | - - - |
1691! +------+----------+----------+----------+----------+
1692! | 5 | 1 9 13 | 1 -1 -1 | 1 9 13 | 1 -1 -1 |
1693! | 6 | 2 10 14 | 2 -1 -1 | 1 10 14 | 1 -1 -1 |
1694! | 7 | 3 11 15 | 3 -1 -1 | - - - | - - - |
1695! | 8 | 4 12 16 | 4 -1 -1 | 1 12 16 | 1 -1 -1 |
1696! +------+----------+----------+----------+----------+
1697! | 9 | 5 -1 17 | - - - | 5 -1 17 | - - - |
1698! | 10 | 6 -1 18 | - - - | 6 -1 18 | - - - |
1699! | 11 | 7 -1 19 | - - - | - - - | - - - |
1700! | 12 | 8 -1 20 | - - - | 8 -1 20 | - - - |
1701! +------+----------+----------+----------+----------+
1702! | 13 | 5 17 -1 | - - - | 5 17 -1 | - - - |
1703! | 14 | 6 18 -1 | - - - | 6 18 -1 | - - - |
1704! | 15 | 7 19 -1 | - - - | - - - | - - - |
1705! | 16 | 8 20 -1 | - - - | 8 20 -1 | - - - |
1706! +------+----------+----------+----------+----------+
1707! | 17 | 9 13 -1 | - - - | 9 13 -1 | - - - |
1708! | 18 | 10 14 -1 | - - - | 10 14 -1 | - - - |
1709! | 19 | 11 15 -1 | - - - | - - - | - - - |
1710! | 20 | 12 16 -1 | - - - | 12 16 -1 | - - - |
1711! +------+----------+----------+----------+----------+
1712!
1713!
1714! connetivite des sous-triangles IC(10,*)
1715!
1716!
1717! +----+----------+----------+-------------+
1718! | ST | noeuds T | noeuds TV| noeuds Quad |
1719! +----+----------+----------+-------------+
1720! | 1 | 5 1 2 | 14 3 4 | 3 4 1 2 |
1721! | 2 | 5 2 3 | 15 4 1 | 4 1 2 3 | 11-------10
1722! | 3 | 5 3 4 | 16 1 2 | 1 2 3 4 | |\ 19 /|
1723! | 4 | 5 4 1 | 17 2 3 | 2 3 4 1 | | \ / |
1724! +----+----------+----------+-------------+ | \ / |
1725! | 5 | 14 2 1 | 5 6 7 | 6 7 2 1 | | 16 |
1726! | 6 | 15 3 2 | 5 8 9 | 8 9 3 2 | |15/ \11|
1727! | 7 | 16 4 3 | 5 10 11 | 10 11 4 3 | | / \ |
1728! | 8 | 17 1 4 | 5 12 13 | 12 13 1 4 | |/ 7 \|
1729! +----+----------+----------+-------------+12-------4-------3-------9
1730! | 9 | 14 1 6 |(13) 7 2 | 7 2 1 6 | |\ 12 /|\ /|\ 14 /|
1731! | 10 | 15 2 8 |( 7) 9 3 | 9 3 2 8 | | \ / | \ 3 / | \ / |
1732! | 11 | 16 3 10 |( 9)11 4 | 11 4 3 10 | | \ / | \ /2 |6 \ /18|
1733! | 12 | 17 4 12 |(11)13 1 | 13 1 4 12 | | 17 | 5 | 15 |
1734! +----+----------+----------+-------------+ |20/ \ 8| 4/ \ | / \ |
1735! | 13 | 14 7 2 |( 8) 1 6 | 1 6 7 2 | | / \ | / 1 \ | / \ |
1736! | 14 | 15 9 3 |(10) 2 8 | 2 8 9 3 | |/ 16 \|/ \|/ 10 \|
1737! | 15 | 16 11 4 |(12) 3 10 | 3 10 11 4 |13-------1-------2-------8
1738! | 16 | 17 13 1 |( 6) 4 12 | 4 12 13 1 | |\ 5 /|
1739! +----+----------+----------+-------------+ | \ / |
1740! | 17 | 14 6 7 | - 2 1 | 2 1 6 7 | |9 \ /13|
1741! | 18 | 15 8 9 | - 3 2 | 3 2 8 9 | | 14 |
1742! | 19 | 16 10 11 | - 4 3 | 4 3 10 11 | | / \ |
1743! | 20 | 17 12 13 | - 1 4 | 1 4 12 13 | | / \ |
1744! +----+----------+----------+-------------+ |/ 17 \|
1745! 6-------7
1746!
1747C=======================================================================
1748c sliding
1749C=======================================================================
1750
1751C--------------------------------------------------------
1752C 0 common initialization node 6 to 17
1753C--------------------------------------------------------
1754
1755 DO i=1,jlt
1756
1757 IF(istep(i) /= 3)cycle
1758
1759
1760 nxx(i,6) =sx21140(i)
1761 nyy(i,6) =sy21140(i)
1762 nzz(i,6) =sz21140(i)
1763
1764 nxx(i,7) =sx21140(i)
1765 nyy(i,7) =sy21140(i)
1766 nzz(i,7) =sz21140(i)
1767
1768 nxx(i,8) =sx32150(i)
1769 nyy(i,8) =sy32150(i)
1770 nzz(i,8) =sz32150(i)
1771
1772 nxx(i,9) =sx32150(i)
1773 nyy(i,9) =sy32150(i)
1774 nzz(i,9) =sz32150(i)
1775
1776 nxx(i,15)=sx32150(i)
1777 nyy(i,15)=sy32150(i)
1778 nzz(i,15)=sz32150(i)
1779
1780 nxx(i,10)=sx43160(i)
1781 nyy(i,10)=sy43160(i)
1782 nzz(i,10)=sz43160(i)
1783
1784 nxx(i,11)=sx43160(i)
1785 nyy(i,11)=sy43160(i)
1786 nzz(i,11)=sz43160(i)
1787
1788 nxx(i,14)=sx21140(i)
1789 nyy(i,14)=sy21140(i)
1790 nzz(i,14)=sz21140(i)
1791
1792 nxx(i,16)=sx43160(i)
1793 nyy(i,16)=sy43160(i)
1794 nzz(i,16)=sz43160(i)
1795
1796 nxx(i,12)=sx14170(i)
1797 nyy(i,12)=sy14170(i)
1798 nzz(i,12)=sz14170(i)
1799
1800 nxx(i,13)=sx14170(i)
1801 nyy(i,13)=sy14170(i)
1802 nzz(i,13)=sz14170(i)
1803
1804 nxx(i,17)=sx14170(i)
1805 nyy(i,17)=sy14170(i)
1806 nzz(i,17)=sz14170(i)
1807
1808c GAP : coordinates correction
1809
1810 n=cand_n(i)
1811
1812c IF(GAPS(I)+GAP_M(CAND_E(I))
1813c . +GAP_NM(5,CAND_E(I))+GAP_NM(6,CAND_E(I))
1814c . +GAP_NM(7,CAND_E(I))+GAP_NM(8,CAND_E(I))
1815c . +GAP_NM(9,CAND_E(I))+GAP_NM(10,CAND_E(I))
1816c . +GAP_NM(11,CAND_E(I))+GAP_NM(12,CAND_E(I)) /= ZERO)THEN
1817 IF(ishel(i)>0) THEN
1818
1819c la normale sur les point 6 a 17 est approche
1820c pour la correction du gap (utilise pendant un seul cycle)
1821
1822 bbb = sqrt(sx21140(i)*sx21140(i)
1823 + +sy21140(i)*sy21140(i)
1824 + +sz21140(i)*sz21140(i))
1825 aaa = gaps(i)+gap_nm(5,cand_e(i))
1826 aaa = aaa/bbb
1827 xx(i,6) = xx0(i,6) + sx21140(i) * aaa
1828 yy(i,6) = yy0(i,6) + sy21140(i) * aaa
1829 zz(i,6) = zz0(i,6) + sz21140(i) * aaa
1830 aaa = gaps(i)+gap_nm(6,cand_e(i))
1831 aaa = aaa/bbb
1832 xx(i,7) = xx0(i,7) + sx21140(i) * aaa
1833 yy(i,7) = yy0(i,7) + sy21140(i) * aaa
1834 zz(i,7) = zz0(i,7) + sz21140(i) * aaa
1835
1836 bbb = sqrt(sx32150(i)*sx32150(i)
1837 + +sy32150(i)*sy32150(i)
1838 + +sz32150(i)*sz32150(i))
1839 aaa = gaps(i)+gap_nm(7,cand_e(i))
1840 aaa = aaa/bbb
1841 xx(i,8) = xx0(i,8) + sx32150(i) * aaa
1842 yy(i,8) = yy0(i,8) + sy32150(i) * aaa
1843 zz(i,8) = zz0(i,8) + sz32150(i) * aaa
1844 aaa = gaps(i)+gap_nm(8,cand_e(i))
1845 aaa = aaa/bbb
1846 xx(i,9) = xx0(i,9) + sx32150(i) * aaa
1847 yy(i,9) = yy0(i,9) + sy32150(i) * aaa
1848 zz(i,9) = zz0(i,9) + sz32150(i) * aaa
1849
1850 IF(ixx(i,3) == ixx(i,4))THEN
1851
1852 xx(i,10) = xx0(i,10)
1853 yy(i,10) = yy0(i,10)
1854 zz(i,10) = zz0(i,10)
1855
1856 xx(i,11) = xx0(i,11)
1857 yy(i,11) = yy0(i,11)
1858 zz(i,11) = zz0(i,11)
1859
1860 ELSE
1861
1862 bbb = sqrt(sx43160(i)*sx43160(i)
1863 + +sy43160(i)*sy43160(i)
1864 + +sz43160(i)*sz43160(i))
1865 aaa = gaps(i)+gap_nm(9,cand_e(i))
1866 aaa = aaa/bbb
1867 xx(i,10) = xx0(i,10) + sx43160(i) * aaa
1868 yy(i,10) = yy0(i,10) + sy43160(i) * aaa
1869 zz(i,10) = zz0(i,10) + sz43160(i) * aaa
1870 aaa = gaps(i)+gap_nm(10,cand_e(i))
1871 aaa = aaa/bbb
1872 xx(i,11) = xx0(i,11) + sx43160(i) * aaa
1873 yy(i,11) = yy0(i,11) + sy43160(i) * aaa
1874 zz(i,11) = zz0(i,11) + sz43160(i) * aaa
1875 ENDIF
1876
1877 bbb = sqrt(sx14170(i)*sx14170(i)
1878 + +sy14170(i)*sy14170(i)
1879 + +sz14170(i)*sz14170(i))
1880 aaa = gaps(i)+gap_nm(11,cand_e(i))
1881 aaa = aaa/bbb
1882 xx(i,12) = xx0(i,12) + sx14170(i) * aaa
1883 yy(i,12) = yy0(i,12) + sy14170(i) * aaa
1884 zz(i,12) = zz0(i,12) + sz14170(i) * aaa
1885 aaa = gaps(i)+gap_nm(12,cand_e(i))
1886 aaa = aaa/bbb
1887 xx(i,13) = xx0(i,13) + sx14170(i) * aaa
1888 yy(i,13) = yy0(i,13) + sy14170(i) * aaa
1889 zz(i,13) = zz0(i,13) + sz14170(i) * aaa
1890
1891c write(iout,*)'i=',i,'i24dst3 29'
1892 IF(ixx(i,6) == ixx(i,7))THEN
1893 xx(i,14) = xx(i,6)
1894 yy(i,14) = yy(i,6)
1895 zz(i,14) = zz(i,6)
1896 ELSE
1897 xx(i,14) = fourth*(xx(i,2)+xx(i,1)+xx(i,6)+xx(i,7))
1898 yy(i,14) = fourth*(yy(i,2)+yy(i,1)+yy(i,6)+yy(i,7))
1899 zz(i,14) = fourth*(zz(i,2)+zz(i,1)+zz(i,6)+zz(i,7))
1900 ENDIF
1901 IF(ixx(i,8)==ixx(i,9))THEN
1902 xx(i,15) = xx(i,8)
1903 yy(i,15) = yy(i,8)
1904 zz(i,15) = zz(i,8)
1905 ELSE
1906 xx(i,15) = fourth*(xx(i,3)+xx(i,2)+xx(i,8)+xx(i,9))
1907 yy(i,15) = fourth*(yy(i,3)+yy(i,2)+yy(i,8)+yy(i,9))
1908 zz(i,15) = fourth*(zz(i,3)+zz(i,2)+zz(i,8)+zz(i,9))
1909 ENDIF
1910 IF(ixx(i,10)==ixx(i,11))THEN
1911 xx(i,16) = xx(i,10)
1912 yy(i,16) = yy(i,10)
1913 zz(i,16) = zz(i,10)
1914 ELSE
1915 xx(i,16) = fourth*(xx(i,4)+xx(i,3)+xx(i,10)+xx(i,11))
1916 yy(i,16) = fourth*(yy(i,4)+yy(i,3)+yy(i,10)+yy(i,11))
1917 zz(i,16) = fourth*(zz(i,4)+zz(i,3)+zz(i,10)+zz(i,11))
1918 ENDIF
1919 IF(ixx(i,12)==ixx(i,13))THEN
1920 xx(i,17) = xx(i,12)
1921 yy(i,17) = yy(i,12)
1922 zz(i,17) = zz(i,12)
1923 ELSE
1924 xx(i,17) = fourth*(xx(i,1)+xx(i,4)+xx(i,12)+xx(i,13))
1925 yy(i,17) = fourth*(yy(i,1)+yy(i,4)+yy(i,12)+yy(i,13))
1926 zz(i,17) = fourth*(zz(i,1)+zz(i,4)+zz(i,12)+zz(i,13))
1927 ENDIF
1928
1929 x141 = xx(i,1) - xx(i,14)
1930 y141 = yy(i,1) - yy(i,14)
1931 z141 = zz(i,1) - zz(i,14)
1932
1933 x142 = xx(i,2) - xx(i,14)
1934 y142 = yy(i,2) - yy(i,14)
1935 z142 = zz(i,2) - zz(i,14)
1936
1937 x152 = xx(i,2) - xx(i,15)
1938 y152 = yy(i,2) - yy(i,15)
1939 z152 = zz(i,2) - zz(i,15)
1940
1941 x153 = xx(i,3) - xx(i,15)
1942 y153 = yy(i,3) - yy(i,15)
1943 z153 = zz(i,3) - zz(i,15)
1944
1945 x163 = xx(i,3) - xx(i,16)
1946 y163 = yy(i,3) - yy(i,16)
1947 z163 = zz(i,3) - zz(i,16)
1948
1949 x164 = xx(i,4) - xx(i,16)
1950 y164 = yy(i,4) - yy(i,16)
1951 z164 = zz(i,4) - zz(i,16)
1952
1953 x174 = xx(i,4) - xx(i,17)
1954 y174 = yy(i,4) - yy(i,17)
1955 z174 = zz(i,4) - zz(i,17)
1956
1957 x171 = xx(i,1) - xx(i,17)
1958 y171 = yy(i,1) - yy(i,17)
1959 z171 = zz(i,1) - zz(i,17)
1960
1961 IF(mvoisin(1,cand_e(i))/=0)THEN
1962 sx2114(i) = y142*z141 - z142*y141
1963 sy2114(i) = z142*x141 - x142*z141
1964 sz2114(i) = x142*y141 - y142*x141
1965 ELSE
1966 sx2114(i) = sx125(i)
1967 sy2114(i) = sy125(i)
1968 sz2114(i) = sz125(i)
1969 ENDIF
1970
1971 IF(mvoisin(2,cand_e(i))/=0)THEN
1972 sx3215(i) = y153*z152 - z153*y152
1973 sy3215(i) = z153*x152 - x153*z152
1974 sz3215(i) = x153*y152 - y153*x152
1975 ELSE
1976 sx3215(i) = sx235(i)
1977 sy3215(i) = sy235(i)
1978 sz3215(i) = sz235(i)
1979 ENDIF
1980
1981 IF(mvoisin(3,cand_e(i))/=0)THEN
1982 sx4316(i) = y164*z163 - z164*y163
1983 sy4316(i) = z164*x163 - x164*z163
1984 sz4316(i) = x164*y163 - y164*x163
1985 ELSE
1986 sx4316(i) = sx345(i)
1987 sy4316(i) = sy345(i)
1988 sz4316(i) = sz345(i)
1989 ENDIF
1990
1991 IF(mvoisin(4,cand_e(i))/=0)THEN
1992 sx1417(i) = y171*z174 - z171*y174
1993 sy1417(i) = z171*x174 - x171*z174
1994 sz1417(i) = x171*y174 - y171*x174
1995 ELSE
1996 sx1417(i) = sx415(i)
1997 sy1417(i) = sy415(i)
1998 sz1417(i) = sz415(i)
1999 ENDIF
2000
2001 ELSE
2002
2003 xx(i,6:17) = xx0(i,6:17)
2004 yy(i,6:17) = yy0(i,6:17)
2005 zz(i,6:17) = zz0(i,6:17)
2006
2007 sx2114(i) = sx21140(i)
2008 sy2114(i) = sy21140(i)
2009 sz2114(i) = sz21140(i)
2010
2011 sx3215(i) = sx32150(i)
2012 sy3215(i) = sy32150(i)
2013 sz3215(i) = sz32150(i)
2014
2015 sx4316(i) = sx43160(i)
2016 sy4316(i) = sy43160(i)
2017 sz4316(i) = sz43160(i)
2018
2019 sx1417(i) = sx14170(i)
2020 sy1417(i) = sy14170(i)
2021 sz1417(i) = sz14170(i)
2022
2023 ENDIF
2024c write(iout,*)'i=',i,'i24dst3 212'
2025 ENDDO
2026
2027
2028
2029C=======================================================================
2030c 3 first sliding
2031C=======================================================================
2032c write(iout,*)'i24dst3 6'
2033 nslid = 0
2034 DO i=1,jlt
2035 IF(istep(i) /= 3)cycle
2036 istep(i) = 5
2037! +------+-------------------------------------------+
2038! | | sous-triangle voisins |
2039! |sous +----------+----------+---------------------+
2040! |trian | n3/=n4 | n3=n4 |
2041! | +----------+----------+----------+----------+
2042! | | IT0 | 2,4,6,8=0| | 2,4,6,8=0|
2043! +------+----------+----------+----------+----------+
2044! | 1 | 5 2 4 | 5 2 4 | 5 6 8 | 5 6 8 |
2045! | 2 | 6 3 1 | 6 3 1 | - - - | - - - |
2046! | 3 | 7 4 2 | 7 4 2 | - - - | - - - |
2047! | 4 | 8 1 3 | 8 1 3 | - - - | - - - |
2048! +------+----------+----------+----------+----------+
2049! | 5 | 1 9 13 | 1 -1 -1 | 1 9 13 | 1 -1 -1 |
2050! | 6 | 2 10 14 | 2 -1 -1 | 1 10 14 | 1 -1 -1 |
2051! | 7 | 3 11 15 | 3 -1 -1 | - - - | - - - |
2052! | 8 | 4 12 16 | 4 -1 -1 | 1 12 16 | 1 -1 -1 |
2053! +------+----------+----------+----------+----------+
2054! | 9 | -1 17 5 | - - - | -1 17 5 | - - - |
2055! | 10 | -1 18 6 | - - - | -1 18 6 | - - - |
2056! | 11 | -1 19 7 | - - - | - - - | - - - |
2057! | 12 | -1 20 8 | - - - | -1 20 8 | - - - |
2058! +------+----------+----------+----------+----------+
2059! | 13 | -1 5 17 | - - - | -1 5 17 | - - - |
2060! | 14 | -1 6 18 | - - - | -1 6 18 | - - - |
2061! | 15 | -1 7 19 | - - - | - - - | - - - |
2062! | 16 | -1 8 20 | - - - | -1 8 20 | - - - |
2063! +------+----------+----------+----------+----------+
2064! | 17 | -1 9 13 | - - - | -1 9 13 | - - - |
2065! | 18 | -1 10 14 | - - - | -1 10 14 | - - - |
2066! | 19 | -1 11 15 | - - - | - - - | - - - |
2067! | 20 | -1 12 16 | - - - | -1 12 16 | - - - |
2068! +------+----------+----------+----------+----------+
2069 it(1:3,1:20,i)=it0(1:3,1:20)
2070
2071 IF(ixx(i,3) /= ixx(i,4))THEN
2072 IF(mvoisin(1,cand_e(i))==0)it( 1, 1,i) = 0
2073 IF(mvoisin(2,cand_e(i))==0)it( 1, 2,i) = 0
2074 IF(mvoisin(3,cand_e(i))==0)it( 1, 3,i) = 0
2075 IF(mvoisin(4,cand_e(i))==0)it( 1, 4,i) = 0
2076 ELSE
2077 it(2,1,i)=6
2078 it(3,1,i)=8
2079 it(1,5,i)=1
2080 it(1,6,i)=1
2081 it(1,8,i)=1
2082 IF(mvoisin(1,cand_e(i))==0)it( 1, 1,i) = 0
2083 IF(mvoisin(2,cand_e(i))==0)it( 2, 1,i) = 0
2084 IF(mvoisin(3,cand_e(i))==0)it( 3, 1,i) = 0
2085c a vrifier 4,IF(MVOISIN(4,CAND_E(I))==0)IT( 3, 1,I) = 0
2086 ENDIF
2087 IF(ixx(i,6)==ixx(i,7))THEN
2088 it(2,5,i)=-1
2089 it(3,5,i)=-1
2090 ENDIF
2091 IF(ixx(i,8)==ixx(i,9))THEN
2092 it(2,6,i)=-1
2093 it(3,6,i)=-1
2094 ENDIF
2095 IF(ixx(i,10)==ixx(i,11))THEN
2096 it(2,7,i)=-1
2097 it(3,7,i)=-1
2098 ENDIF
2099 IF(ixx(i,12)==ixx(i,13))THEN
2100 it(2,8,i)=-1
2101 it(3,8,i)=-1
2102 ENDIF
2103
2104 IF(ixx(i,3) /= ixx(i,4))THEN
2105
2106 IF(subtria(i)==1)THEN
2107 xa(i) = xx(i,5) ! f-----------a-----------e
2108 ya(i) = yy(i,5) ! \ / \ /
2109 za(i) = zz(i,5) ! \ T / \ S /
2110 xb(i) = xx(i,1) ! \ / \ /
2111 yb(i) = yy(i,1) ! \ / Q \ /
2112 zb(i) = zz(i,1) ! \ / \ /
2113 xc(i) = xx(i,2) ! b-----------c
2114 yc(i) = yy(i,2) ! \ /
2115 zc(i) = zz(i,2) ! \ R /
2116 ! \ /
2117 nqx = sx125(i) ! \ /
2118 nqy = sy125(i) ! \ /
2119 nqz = sz125(i) ! d
2120 nsx = sx235(i)
2121 nsy = sy235(i)
2122 nsz = sz235(i)
2123 ntx = sx415(i)
2124 nty = sy415(i)
2125 ntz = sz415(i)
2126 nrx = sx2114(i)
2127 nry = sy2114(i)
2128 nrz = sz2114(i)
2129 itr = 5
2130 IF(mvoisin(1,cand_e(i))==0)THEN
2131 itr = -itr
2132 nrx = nqx
2133 nry = nqy
2134 nrz = nqz
2135 ENDIF
2136 its = 2
2137 itt = 4
2138 ELSEIF(subtria(i)==2)THEN
2139 xa(i) = xx(i,5)
2140 ya(i) = yy(i,5)
2141 za(i) = zz(i,5)
2142 xb(i) = xx(i,2)
2143 yb(i) = yy(i,2)
2144 zb(i) = zz(i,2)
2145 xc(i) = xx(i,3)
2146 yc(i) = yy(i,3)
2147 zc(i) = zz(i,3)
2148 nqx = sx235(i)
2149 nqy = sy235(i)
2150 nqz = sz235(i)
2151 nsx = sx345(i)
2152 nsy = sy345(i)
2153 nsz = sz345(i)
2154 ntx = sx125(i)
2155 nty = sy125(i)
2156 ntz = sz125(i)
2157 nrx = sx3215(i)
2158 nry = sy3215(i)
2159 nrz = sz3215(i)
2160 itr = 6
2161 IF(mvoisin(2,cand_e(i))==0)THEN
2162 itr = -itr
2163 nrx = nqx
2164 nry = nqy
2165 nrz = nqz
2166 ENDIF
2167 its = 3
2168 itt = 1
2169 ELSEIF(subtria(i)==3)THEN
2170 xa(i) = xx(i,5)
2171 ya(i) = yy(i,5)
2172 za(i) = zz(i,5)
2173 xb(i) = xx(i,3)
2174 yb(i) = yy(i,3)
2175 zb(i) = zz(i,3)
2176 xc(i) = xx(i,4)
2177 yc(i) = yy(i,4)
2178 zc(i) = zz(i,4)
2179 nqx = sx345(i)
2180 nqy = sy345(i)
2181 nqz = sz345(i)
2182 nsx = sx415(i)
2183 nsy = sy415(i)
2184 nsz = sz415(i)
2185 ntx = sx235(i)
2186 nty = sy235(i)
2187 ntz = sz235(i)
2188 nrx = sx4316(i)
2189 nry = sy4316(i)
2190 nrz = sz4316(i)
2191 itr = 7
2192 IF(mvoisin(3,cand_e(i))==0)THEN
2193 itr = -itr
2194 nrx = nqx
2195 nry = nqy
2196 nrz = nqz
2197 ENDIF
2198 its = 4
2199 itt = 2
2200 ELSEIF(subtria(i)==4)THEN
2201 xa(i) = xx(i,5)
2202 ya(i) = yy(i,5)
2203 za(i) = zz(i,5)
2204 xb(i) = xx(i,4)
2205 yb(i) = yy(i,4)
2206 zb(i) = zz(i,4)
2207 xc(i) = xx(i,1)
2208 yc(i) = yy(i,1)
2209 zc(i) = zz(i,1)
2210 nqx = sx415(i)
2211 nqy = sy415(i)
2212 nqz = sz415(i)
2213 nsx = sx125(i)
2214 nsy = sy125(i)
2215 nsz = sz125(i)
2216 ntx = sx345(i)
2217 nty = sy345(i)
2218 ntz = sz345(i)
2219 nrx = sx1417(i)
2220 nry = sy1417(i)
2221 nrz = sz1417(i)
2222 itr = 8
2223 IF(mvoisin(4,cand_e(i))==0)THEN
2224 itr = -itr
2225 nrx = nqx
2226 nry = nqy
2227 nrz = nqz
2228 ENDIF
2229 its = 1
2230 itt = 3
2231 ENDIF
2232 else!IF(IXX(I,3) /= IXX(I,4)))THEN
2233 xa(i) = xx(i,3)
2234 ya(i) = yy(i,3)
2235 za(i) = zz(i,3)
2236 xb(i) = xx(i,1) ! f-----------a-----------e
2237 yb(i) = yy(i,1) ! \ / \ /
2238 zb(i) = zz(i,1) ! \ T / \ S /
2239 xc(i) = xx(i,2) ! \ / \ /
2240 yc(i) = yy(i,2) ! \ / Q \ /
2241 zc(i) = zz(i,2) ! \ / \ /
2242 nqx = sx125(i) ! b-----------c
2243 nqy = sy125(i) ! \ /
2244 nqz = sz125(i) ! \ R /
2245 nsx = sx3215(i) ! \ /
2246 nsy = sy3215(i) ! \ /
2247 nsz = sz3215(i) ! \ /
2248 ntx = sx1417(i) ! d
2249 nty = sy1417(i)
2250 ntz = sz1417(i)
2251 nrx = sx2114(i)
2252 nry = sy2114(i)
2253 nrz = sz2114(i)
2254 itr = 5
2255 IF(mvoisin(1,cand_e(i))==0)THEN
2256 itr = -itr
2257 nrx = nqx
2258 nry = nqy
2259 nrz = nqz
2260 ENDIF
2261 its = 6
2262 IF(mvoisin(2,cand_e(i))==0)THEN
2263 its = -its
2264 nsx = nqx
2265 nsy = nqy
2266 nsz = nqz
2267 ENDIF
2268 itt = 8
2269 IF(mvoisin(4,cand_e(i))==0)THEN
2270 itt = -itt
2271 ntx = nqx
2272 nty = nqy
2273 ntz = nqz
2274 ENDIF
2275 ENDIF
2276
2277 nx(i) = nqx
2278 ny(i) = nqy
2279 nz(i) = nqz
2280
2281c calcul des plans bisecteurs
2282 nqrx = nrx + nqx
2283 nqry = nry + nqy
2284 nqrz = nrz + nqz
2285 nqsx = nsx + nqx
2286 nqsy = nsy + nqy
2287 nqsz = nsz + nqz
2288 nqtx = ntx + nqx
2289 nqty = nty + nqy
2290 nqtz = ntz + nqz
2291 IF(abs(itr)==subtria_n(i))THEN
2292 prx = nqry*(zc(i)-zb(i)) - nqrz*(yc(i)-yb(i))
2293 pry = nqrz*(xc(i)-xb(i)) - nqrx*(zc(i)-zb(i))
2294 prz = nqrx*(yc(i)-yb(i)) - nqry*(xc(i)-xb(i))
2295 vr = prx*(xi(i)-xb(i))+ pry*(yi(i)-yb(i))+ prz*(zi(i)-zb(i))
2296 IF(vr < zero)THEN
2297 IF(itr>0)THEN
2298 istep(i) = 5
2299 nslid = nslid+1
2300 subtria(i) = itr
2301 ELSE
2302c sliding out of the surface
2303 icontact(i) = 0
2304 istep(i) = 0
2305 subtria(i) = 0
2306 ENDIF
2307 nx(i) = nrx
2308 ny(i) = nry
2309 nz(i) = nrz
2310 ENDIF
2311 ELSEIF(abs(its)==subtria_n(i))THEN
2312 psx = nqsy*(za(i)-zc(i)) - nqsz*(ya(i)-yc(i))
2313 psy = nqsz*(xa(i)-xc(i)) - nqsx*(za(i)-zc(i))
2314 psz = nqsx*(ya(i)-yc(i)) - nqsy*(xa(i)-xc(i))
2315 vs = psx*(xi(i)-xc(i))+ psy*(yi(i)-yc(i))+ psz*(zi(i)-zc(i))
2316 IF(vs < zero)THEN
2317 IF(its>0)THEN
2318 istep(i) = 5
2319 nslid = nslid+1
2320 subtria(i) = its
2321 ELSE
2322c sliding out of the surface
2323 icontact(i) = 0
2324 istep(i) = 0
2325 subtria(i) = 0
2326 ENDIF
2327 nx(i) = nsx
2328 ny(i) = nsy
2329 nz(i) = nsz
2330 ENDIF
2331 ELSEIF(abs(itt)==subtria_n(i))THEN
2332 ptx = nqty*(zb(i)-za(i)) - nqtz*(yb(i)-ya(i))
2333 pty = nqtz*(xb(i)-xa(i)) - nqtx*(zb(i)-za(i))
2334 ptz = nqtx*(yb(i)-ya(i)) - nqty*(xb(i)-xa(i))
2335 vt = ptx*(xi(i)-xa(i))+ pty*(yi(i)-ya(i))+ ptz*(zi(i)-za(i))
2336 IF(vt < zero)THEN
2337 IF(itt>0)THEN
2338 istep(i) = 5
2339 nslid = nslid+1
2340 subtria(i) = itt
2341 ELSE
2342c sliding out of the surface
2343 icontact(i) = 0
2344 istep(i) = 0
2345 subtria(i) = 0
2346 ENDIF
2347 nx(i) = ntx
2348 ny(i) = nty
2349 nz(i) = ntz
2350 ENDIF
2351 ENDIF
2352 ENDDO
2353C=======================================================================
2354c 4 next sliding
2355C=======================================================================
2356c write(iout,*)'NSLID=',NSLID,'i24dst3 7'
2357 nss = 1
2358!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2359 nslid=0 ! debranchement
2360!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2361C=======================================================================
2362c 5 surface coordinates
2363C=======================================================================
2364#include "lockon.inc"
2365c write(iout,*)'i24dst3 8'
2366 DO i=1,jlt
2367 IF(istep(i) /= 1)pene(i) = zero
2368 IF(istep(i) < 4 .OR.istep(i) == 6)cycle
2369 itq = subtria(i)
2370 k = ic(1,itq)
2371 xa(i) = xx(i,k)
2372 ya(i) = yy(i,k)
2373 za(i) = zz(i,k)
2374 nax = nxx(i,k)
2375 nay = nyy(i,k)
2376 naz = nzz(i,k)
2377
2378 k = ic(2,itq)
2379 xb(i) = xx(i,k)
2380 yb(i) = yy(i,k)
2381 zb(i) = zz(i,k)
2382 nbx = nxx(i,k)
2383 nby = nyy(i,k)
2384 nbz = nzz(i,k)
2385
2386 k = ic(3,itq)
2387 xc(i) = xx(i,k)
2388 yc(i) = yy(i,k)
2389 zc(i) = zz(i,k)
2390 ncx = nxx(i,k)
2391 ncy = nyy(i,k)
2392 ncz = nzz(i,k)
2393
2394 nqx = nx(i)
2395 nqy = ny(i)
2396 nqz = nz(i)
2397
2398 aaa = max(em30,sqrt(nqx*nqx + nqy*nqy + nqz*nqz))
2399 area(i) = aaa
2400 s2 = one/aaa
2401 nqx = nqx*s2
2402 nqy = nqy*s2
2403 nqz = nqz*s2
2404
2405 bbb = (xi(i)-xa(i))*nqx+(yi(i)-ya(i))*nqy+(zi(i)-za(i))*nqz
2406 ikeep=0
2407C---------------change NSNE to global flag after---
2408C------debug P/ON RD-5326 limiting this no-physical case (remaining contact)----
2409 IF (bbb>=zero .AND.bbb<marge025.AND.eps>em02 .AND. impl_s==0
2410 + .AND. nsne==0) THEN
2411 pene(i) = em15
2412 ikeep=1
2413 ELSE
2414 pene(i) = -min(zero,bbb)
2415 END IF
2416c PMAX_GAP = MAX(PMAX_GAP,ABS(PENE(I)))
2417
2418C--------exclude high pene in cas of ICONT_R
2419 IF(icont_r(i) >0)THEN
2420 tole =eps02*aaa
2421 gap2 = gaps(i)+gap_m(cand_e(i))
2422 IF (gap2 >zero) tole =min(tole,gap2*gap2)
2423 IF(pene(i)*pene(i) > tole) THEN
2424 pene(i) = zero
2425 END IF
2426 END if!(ICONT_R(I) >0)THEN
2427
2428 IF(pene(i) == zero)THEN
2429 icontact(i) = 0
2430 cycle
2431 ENDIF
2432
2433c N1(I) = NQX
2434c N2(I) = NQY
2435c N3(I) = NQZ
2436
2437c project triangle on SECONDARY node
2438
2439
2440 ! a
2441 ! / \
2442 ! / \ a,b,c position at current time
2443 ! / A \ A,B,C projected on SECONDARY node
2444 ! / /\ \
2445 ! / / \ \
2446 ! / / \ \ normal at point a:
2447 ! / / \ \ Na = a - A
2448 ! / / \ \
2449 ! / / \ \
2450 ! / / \ \
2451 ! / b--------------c \
2452 ! / \
2453 ! b-------------------------c
2454 !
2455
2456 nni = (nqx*nax + nqy*nay + nqz*naz)
2457 ni2 = nax*nax + nay*nay + naz*naz
2458 IF(two*nni*nni < ni2)THEN
2459c scharp angle bound nodal normal to 45 from segment normal
2460 aaa = sqrt(ni2-nni*nni) - nni
2461 nax = nax + aaa*nqx
2462 nay = nay + aaa*nqy
2463 naz = naz + aaa*nqz
2464 nni = (nqx*nax + nqy*nay + nqz*naz)
2465 ENDIF
2466
2467 aaa = pene(i)/nni
2468 nax = nax*aaa
2469 nay = nay*aaa
2470 naz = naz*aaa
2471
2472 xpa = xa(i) - nax
2473 ypa = ya(i) - nay
2474 zpa = za(i) - naz
2475
2476 nni = (nqx*nbx + nqy*nby + nqz*nbz)
2477 ni2 = nbx*nbx + nby*nby + nbz*nbz
2478 IF(two*nni*nni < ni2)THEN
2479c scharp angle bound nodal normal to 45 from segment normal
2480 aaa = sqrt(ni2-nni*nni) - nni
2481 nbx = nbx + aaa*nqx
2482 nby = nby + aaa*nqy
2483 nbz = nbz + aaa*nqz
2484 nni = (nqx*nbx + nqy*nby + nqz*nbz)
2485 ENDIF
2486
2487 aaa = pene(i)/nni
2488 nbx = nbx*aaa
2489 nby = nby*aaa
2490 nbz = nbz*aaa
2491
2492 xpb = xb(i) - nbx
2493 ypb = yb(i) - nby
2494 zpb = zb(i) - nbz
2495
2496 nni = (nqx*ncx + nqy*ncy + nqz*ncz)
2497 ni2 = ncx*ncx + ncy*ncy + ncz*ncz
2498 IF(two*nni*nni < ni2)THEN
2499c scharp angle bound nodal normal to 45 from segment normal
2500 aaa = sqrt(ni2-nni*nni) - nni
2501 ncx = ncx + aaa*nqx
2502 ncy = ncy + aaa*nqy
2503 ncz = ncz + aaa*nqz
2504 nni = (nqx*ncx + nqy*ncy + nqz*ncz)
2505 ENDIF
2506
2507 aaa = pene(i)/nni
2508 ncx = ncx*aaa
2509 ncy = ncy*aaa
2510 ncz = ncz*aaa
2511
2512 xpc = xc(i) - ncx
2513 ypc = yc(i) - ncy
2514 zpc = zc(i) - ncz
2515
2516 nnx=(ypb-ypa)*(zpc-zpa) - (zpb-zpa)*(ypc-ypa)
2517 nny=(zpb-zpa)*(xpc-xpa) - (xpb-xpa)*(zpc-zpa)
2518 nnz=(xpb-xpa)*(ypc-ypa) - (ypb-ypa)*(xpc-xpa)
2519 aaa = nnx*nqx + nny*nqy + nnz*nqz
2520 IF(aaa <= zero)THEN
2521 nax = nqx*pene(i)
2522 nay = nqy*pene(i)
2523 naz = nqz*pene(i)
2524 nbx = nax
2525 nby = nay
2526 nbz = naz
2527 ncx = nax
2528 ncy = nay
2529 ncz = naz
2530 xpa = xa(i) - nax
2531 ypa = ya(i) - nay
2532 zpa = za(i) - naz
2533 xpb = xb(i) - nbx
2534 ypb = yb(i) - nby
2535 zpb = zb(i) - nbz
2536 xpc = xc(i) - ncx
2537 ypc = yc(i) - ncy
2538 zpc = zc(i) - ncz
2539 ENDIF
2540
2541 xab = xpb-xpa
2542 yab = ypb-ypa
2543 zab = zpb-zpa
2544 xbc = xpc-xpb
2545 ybc = ypc-ypb
2546 zbc = zpc-zpb
2547 xca = xpa-xpc
2548 yca = ypa-ypc
2549 zca = zpa-zpc
2550
2551 xia = xpa-xi(i)
2552 yia = ypa-yi(i)
2553 zia = zpa-zi(i)
2554 xib = xpb-xi(i)
2555 yib = ypb-yi(i)
2556 zib = zpb-zi(i)
2557 xic = xpc-xi(i)
2558 yic = ypc-yi(i)
2559 zic = zpc-zi(i)
2560 sx = - yab*zca + zab*yca
2561 sy = - zab*xca + xab*zca
2562 sz = - xab*yca + yab*xca
2563 s2 = sx*sx+sy*sy+sz*sz
2564 sax = yib*zic - zib*yic
2565 say = zib*xic - xib*zic
2566 saz = xib*yic - yib*xic
2567 la(i) = (sx*sax+sy*say+sz*saz)/s2
2568 sbx = yic*zia - zic*yia
2569 sby = zic*xia - xic*zia
2570 sbz = xic*yia - yic*xia
2571 lb(i) = (sx*sbx+sy*sby+sz*sbz)/s2
2572 lc(i) = one - la(i) - lb(i)
2573
2574 IF(la(i)>=zero.and.lb(i)>=zero.and.lc(i)>=zero)THEN
2575c interpolation des normales
2576c avec surponderation du point central (quadrangles)
2577c avec surponderation de la facette (triangles)
2578c goto 2345
2579c add ici
2580 2345 continue
2581 IF (ikeep==1) THEN
2582 pene(i) = zero
2583 icontact(i) = 0
2584 cycle
2585 END IF
2586 n1(i) = nqx
2587 n2(i) = nqy
2588 n3(i) = nqz
2589 ELSE
2590C--------return to seg before sliding +
2591 bbb = pene_tm(i)-pene(i)
2592 IF (bbb>tol_b.AND.la(i)<zero) THEN
2593 pene(i) = pene_tm(i)
2594 itq = subtria_tm(i)
2595C---------for implicit
2596 subtria(i) = itq
2597 ix1(i) = ixx(i,ic( 7,itq))
2598 ix2(i) = ixx(i,ic( 8,itq))
2599 ix3(i) = ixx(i,ic( 9,itq))
2600 ix4(i) = ixx(i,ic(10,itq))
2601 x1(i) = xx0(i,ic( 7,itq))
2602 x2(i) = xx0(i,ic( 8,itq))
2603 x3(i) = xx0(i,ic( 9,itq))
2604 x4(i) = xx0(i,ic(10,itq))
2605 y1(i) = yy0(i,ic( 7,itq))
2606 y2(i) = yy0(i,ic( 8,itq))
2607 y3(i) = yy0(i,ic( 9,itq))
2608 y4(i) = yy0(i,ic(10,itq))
2609 z1(i) = zz0(i,ic( 7,itq))
2610 z2(i) = zz0(i,ic( 8,itq))
2611 z3(i) = zz0(i,ic( 9,itq))
2612 z4(i) = zz0(i,ic(10,itq))
2613
2614 vx1(i) = vx(i,ic( 7,itq))
2615 vx2(i) = vx(i,ic( 8,itq))
2616 vx3(i) = vx(i,ic( 9,itq))
2617 vx4(i) = vx(i,ic(10,itq))
2618 vy1(i) = vy(i,ic( 7,itq))
2619 vy2(i) = vy(i,ic( 8,itq))
2620 vy3(i) = vy(i,ic( 9,itq))
2621 vy4(i) = vy(i,ic(10,itq))
2622 vz1(i) = vz(i,ic( 7,itq))
2623 vz2(i) = vz(i,ic( 8,itq))
2624 vz3(i) = vz(i,ic( 9,itq))
2625 vz4(i) = vz(i,ic(10,itq))
2626
2627 la(i) = la_tm(i)
2628 lb(i) = lb_tm(i)
2629 lc(i) = one - la(i) - lb(i)
2630 IF(la(i)<zero)THEN
2631 IF(lb(i)<zero)THEN
2632 la(i) = zero
2633 lb(i) = zero
2634 lc(i) = one
2635 ELSEIF(lc(i)<zero)THEN
2636 lc(i) = zero
2637 la(i) = zero
2638 lb(i) = one
2639 ELSE
2640 la(i) = zero
2641 aaa = lb(i) + lc(i)
2642 lb(i) = lb(i)/aaa
2643 lc(i) = lc(i)/aaa
2644 ENDIF
2645 ELSEIF(lb(i)<zero)THEN
2646 IF(lc(i)<zero)THEN
2647 lb(i) = zero
2648 lc(i) = zero
2649 la(i) = one
2650 ELSE
2651 lb(i) = zero
2652 aaa = lc(i) + la(i)
2653 lc(i) = lc(i)/aaa
2654 la(i) = la(i)/aaa
2655 ENDIF
2656 ELSEIF(lc(i)<zero)THEN
2657 lc(i) = zero
2658 aaa = la(i) + lb(i)
2659 la(i) = la(i)/aaa
2660 lb(i) = lb(i)/aaa
2661 ENDIF
2662
2663 IF (ix1(i) == ix2(i))THEN
2664 h1(i) = la(i)
2665 h2(i) = zero
2666 h3(i) = lb(i)
2667 h4(i) = lc(i)
2668 ELSEIF(ix2(i) == ix3(i))THEN
2669 h1(i) = la(i)
2670 h2(i) = lb(i)
2671 h3(i) = zero
2672 h4(i) = lc(i)
2673c ELSEIF(IX3(I) == IX4(I))THEN impossible
2674 ELSEIF(ix4(i) == ix1(i))THEN
2675 h1(i) = lc(i)
2676 h2(i) = la(i)
2677 h3(i) = lb(i)
2678 h4(i) = zero
2679 ELSE
2680 h0 = fourth * la(i)
2681 h1(i) = h0
2682 h2(i) = h0
2683 h3(i) = h0 + lb(i)
2684 h4(i) = h0 + lc(i)
2685 ENDIF
2686 n1(i) = n_tm(1,i)
2687 n2(i) = n_tm(2,i)
2688 n3(i) = n_tm(3,i)
2689 cycle
2690 END IF
2691 n1(i) = nqx
2692 n2(i) = nqy
2693 n3(i) = nqz
2694 IF(la(i)<zero.and.lb(i)<zero)THEN
2695 la(i) = zero
2696 lb(i) = zero
2697 lc(i) = one
2698 ELSEIF(lb(i)<zero.and.lc(i)<zero)THEN
2699 lb(i) = zero
2700 lc(i) = zero
2701 la(i) = one
2702 ELSEIF(lc(i)<zero.and.la(i)<zero)THEN
2703 lc(i) = zero
2704 la(i) = zero
2705 lb(i) = one
2706 ELSEIF(la(i)<zero)THEN
2707 la(i) = zero
2708 aaa = lb(i) + lc(i)
2709 lb(i) = lb(i)/aaa
2710 lc(i) = lc(i)/aaa
2711 ELSEIF(lb(i)<zero)THEN
2712 lb(i) = zero
2713 aaa = lc(i) + la(i)
2714 lc(i) = lc(i)/aaa
2715 la(i) = la(i)/aaa
2716 ELSEIF(lc(i)<zero)THEN
2717 lc(i) = zero
2718 aaa = la(i) + lb(i)
2719 la(i) = la(i)/aaa
2720 lb(i) = lb(i)/aaa
2721 ENDIF
2722 goto 3456
2723 xpi = la(i)*xa(i) + lb(i)*xb(i) + lc(i)*xc(i)
2724 ypi = la(i)*ya(i) + lb(i)*yb(i) + lc(i)*yc(i)
2725 zpi = la(i)*za(i) + lb(i)*zb(i) + lc(i)*zc(i)
2726 n1(i) = xpi - xi(i)
2727 n2(i) = ypi - yi(i)
2728 n3(i) = zpi - zi(i)
2729 pene(i) = sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2)
2730 s2 = one/max(em30,pene(i))
2731 n1(i) = n1(i)*s2
2732 n2(i) = n2(i)*s2
2733 n3(i) = n3(i)*s2
2734 3456 continue
2735 ENDIF
2736c interpolation des normales seulement en glissement et en option avec implicite
2737C------issue w/ sp too small value of S2----
2738 IF ((impl_s>0 .AND. ittoff>0)) THEN
2739C IF (ISPT2(I)>0.OR.(IMPL_S>0 .AND. ITTOFF>0)) THEN
2740c IF ((ITQ==1 .or. ITQ==2 .or. ITQ==3 .or. ITQ==4).and.
2741c . (IXX(I,3) /= IXX(I,4)))THEN
2742 overw = overw0
2743 n1(i) = la(i)*nax + lb(i)*nbx + lc(i)*ncx
2744 n2(i) = la(i)*nay + lb(i)*nby + lc(i)*ncy
2745 n3(i) = la(i)*naz + lb(i)*nbz + lc(i)*ncz
2746c ELSE
2747c OVERW = OVERW0*LA(I)*LB(I)*LC(I)
2748c N1(I) = OVERW*NQX + LA(I)*NAX + LB(I)*NBX + LC(I)*NCX
2749c N2(I) = OVERW*NQY + LA(I)*NAY + LB(I)*NBY + LC(I)*NCY
2750c N3(I) = OVERW*NQZ + LA(I)*NAZ + LB(I)*NBZ + LC(I)*NCZ
2751c ENDIF
2752
2753 s2 = one/max(em30,sqrt(n1(i)**2 + n2(i)**2 + n3(i)**2))
2754 n1(i) = n1(i)*s2
2755 n2(i) = n2(i)*s2
2756 n3(i) = n3(i)*s2
2757
2758 aaa = nqx*n1(i) + nqy*n2(i) + nqz*n3(i)
2759 aaa = max(one,one/max(aaa,em20))
2760 pene(i) = pene(i)*aaa
2761
2762c if(pene(i) > 1.e17)then
2763c ibug=1
2764c stop 9876
2765c endif
2766
2767 END IF !(IMPL_S>0 .AND. ITTOFF>0) THEN
2768
2769 ix1(i) = ixx(i,ic( 7,itq))
2770 ix2(i) = ixx(i,ic( 8,itq))
2771 ix3(i) = ixx(i,ic( 9,itq))
2772 ix4(i) = ixx(i,ic(10,itq))
2773
2774 x1(i) = xx0(i,ic( 7,itq))
2775 x2(i) = xx0(i,ic( 8,itq))
2776 x3(i) = xx0(i,ic( 9,itq))
2777 x4(i) = xx0(i,ic(10,itq))
2778 y1(i) = yy0(i,ic( 7,itq))
2779 y2(i) = yy0(i,ic( 8,itq))
2780 y3(i) = yy0(i,ic( 9,itq))
2781 y4(i) = yy0(i,ic(10,itq))
2782 z1(i) = zz0(i,ic( 7,itq))
2783 z2(i) = zz0(i,ic( 8,itq))
2784 z3(i) = zz0(i,ic( 9,itq))
2785 z4(i) = zz0(i,ic(10,itq))
2786
2787 vx1(i) = vx(i,ic( 7,itq))
2788 vx2(i) = vx(i,ic( 8,itq))
2789 vx3(i) = vx(i,ic( 9,itq))
2790 vx4(i) = vx(i,ic(10,itq))
2791 vy1(i) = vy(i,ic( 7,itq))
2792 vy2(i) = vy(i,ic( 8,itq))
2793 vy3(i) = vy(i,ic( 9,itq))
2794 vy4(i) = vy(i,ic(10,itq))
2795 vz1(i) = vz(i,ic( 7,itq))
2796 vz2(i) = vz(i,ic( 8,itq))
2797 vz3(i) = vz(i,ic( 9,itq))
2798 vz4(i) = vz(i,ic(10,itq))
2799
2800 IF (ix1(i) == ix2(i))THEN
2801 h1(i) = la(i)
2802 h2(i) = zero
2803 h3(i) = lb(i)
2804 h4(i) = lc(i)
2805 ELSEIF(ix2(i) == ix3(i))THEN
2806 h1(i) = la(i)
2807 h2(i) = lb(i)
2808 h3(i) = zero
2809 h4(i) = lc(i)
2810c ELSEIF(IX3(I) == IX4(I))THEN impossible
2811 ELSEIF(ix4(i) == ix1(i))THEN
2812 h1(i) = lc(i)
2813 h2(i) = la(i)
2814 h3(i) = lb(i)
2815 h4(i) = zero
2816 ELSE
2817 h0 = fourth * la(i)
2818 h1(i) = h0
2819 h2(i) = h0
2820 h3(i) = h0 + lb(i)
2821 h4(i) = h0 + lc(i)
2822 ENDIF
2823C IF (INCONV /= 1) CYCLE
2824
2825 n=cand_n(i)
2826
2827! if(nsvg(i) == ix1(i).or.nsvg(i) == ix2(i).or.
2828! . nsvg(i) == ix3(i).or.nsvg(i) == ix4(i))then
2829! write(iout,*)'1 t=',tt,'i=',i,'nsvg(i)',nsvg(i)
2830! write(iout,*)'ix1=',ix1(i),'ix2=',ix2(i)
2831! write(iout,*)'ix3=',ix3(i),'ix4=',ix4(i)
2832! pene(i) = zero
2833! H1(I) = zero
2834! H2(I) = zero
2835! H3(I) = zero
2836! H4(I) = zero
2837! ICONTACT(I) = 0
2838! cycle
2839! endif
2840
2841c IF(N <= NSN)THEN
2842c MGLOB = IRTLM(1,N)
2843c ELSE
2844c MGLOB = IRTLM_FI(NIN)%P(1,N-NSN)
2845c ENDIF
2846c IF(ICONTACT(I)==MGLOB)THEN
2847c ICONTACT(I) == IRTLM > 0 old contact, new segment
2848c ICONTACT(I) == IRTLM < 0 new contact, new segment, retained provisional
2849c ICONTACT(I) /= IRTLM < 0 new contact, new segment, not retained
2850 mglob = mseglo(cand_e(i))
2851 IF (itq==5 .or. itq==9 .or. itq==13 .or. itq==17)THEN
2852 mglob = mvoisin(1,cand_e(i))
2853 itq = itriv(1,i)
2854 ELSEIF(itq==6 .or. itq==10 .or. itq==14 .or. itq==18)THEN
2855 mglob = mvoisin(2,cand_e(i))
2856 itq = itriv(2,i)
2857 ELSEIF(itq==7 .or. itq==11 .or. itq==15 .or. itq==19)THEN
2858 mglob = mvoisin(3,cand_e(i))
2859 itq = itriv(3,i)
2860 ELSEIF(itq==8 .or. itq==12 .or. itq==16 .or. itq==20)THEN
2861 mglob = mvoisin(4,cand_e(i))
2862 itq = itriv(4,i)
2863 ENDIF
2864 IF(n <= nsn)THEN
2865 IF(irtlm(1,n)>0) THEN
2866 irtlm(1,n) = mglob
2867 time_s(n) = -ep20
2868 irtlm(2,n)= itq !sub triangle
2869c ELSEIF((TIME_S(N) < ADT0(I)) .or. ! multi or new intersection
2870c . (TIME_S(N) == ADT0(I) .and. ! multi at same time
2871c . -IRTLM(1,N) < MGLOB ) )THEN
2872c IRTLM(2,N)= ITQ !sub triangle
2873c IRTLM(1,N) = -MGLOB
2874c TIME_S(N) = ADT0(I)
2875 ELSEIF((time_s(n) < pene(i)) .or. ! multi or new intersection
2876 . (time_s(n) == pene(i) .and. ! multi at same time
2877 . -irtlm(1,n) < mglob ) )THEN
2878 irtlm(2,n)= itq !sub triangle
2879 irtlm(1,n) = -mglob
2880 time_s(n) = pene(i)
2881 ENDIF
2882 ELSE ! Remote nodes in SPMD
2883 ii=n-nsn
2884 IF(irtlm_fi(nin)%P(1,ii) > 0)THEN
2885 irtlm_fi(nin)%P(1,ii) = mglob
2886 time_sfi(nin)%P(ii) = -ep20
2887 irtlm_fi(nin)%P(2,ii) = itq
2888 ELSEIF((time_sfi(nin)%P(ii) < pene(i)).or.
2889 . (time_sfi(nin)%P(ii)== pene(i).and.
2890 . -irtlm_fi(nin)%P(1,ii) < mglob) )THEN
2891 irtlm_fi(nin)%P(2,ii) = itq
2892 irtlm_fi(nin)%P(1,ii) = -mglob
2893 time_sfi(nin)%P(ii) = pene(i)
2894 ENDIF
2895 ENDIF
2896 ENDDO
2897
2898c debug
2899!!! goto 123
2900 123 continue
2901C=======================================================================
2902c special plan/plan: node/edge
2903C=======================================================================
2904! old contacts look at if sliding on another edge
2905! +----+----------+----------+-------------+
2906! | IC |
2907! +----+----------+----------+-------------+
2908! | ST | noeuds T | noeuds TV| noeuds Quad |
2909! +----+----------+----------+-------------+ 11-------10
2910! | 1 | 5 1 2 | 14 3 4 | 3 4 1 2 | |\ 19 /|
2911! | 2 | 5 2 3 | 15 4 1 | 4 1 2 3 | | \ / |
2912! | 3 | 5 3 4 | 16 1 2 | 1 2 3 4 | | \ / |
2913! | 4 | 5 4 1 | 17 2 3 | 2 3 4 1 | | 16 |
2914! +----+----------+----------+-------------+ |15/ \11|
2915! | 5 | 14 2 1 | 5 6 7 | 6 7 2 1 | | / \ |
2916! | 6 | 15 3 2 | 5 8 9 | 8 9 3 2 | |/ 7 \|
2917! | 7 | 16 4 3 | 5 10 11 | 10 11 4 3 |12-------4-------3-------9
2918! | 8 | 17 1 4 | 5 12 13 | 12 13 1 4 | |\ 12 /|\ /|\ 14 /|
2919! +----+----------+----------+-------------+ | \ / | \ 3 / | \ / |
2920! | 9 | 14 1 6 |(13) 7 2 | 7 2 1 6 | | \ / | \ /2 |6 \ /18|
2921! | 10 | 15 2 8 |( 7) 9 3 | 9 3 2 8 | | 17 | 5 | 15 |
2922! | 11 | 16 3 10 |( 9)11 4 | 11 4 3 10 | |20/ \ 8| 4/ \ | / \ |
2923! | 12 | 17 4 12 |(11)13 1 | 13 1 4 12 | | / \ | / 1 \ | / \ |
2924! +----+----------+----------+-------------+ |/ 16 \|/ \|/ 10 \|
2925! | 13 | 14 7 2 |( 8) 1 6 | 1 6 7 2 |13-------1-------2-------8
2926! | 14 | 15 9 3 |(10) 2 8 | 2 8 9 3 | |\ 5 /|
2927! | 15 | 16 11 4 |(12) 3 10 | 3 10 11 4 | | \ / |
2928! | 16 | 17 13 1 |( 6) 4 12 | 4 12 13 1 | |9 \ /13|
2929! +----+----------+----------+-------------+ | 14 |
2930! | 17 | 14 6 7 | - 2 1 | 2 1 6 7 | | / \ |
2931! | 18 | 15 8 9 | - 3 2 | 3 2 8 9 | | / \ |
2932! | 19 | 16 10 11 | - 4 3 | 4 3 10 11 | |/ 17 \|
2933! | 20 | 17 12 13 | - 1 4 | 1 4 12 13 | 6-------7
2934! +----+----------+----------+-------------+
2935! sliding of old contacts, SUBTRIA(I)>20: new ones
2936 DO i=1,jlt
2937 IF(istep(i) /= 6 .OR. subtria(i)>20)cycle
2938C old contact ITQ =1-4, edge BC determine if FB BC CG as final one
2939 itq = subtria(i)
2940 IF (ixx(i,4)==ixx(i,3)) itq =1
2941 l = cand_e(i)
2942 k = ic(1,itq) ! 5
2943 xa(i) = xx0(i,k)
2944 ya(i) = yy0(i,k)
2945 za(i) = zz0(i,k)
2946
2947 k = ic(2,itq) ! N1
2948 xb(i) = xx0(i,k)
2949 yb(i) = yy0(i,k)
2950 zb(i) = zz0(i,k)
2951
2952 k = ic(3,itq) ! N2
2953 xc(i) = xx0(i,k)
2954 yc(i) = yy0(i,k)
2955 zc(i) = zz0(i,k)
2956C
2957 k = ic(5,itq) ! N3
2958 xd(i) = xx0(i,k)
2959 yd(i) = yy0(i,k)
2960 zd(i) = zz0(i,k)
2961
2962 k = ic(6,itq) ! N4
2963 xe(i) = xx0(i,k)
2964 ye(i) = yy0(i,k)
2965 ze(i) = zz0(i,k)
2966 xib = xb(i)-xi(i)
2967 yib = yb(i)-yi(i)
2968 zib = zb(i)-zi(i)
2969 xbc = xc(i)-xb(i)
2970 ybc = yc(i)-yb(i)
2971 zbc = zc(i)-zb(i)
2972 aaa = xib*xbc+yib*ybc+zib*zbc
2973 xbd = xd(i)-xb(i)
2974 ybd = yd(i)-yb(i)
2975 zbd = zd(i)-zb(i)
2976 xce = xe(i)-xc(i)
2977 yce = ye(i)-yc(i)
2978 zce = ze(i)-zc(i)
2979 bbb = xib*xbd+yib*ybd+zib*zbd
2980 subtria_n(i) = 0
2981C----- sliding around B
2982 IF (aaa>zero) THEN
2983 IF(itq ==1.AND.mvoisin(4,l)>0)THEN
2984 subtria_n(i) = 16
2985 ELSEIF(itq ==2.AND.mvoisin(1,l)>0)THEN
2986 subtria_n(i) = 13
2987 ELSEIF(itq ==3.AND.mvoisin(2,l)>0)THEN
2988 subtria_n(i) = 14
2989 ELSEIF(itq ==4.AND.mvoisin(3,l)>0)THEN
2990 subtria_n(i) = 15
2991C----- sliding out
2992 ELSE
2993 icontact(i) = 0
2994 istep(i) = 0
2995 cycle
2996 ENDIF
2997C----- internal sliding
2998 ELSEIF (bbb<zero) THEN
2999 nni =zero
3000 IF(itq ==1.AND.mvoisin(4,l)==0)THEN
3001 nqx = sx4150(i)
3002 nqy = sy4150(i)
3003 nqz = sz4150(i)
3004 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3005 . (nqx*ybd - nqy*xbd)*zib
3006 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 4
3007 ELSEIF(itq ==2.AND.mvoisin(1,l)==0)THEN
3008 nqx = sx1250(i)
3009 nqy = sy1250(i)
3010 nqz = sz1250(i)
3011 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3012 . (nqx*ybd - nqy*xbd)*zib
3013 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 1
3014 ELSEIF(itq ==3.AND.mvoisin(2,l)==0)THEN
3015 nqx = sx3450(i)
3016 nqy = sy3450(i)
3017 nqz = sz3450(i)
3018 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3019 . (nqx*ybd - nqy*xbd)*zib
3020 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 2
3021 ELSEIF(itq ==4.AND.mvoisin(3,l)==0)THEN
3022 nqx = sx2350(i)
3023 nqy = sy2350(i)
3024 nqz = sz2350(i)
3025 nni = (nqy*zbd - nqz*ybd)*xib+(nqz*xbd - nqx*zbd)*yib+
3026 . (nqx*ybd - nqy*xbd)*zib
3027 IF (nni<zero.AND.bbb<aaa) subtria_n(i) = 3
3028 ENDIF
3029 IF (nni>zero) THEN
3030 icontact(i) = 0
3031 istep(i) = 0
3032 cycle
3033 END IF
3034 END IF
3035C-
3036 xic = xc(i)-xi(i)
3037 yic = yc(i)-yi(i)
3038 zic = zc(i)-zi(i)
3039 aaa = xic*xbc+yic*ybc+zic*zbc
3040 bbb = xic*xce+yic*yce+zic*zce
3041C----- sliding around C
3042 IF (aaa<zero.AND.subtria_n(i)==0) THEN
3043 IF(itq ==1.AND.mvoisin(2,l)>0)THEN
3044 subtria_n(i) = 10
3045 ELSEIF(itq ==2.AND.mvoisin(3,l)>0)THEN
3046 subtria_n(i) = 11
3047 ELSEIF(itq ==3.AND.mvoisin(4,l)>0)THEN
3048 subtria_n(i) = 12
3049 ELSEIF(itq ==4.AND.mvoisin(1,l)>0)THEN
3050 subtria_n(i) = 9
3051C----- sliding out
3052 ELSE
3053 icontact(i) = 0
3054 istep(i) = 0
3055 cycle
3056 ENDIF
3057C----- internal sliding
3058 ELSEIF (bbb>zero) THEN
3059 nni =zero
3060 IF(itq ==1.AND.mvoisin(2,l)==0)THEN
3061 nqx = sx2350(i)
3062 nqy = sy2350(i)
3063 nqz = sz2350(i)
3064 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3065 . (nqx*yce - nqy*xce)*zic
3066 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 2
3067 ELSEIF(itq ==2.AND.mvoisin(1,l)==0)THEN
3068 nqx = sx3450(i)
3069 nqy = sy3450(i)
3070 nqz = sz3450(i)
3071 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3072 . (nqx*yce - nqy*xce)*zic
3073 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 3
3074 ELSEIF(itq ==3.AND.mvoisin(2,l)==0)THEN
3075 nqx = sx4150(i)
3076 nqy = sy4150(i)
3077 nqz = sz4150(i)
3078 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3079 . (nqx*yce - nqy*xce)*zic
3080 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 4
3081 ELSEIF(itq ==4.AND.mvoisin(3,l)==0)THEN
3082 nqx = sx1250(i)
3083 nqy = sy1250(i)
3084 nqz = sz1250(i)
3085 nni = (nqy*zce - nqz*yce)*xic+(nqz*xce - nqx*zce)*yic+
3086 . (nqx*yce - nqy*xce)*zic
3087 IF (nni<zero.AND.bbb>aaa) subtria_n(i) = 1
3088 ENDIF
3089 IF (nni>zero) THEN
3090 icontact(i) = 0
3091 istep(i) = 0
3092 cycle
3093 END IF
3094 END IF
3095 IF (subtria_n(i)>0) subtria(i) = subtria_n(i)
3096 subtria_n(i) = itq
3097 END DO
3098C------possible nodal normal with interpolation after
3099 DO i=1,jlt
3100 IF(istep(i) /= 6 )cycle
3101C---- finalizing penetration for node/edge possible multi-contactcase for new contact
3102 itq = subtria(i)
3103 IF (itq>20) itq = itq -20
3104C EDGE : BC
3105 k = ic(1,itq) ! 5
3106 xa(i) = xx0(i,k)
3107 ya(i) = yy0(i,k)
3108 za(i) = zz0(i,k)
3109
3110 k = ic(2,itq) ! N1
3111 xb(i) = xx0(i,k)
3112 yb(i) = yy0(i,k)
3113 zb(i) = zz0(i,k)
3114
3115 k = ic(3,itq) ! N2
3116 xc(i) = xx0(i,k)
3117 yc(i) = yy0(i,k)
3118 zc(i) = zz0(i,k)
3119C--------
3120 xbc = xc(i)-xb(i)
3121 ybc = yc(i)-yb(i)
3122 zbc = zc(i)-zb(i)
3123C--- edge normal
3124 IF(itq ==1)THEN
3125 nqx = sx1250(i)
3126 nqy = sy1250(i)
3127 nqz = sz1250(i)
3128 ELSEIF(itq ==2)THEN
3129 nqx = sx2350(i)
3130 nqy = sy2350(i)
3131 nqz = sz2350(i)
3132 ELSEIF(itq ==3)THEN
3133 nqx = sx3450(i)
3134 nqy = sy3450(i)
3135 nqz = sz3450(i)
3136 ELSEIF(itq ==4)THEN
3137 nqx = sx4150(i)
3138 nqy = sy4150(i)
3139 nqz = sz4150(i)
3140 ELSE
3141C--- re-compute seg normal for sliding case
3142 xab = xb(i)-xa(i)
3143 yab = yb(i)-ya(i)
3144 zab = zb(i)-za(i)
3145 nqx = yab*zbc - zab*ybc
3146 nqy = zab*xbc - xab*zbc
3147 nqz = xab*ybc - yab*xbc
3148 ENDIF
3149 nx(i) = nqy*zbc - nqz*ybc
3150 ny(i) = nqz*xbc - nqx*zbc
3151 nz(i) = nqx*ybc - nqy*xbc
3152C in right direction
3153 nqx = -nx(i)
3154 nqy = -ny(i)
3155 nqz = -nz(i)
3156
3157 aaa = max(em30,sqrt(nqx*nqx + nqy*nqy + nqz*nqz))
3158 s2 = one/aaa
3159 nqx = nqx*s2
3160 nqy = nqy*s2
3161 nqz = nqz*s2
3162
3163 bbb = (xi(i)-xb(i))*nqx+(yi(i)-yb(i))*nqy+(zi(i)-zb(i))*nqz
3164 pene(i) = -min(zero,bbb)
3165 IF(pene(i) == zero)THEN
3166 icontact(i) = 0
3167 cycle
3168 ENDIF
3169 n1(i) = nqx
3170 n2(i) = nqy
3171 n3(i) = nqz
3172 ix1(i) = ixx(i,ic( 7,itq))
3173 ix2(i) = ixx(i,ic( 8,itq))
3174 ix3(i) = ixx(i,ic( 9,itq))
3175 ix4(i) = ixx(i,ic(10,itq))
3176
3177 x1(i) = xx0(i,ic( 7,itq))
3178 x2(i) = xx0(i,ic( 8,itq))
3179 x3(i) = xx0(i,ic( 9,itq))
3180 x4(i) = xx0(i,ic(10,itq))
3181 y1(i) = yy0(i,ic( 7,itq))
3182 y2(i) = yy0(i,ic( 8,itq))
3183 y3(i) = yy0(i,ic( 9,itq))
3184 y4(i) = yy0(i,ic(10,itq))
3185 z1(i) = zz0(i,ic( 7,itq))
3186 z2(i) = zz0(i,ic( 8,itq))
3187 z3(i) = zz0(i,ic( 9,itq))
3188 z4(i) = zz0(i,ic(10,itq))
3189
3190 vx1(i) = vx(i,ic( 7,itq))
3191 vx2(i) = vx(i,ic( 8,itq))
3192 vx3(i) = vx(i,ic( 9,itq))
3193 vx4(i) = vx(i,ic(10,itq))
3194 vy1(i) = vy(i,ic( 7,itq))
3195 vy2(i) = vy(i,ic( 8,itq))
3196 vy3(i) = vy(i,ic( 9,itq))
3197 vy4(i) = vy(i,ic(10,itq))
3198 vz1(i) = vz(i,ic( 7,itq))
3199 vz2(i) = vz(i,ic( 8,itq))
3200 vz3(i) = vz(i,ic( 9,itq))
3201 vz4(i) = vz(i,ic(10,itq))
3202 h1(i) = zero
3203 h2(i) = zero
3204 h3(i) = half
3205 h4(i) = half
3206
3207 n=cand_n(i)
3208 mglob = mseglo(cand_e(i))
3209 IF (itq==5 .or. itq==9 .or. itq==13 .or. itq==17)THEN
3210 mglob = mvoisin(1,cand_e(i))
3211c ITQ = ITRIV(1,I)
3212 ELSEIF(itq==6 .or. itq==10 .or. itq==14 .or. itq==18)THEN
3213 mglob = mvoisin(2,cand_e(i))
3214c ITQ = ITRIV(2,I)
3215 ELSEIF(itq==7 .or. itq==11 .or. itq==15 .or. itq==19)THEN
3216 mglob = mvoisin(3,cand_e(i))
3217c ITQ = ITRIV(3,I)
3218 ELSEIF(itq==8 .or. itq==12 .or. itq==16 .or. itq==20)THEN
3219 mglob = mvoisin(4,cand_e(i))
3220c ITQ = ITRIV(4,I)
3221 ENDIF
3222 IF (itq>8) itq = subtria_n(i)
3223
3224 IF(n <= nsn)THEN
3225 IF(irtlm(1,n)>0) THEN
3226 irtlm(1,n) = mglob
3227 time_s(n) = -ep20
3228 irtlm(2,n)= itq+20 !sub triangle
3229 ELSEIF((time_s(n) < pene(i)) .or. ! multi or new intersection
3230 . (time_s(n) == pene(i) .and. ! multi at same time
3231 . -irtlm(1,n) < mglob ) )THEN
3232 irtlm(2,n)= itq+20 !sub triangle
3233 irtlm(1,n) = -mglob
3234 time_s(n) = pene(i)
3235 ENDIF
3236 ELSE ! Remote nodes in SPMD
3237 ii=n-nsn
3238 IF(irtlm_fi(nin)%P(1,ii) > 0)THEN
3239 irtlm_fi(nin)%P(1,ii) = mglob
3240 time_sfi(nin)%P(ii) = -ep20
3241 irtlm_fi(nin)%P(2,ii) = itq+20
3242 ELSEIF((time_sfi(nin)%P(ii) < pene(i)).or.
3243 . (time_sfi(nin)%P(ii)== pene(i).and.
3244 . -irtlm_fi(nin)%P(1,ii) < mglob) )THEN
3245 irtlm_fi(nin)%P(2,ii) = itq+20
3246 irtlm_fi(nin)%P(1,ii) = -mglob
3247 time_sfi(nin)%P(ii) = pene(i)
3248 ENDIF
3249 ENDIF
3250 END DO
3251c write(iout,*)'i24dst3 9'
3252C=======================================================================
3253c 6 Reset IRTLM (only for old contact: IRTLM(1,N)>0)
3254c and define PENE_OLD(3,N) = PENE-GAP or dist form SECONDARY to surf
3255C=======================================================================
3256C IF (INCONV==1) THEN
3257 DO i=1,jlt
3258 ce_loc(i) = cand_e(i)
3259 cn_loc(i) = cand_n(i)
3260 IF(pene(i) == 0 )THEN
3261
3262 n=cand_n(i)
3263 IF(n <= nsn)THEN
3264c IF(IRTLM(1,N)==MSEGLO(CAND_E(I)).and.
3265c . IRTLM(1,N)==IRTLM_OLD(I))THEN
3266 IF(irtlm(1,n) > 0)THEN
3267 irtlm(1,n)=0
3268 time_s(n) = -ep20
3269 secnd_fr(4,n)= zero
3270 secnd_fr(5,n)= zero
3271 secnd_fr(6,n)= zero
3272 IF (imconv==1) pene_old(2,n) = zero
3273 pene_old(5,n) = zero
3274 IF (imconv==1) stif_old(2,n) = zero
3275 ENDIF
3276 ELSE
3277 IF(irtlm_fi(nin)%P(1,n-nsn) > 0)THEN
3278 irtlm_fi(nin)%P(1,n-nsn)=0
3279 time_sfi(nin)%P(n-nsn) = -ep20
3280 secnd_frfi(nin)%P(4,n-nsn)= zero
3281 secnd_frfi(nin)%P(5,n-nsn)= zero
3282 secnd_frfi(nin)%P(6,n-nsn)= zero
3283 IF (imconv==1) pene_oldfi(nin)%P(2,n-nsn) = zero
3284 pene_oldfi(nin)%P(5,n-nsn) = zero
3285 IF (imconv==1) stif_oldfi(nin)%P(2,n-nsn) = zero
3286 ENDIF
3287 ENDIF
3288 ELSE
3289 ns=nsvg(i)
3290 IF (nm1(i)==zero.AND.nm2(i)==zero.AND.nm3(i)==zero) THEN
3291 nm1(i) = n1(i)
3292 nm2(i) = n2(i)
3293 nm3(i) = n3(i)
3294 END IF
3295C------case huge warped element
3296 IF (ix1(i) == ix2(i))THEN
3297 ELSEIF(ix2(i) == ix3(i))THEN
3298c ELSEIF(IX3(I) == IX4(I))THEN impossible
3299 ELSEIF(ix4(i) == ix1(i))THEN
3300 ELSE
3301C -------Calculate Nz,Area,Z1 IF Z1*Z1/AREA > f(angle 45) CYCLE
3302 ENDIF
3303 xs = h1(i)*x1(i) + h2(i)*x2(i) + h3(i)*x3(i) + h4(i)*x4(i)
3304 ys = h1(i)*y1(i) + h2(i)*y2(i) + h3(i)*y3(i) + h4(i)*y4(i)
3305 zs = h1(i)*z1(i) + h2(i)*z2(i) + h3(i)*z3(i) + h4(i)*z4(i)
3306 xs = xs - xi(i)
3307 ys = ys - yi(i)
3308 zs = zs - zi(i)
3309
3310 aaa = xs*n1(i) + ys*n2(i) + zs*n3(i)
3311
3312 IF(aaa > zero)THEN
3313
3314 aaa = xs**2 + ys**2 + zs**2
3315 aaa = onep01*sqrt(aaa)
3316 pmax_gap = max(pmax_gap,aaa)
3317
3318 n=cand_n(i)
3319 IF(n <= nsn)THEN
3320 pene_old(3,n) = max(pene_old(3,n),aaa)
3321 ELSE
3322 pene_oldfi(nin)%P(3,n-nsn) =
3323 . max(pene_oldfi(nin)%P(3,n-nsn),aaa)
3324 ENDIF
3325 ENDIF
3326 ENDIF
3327 ENDDO
3328#include "lockoff.inc"
3329C END IF !(INCONV==1) THEN
3330C-----Correction due to difference between starter/engine
3331 IF (inacti==5) THEN
3332 IF (tt==zero) THEN
3333 DO i=1,jlt
3334 IF(pene(i) == zero)cycle
3335 jg = nsvg(i)
3336 n = cand_n(i)
3337 IF(jg > 0)THEN
3338#include "lockon.inc"
3339 pene_old(5,n) = max( pene(i) ,pene_old(5,n) )
3340 stif_old(1,n) = max( stif(i) ,stif_old(1,n) )
3341#include "lockoff.inc"
3342 ELSE
3343 jg = -jg
3344#include "lockon.inc"
3345 pene_oldfi(nin)%P(5,jg) = max( pene(i),pene_oldfi(nin)%P(5,jg))
3346 stif_oldfi(nin)%P(1,jg) = max( stif(i),stif_oldfi(nin)%P(1,jg))
3347#include "lockoff.inc"
3348 ENDIF
3349 ENDDO
3350C-------zero force for all first contact
3351 ELSE
3352 DO i=1,jlt
3353 IF(pene(i) == zero)cycle
3354 jg = nsvg(i)
3355 n = cand_n(i)
3356 IF(jg > 0)THEN
3357 IF (ipen0(i)==0)THEN
3358#include "lockon.inc"
3359 pene_old(5,n) = max( pene(i) ,pene_old(5,n) )
3360 stif_old(1,n) = max( stif(i) ,stif_old(1,n) )
3361#include "lockoff.inc"
3362 END IF
3363 ELSE
3364 jg = -jg
3365 IF (ipen0(i)==0)THEN
3366#include "lockon.inc"
3367 pene_oldfi(nin)%P(5,jg) = max( pene(i),pene_oldfi(nin)%P(5,jg))
3368 stif_oldfi(nin)%P(1,jg) = max( stif(i),stif_oldfi(nin)%P(1,jg))
3369#include "lockoff.inc"
3370 END IF
3371 ENDIF
3372 ENDDO
3373 END IF !(TT==ZERO) THEN
3374C---------PENE -> PENE_Offset
3375 DO i=1,jlt
3376 IF(pene(i) == zero)cycle
3377 jg = nsvg(i)
3378 n = cand_n(i)
3379 IF(jg > 0)THEN
3380 pene_sh = max(zero,pene(i)-pene_old(5,n))
3381 ELSE
3382 jg = -jg
3383 pene_sh = max(zero,pene(i)-pene_oldfi(nin)%P(5,jg))
3384 ENDIF
3385C---------PENE -> updated bt PENE_Offset
3386 pene(i) = pene_sh
3387 ENDDO
3388 END IF !IF (INACTI==5)
3389C-----peneref for nonlinear stif (quadratic...)
3390 penref(1:jlt) = zero
3391 IF (inacti==-1) THEN
3392 DO i=1,jlt
3393 IF(pene(i) == zero.OR.iknon(i)==0) cycle
3394 penref(i) = sqrt(eps2*area(i))
3395 jg = nsvg(i)
3396 n = cand_n(i)
3397 IF(jg > 0)THEN
3398 pene_sh = em01*pene_old(5,n)
3399 ELSE
3400 jg = -jg
3401 pene_sh = em01*pene_oldfi(nin)%P(5,jg)
3402 ENDIF
3403 penref(i) = max(penref(i),pene_sh)
3404 IF (iknon(i)<0) THEN
3405 f_pene = em01*pene(i)/pene_sh
3406 penref(i) = zero ! not to do quadratic stif
3407 IF (f_pene>two_third) THEN
3408 fac_k = min(twenty,six*f_pene)
3409 penref(i) = pene_sh/sqrt(fac_k)
3410 END IF
3411 END IF
3412 ENDDO
3413 ELSE
3414 DO i=1,jlt
3415 IF(pene(i) == zero.OR.iknon(i)==0) cycle
3416 n=cand_n(i)
3417 l=cand_e(i)
3418 penref(i) = one_fifth*sqrt(area(i))
3419 gap2 = one_fifth*(gaps(i)+gap_m(l))
3420 IF (gap2 > em04) penref(i)=min(penref(i),gap2)
3421 jg = nsvg(i)
3422 IF(jg > 0)THEN
3423 pene_sh = pene_old(5,n)
3424 ELSE
3425 jg = -jg
3426 pene_sh = pene_oldfi(nin)%P(5,jg)
3427 ENDIF
3428 penref(i) = max(penref(i),pene_sh)
3429 IF(iknon(i)==1) penref(i) = ten*penref(i)
3430 IF(iknon(i)>2) penref(i) = one_fifth*penref(i)
3431 ENDDO
3432 END IF
3433C
3434 RETURN
3435 END
3436!||====================================================================
3437!|| i24nexttria ../engine/source/interfaces/int24/i24dst3.F
3438!||--- called by ------------------------------------------------------
3439!|| i24dst3 ../engine/source/interfaces/int24/i24dst3.F
3440!||--- calls -----------------------------------------------------
3441!|| i24nexttria2 ../engine/source/interfaces/int24/i24dst3.F
3442!||--- uses -----------------------------------------------------
3443!|| tri7box ../engine/share/modules/tri7box.F
3444!||====================================================================
3445 SUBROUTINE i24nexttria(
3446 1 IZLIM,ISTEP,SUBTRIA_N,SUBTRIA,
3447 2 LARGEP,PENE ,XXL ,YYL ,ZZL ,
3448 3 SX125,SY125,SZ125,SX235,SY235,
3449 4 SZ235,SX345,SY345,SZ345,SX415,
3450 5 SY415,SZ415,XI ,YI ,ZI ,
3451 6 N1 ,N2 ,N3 ,LA ,LB ,
3452 7 LC ,GAP2 ,EPS0 )
3453C-----------------------------------------------
3454C M o d u l e s
3455C-----------------------------------------------
3456 USE tri7box
3457C-----------------------------------------------
3458C I m p l i c i t T y p e s
3459C-----------------------------------------------
3460#include "implicit_f.inc"
3461C-----------------------------------------------
3462C G l o b a l P a r a m e t e r s
3463C-----------------------------------------------
3464#include "comlock.inc"
3465C-----------------------------------------------
3466C D u m m y A r g u m e n t s
3467C-----------------------------------------------
3468 INTEGER ISTEP,SUBTRIA_N,SUBTRIA,IZLIM,LARGEP
3469 my_real
3470 . PENE,XXL(17),YYL(17),ZZL(17),XI ,YI ,ZI,
3471 . sx125,sy125,sz125,sx235,sy235,
3472 . sz235,sx345,sy345,sz345,sx415,
3473 . sy415,sz415,n1,n2,n3,la,lb,lc,gap2,eps0
3474C-----------------------------------------------
3475C L o c a l V a r i a b l e s
3476C-----------------------------------------------
3477 INTEGER I, J, K,ITQ
3478 my_real
3479 . AAA,BBB,VOL,XA,YA,ZA,XB,YB,ZB,XC,YC,ZC,XD,YD,ZD,
3480 . XAB,YAB,ZAB,XBC,YBC,ZBC,XCA,YCA,ZCA,XBD,YBD,ZBD,
3481 . XIA,YIA,ZIA,XIB,YIB,ZIB,XIC,YIC,ZIC,SX ,SY, SZ,
3482 . s2 ,sax,say,saz,sbx,sby,sbz,unp,zerom,eps
3483 INTEGER IT0(3,20),IC(10,20)
3484! +----+----------+----------+-------------+
3485! | IC |
3486! +----+----------+----------+-------------+
3487! | ST | noeuds T | noeuds TV| noeuds Quad |
3488! +----+----------+----------+-------------+ 11-------10
3489! | 1 | 5 1 2 | 14 3 4 | 3 4 1 2 | |\ 19 /|
3490! | 2 | 5 2 3 | 15 4 1 | 4 1 2 3 | | \ / |
3491! | 3 | 5 3 4 | 16 1 2 | 1 2 3 4 | | \ / |
3492! | 4 | 5 4 1 | 17 2 3 | 2 3 4 1 | | 16 |
3493! +----+----------+----------+-------------+ |15/ \11|
3494! | 5 | 14 2 1 | 5 6 7 | 6 7 2 1 | | / \ |
3495! | 6 | 15 3 2 | 5 8 9 | 8 9 3 2 | |/ 7 \|
3496! | 7 | 16 4 3 | 5 10 11 | 10 11 4 3 |12-------4-------3-------9
3497! | 8 | 17 1 4 | 5 12 13 | 12 13 1 4 | |\ 12 /|\ /|\ 14 /|
3498! +----+----------+----------+-------------+ | \ / | \ 3 / | \ / |
3499! | 9 | 14 1 6 |(13) 7 2 | 7 2 1 6 | | \ / | \ /2 |6 \ /18|
3500! | 10 | 15 2 8 |( 7) 9 3 | 9 3 2 8 | | 17 | 5 | 15 |
3501! | 11 | 16 3 10 |( 9)11 4 | 11 4 3 10 | |20/ \ 8| 4/ \ | / \ |
3502! | 12 | 17 4 12 |(11)13 1 | 13 1 4 12 | | / \ | / 1 \ | / \ |
3503! +----+----------+----------+-------------+ |/ 16 \|/ \|/ 10 \|
3504! | 13 | 14 7 2 |( 8) 1 6 | 1 6 7 2 |13-------1-------2-------8
3505! | 14 | 15 9 3 |(10) 2 8 | 2 8 9 3 | |\ 5 /|
3506! | 15 | 16 11 4 |(12) 3 10 | 3 10 11 4 | | \ / |
3507! | 16 | 17 13 1 |( 6) 4 12 | 4 12 13 1 | |9 \ /13|
3508! +----+----------+----------+-------------+ | 14 |
3509! | 17 | 14 6 7 | - 2 1 | 2 1 6 7 | | / \ |
3510! | 18 | 15 8 9 | - 3 2 | 3 2 8 9 | | / \ |
3511! | 19 | 16 10 11 | - 4 3 | 4 3 10 11 | |/ 17 \|
3512! | 20 | 17 12 13 | - 1 4 | 1 4 12 13 | 6-------7
3513! +----+----------+----------+-------------+
3514 DATA ic /
3515 1 5, 1, 2 , 14, 3, 4 , 3, 4, 1, 2,
3516 2 5, 2, 3 , 15, 4, 1 , 4, 1, 2, 3,
3517 3 5, 3, 4 , 16, 1, 2 , 1, 2, 3, 4,
3518 4 5, 4, 1 , 17, 2, 3 , 2, 3, 4, 1,
3519 5 14, 2, 1 , 5, 6, 7 , 6, 7, 2, 1,
3520 6 15, 3, 2 , 5, 8, 9 , 8, 9, 3, 2,
3521 7 16, 4, 3 , 5,10,11 , 10,11, 4, 3,
3522 8 17, 1, 4 , 5,12,13 , 12,13, 1, 4,
3523 9 14, 1, 6 , 0, 7, 2 , 7, 2, 1, 6,
3524 . 15, 2, 8 , 0, 9, 3 , 9, 3, 2, 8,
3525 1 16, 3,10 , 0,11, 4 , 11, 4, 3,10,
3526 2 17, 4,12 , 0,13, 1 , 13, 1, 4,12,
3527 3 14, 7, 2 , 0, 1, 6 , 1, 6, 7, 2,
3528 4 15, 9, 3 , 0, 2, 8 , 2, 8, 9, 3,
3529 5 16,11, 4 , 0, 3,10 , 3,10,11, 4,
3530 6 17,13, 1 , 0, 4,12 , 4,12,13, 1,
3531 7 14, 6, 7 , 0, 2, 1 , 2, 1, 6, 7,
3532 8 15, 8, 9 , 0, 3, 2 , 3, 2, 8, 9,
3533 9 16,10,11 , 0, 4, 3 , 4, 3,10,11,
3534 . 17,12,13 , 0, 1, 4 , 1, 4,12,13/
3535
3536! +------+-------------------------------------------+
3537! | | sous-triangle voisins |
3538! |sous +----------+----------+---------------------+
3539! |trian | n3/=n4 | n3=n4 |
3540! | +----------+----------+----------+----------+
3541! | | IT0 | 2,4,6,8=0| | 2,4,6,8=0|
3542! +------+----------+----------+----------+----------+
3543! | 1 | 5 2 4 | 5 2 4 | 5 6 8 | 5 6 8 |
3544! | 2 | 6 3 1 | 6 3 1 | - - - | - - - |
3545! | 3 | 7 4 2 | 7 4 2 | - - - | - - - |
3546! | 4 | 8 1 3 | 8 1 3 | - - - | - - - |
3547! +------+----------+----------+----------+----------+
3548! | 5 | 1 9 13 | 1 -1 -1 | 1 9 13 | 1 -1 -1 |
3549! | 6 | 2 10 14 | 2 -1 -1 | 1 10 14 | 1 -1 -1 |
3550! | 7 | 3 11 15 | 3 -1 -1 | - - - | - - - |
3551! | 8 | 4 12 16 | 4 -1 -1 | 1 12 16 | 1 -1 -1 |
3552! +------+----------+----------+----------+----------+
3553! | 9 | -1 17 5 | - - - | -1 17 5 | - - - |
3554! | 10 | -1 18 6 | - - - | -1 18 6 | - - - |
3555! | 11 | -1 19 7 | - - - | - - - | - - - |
3556! | 12 | -1 20 8 | - - - | -1 20 8 | - - - |
3557! +------+----------+----------+----------+----------+
3558! | 13 | -1 5 17 | - - - | -1 5 17 | - - - |
3559! | 14 | -1 6 18 | - - - | -1 6 18 | - - - |
3560! | 15 | -1 7 19 | - - - | - - - | - - - |
3561! | 16 | -1 8 20 | - - - | -1 8 20 | - - - |
3562! +------+----------+----------+----------+----------+
3563! | 17 | -1 9 13 | - - - | -1 9 13 | - - - |
3564! | 18 | -1 10 14 | - - - | -1 10 14 | - - - |
3565! | 19 | -1 11 15 | - - - | - - - | - - - |
3566! | 20 | -1 12 16 | - - - | -1 12 16 | - - - |
3567! +------+----------+----------+----------+----------+
3568 DATA it0 /
3569 1 5, 2, 4,
3570 2 6, 3, 1,
3571 3 7, 4, 2,
3572 4 8, 1, 3,
3573 5 1, 9,13,
3574 6 2,10,14,
3575 7 3,11,15,
3576 8 4,12,16,
3577 9 -1,17, 5,
3578 . -1,18, 6,
3579 1 -1,19, 7,
3580 2 -1,20, 8,
3581 3 -1, 5,17,
3582 4 -1, 6,18,
3583 5 -1, 7,19,
3584 6 -1, 8,20,
3585 7 -1, 9,13,
3586 8 -1,10,14,
3587 9 -1,11,15,
3588 . -1,12,16/
3589
3590 unp = one + em4
3591 zerom = zero - em4
3592 eps = eps0
3593
3594 itq = subtria
3595 k = ic(1,itq)
3596 xa = xxl(k)
3597 ya = yyl(k)
3598 za = zzl(k)
3599
3600 k = ic(2,itq)
3601 xb = xxl(k)
3602 yb = yyl(k)
3603 zb = zzl(k)
3604
3605 k = ic(3,itq)
3606 xc = xxl(k)
3607 yc = yyl(k)
3608 zc = zzl(k)
3609
3610 k = ic(4,itq)
3611 xd = xxl(k)
3612 yd = yyl(k)
3613 zd = zzl(k)
3614
3615 IF(itq ==1)THEN
3616 n1 = sx125
3617 n2 = sy125
3618 n3 = sz125
3619 ELSEIF(itq ==2)THEN
3620 n1 = sx235
3621 n2 = sy235
3622 n3 = sz235
3623 ELSEIF(itq ==3)THEN
3624 n1 = sx345
3625 n2 = sy345
3626 n3 = sz345
3627 ELSEIF(itq ==4)THEN
3628 n1 = sx415
3629 n2 = sy415
3630 n3 = sz415
3631 ENDIF
3632
3633 s2 = one/max(em30,sqrt(n1*n1 + n2*n2 + n3*n3))
3634 n1 = n1*s2
3635 n2 = n2*s2
3636 n3 = n3*s2
3637
3638 bbb = (xi-xa)*n1+(yi-ya)*n2+(zi-za)*n3
3639 pene = -min(zero,bbb)
3640
3641 xab = xb-xa
3642 yab = yb-ya
3643 zab = zb-za
3644 xbc = xc-xb
3645 ybc = yc-yb
3646 zbc = zc-zb
3647 xca = xa-xc
3648 yca = ya-yc
3649 zca = za-zc
3650
3651 xbd = xd-xb
3652 ybd = yd-yb
3653 zbd = zd-zb
3654
3655 xia = xa-xi
3656 yia = ya-yi
3657 zia = za-zi
3658 xib = xb-xi
3659 yib = yb-yi
3660 zib = zb-zi
3661 xic = xc-xi
3662 yic = yc-yi
3663 zic = zc-zi
3664 sx = - yab*zca + zab*yca
3665 sy = - zab*xca + xab*zca
3666 sz = - xab*yca + yab*xca
3667 s2 = sx*sx+sy*sy+sz*sz
3668 sax = yib*zic - zib*yic
3669 say = zib*xic - xib*zic
3670 saz = xib*yic - yib*xic
3671 la = (sx*sax+sy*say+sz*saz)/s2
3672 sbx = yic*zia - zic*yia
3673 sby = zic*xia - xic*zia
3674 sbz = xic*yia - yic*xia
3675 lb = (sx*sbx+sy*sby+sz*sbz)/s2
3676 lc = one - la - lb
3677
3678c------------------------------------------------------
3679c quadrangles
3680c check if outside side CA
3681 aaa = zerom
3682 IF(lb<aaa)THEN
3683 istep=3
3684 subtria_n=it0(2,itq)
3685 subtria =it0(2,itq)
3686 izlim=-1
3687 CALL i24nexttria2(
3688 1 izlim,istep,subtria_n,subtria,
3689 2 largep,pene ,xxl ,yyl ,zzl ,
3690 3 sx125,sy125,sz125,sx235,sy235,
3691 4 sz235,sx345,sy345,sz345,sx415,
3692 5 sy415,sz415,xi ,yi ,zi ,
3693 6 n1 ,n2 ,n3 ,la ,lb ,
3694 7 lc ,gap2 ,eps )
3695c check if outside side AB
3696 ELSEIF(lc<aaa)THEN
3697 istep=3
3698 subtria_n=it0(3,itq)
3699 subtria =it0(3,itq)
3700 izlim=-1
3701 CALL i24nexttria2(
3702 1 izlim,istep,subtria_n,subtria,
3703 2 largep,pene ,xxl ,yyl ,zzl ,
3704 3 sx125,sy125,sz125,sx235,sy235,
3705 4 sz235,sx345,sy345,sz345,sx415,
3706 5 sy415,sz415,xi ,yi ,zi ,
3707 6 n1 ,n2 ,n3 ,la ,lb ,
3708 7 lc ,gap2 ,eps )
3709c check if outside side BC
3710 ELSEIF(la<-eps)THEN
3711c hors zone d'extension de surface
3712 istep=3
3713 subtria_n=it0(1,itq)
3714 izlim=-1
3715 ELSEIF(la<eps)THEN
3716c zone limite interpolation des normales si convex
3717 vol = n1*xbd + n2*ybd + n3*zbd
3718 IF (gap2<eps.AND.(lb<eps.OR.lc<eps)) vol=-abs(vol)
3719 IF(largep == 1)THEN
3720c large penetration inside +- extension zone
3721 istep=3
3722 subtria_n=it0(1,itq)
3723 izlim = -1
3724 ELSEIF(vol < zero)THEN
3725c convex
3726 izlim=1
3727 ELSEIF(la<zerom)THEN
3728 istep=3
3729 subtria_n=it0(1,itq)
3730 izlim=-1
3731 ENDIF
3732 ENDIF
3733
3734 RETURN
3735 END
3736!||====================================================================
3737!|| i24nexttria2 ../engine/source/interfaces/int24/i24dst3.F
3738!||--- called by ------------------------------------------------------
3739!|| i24nexttria ../engine/source/interfaces/int24/i24dst3.F
3740!||--- uses -----------------------------------------------------
3741!|| tri7box ../engine/share/modules/tri7box.F
3742!||====================================================================
3743 SUBROUTINE i24nexttria2(
3744 1 IZLIM,ISTEP,SUBTRIA_N,SUBTRIA,
3745 2 LARGEP,PENE ,XXL ,YYL ,ZZL ,
3746 3 SX125,SY125,SZ125,SX235,SY235,
3747 4 SZ235,SX345,SY345,SZ345,SX415,
3748 5 SY415,SZ415,XI ,YI ,ZI ,
3749 6 N1 ,N2 ,N3 ,LA ,LB ,
3750 7 LC ,GAP2 ,EPS )
3751C-----------------------------------------------
3752C*******************************************
3753C
3754C THIS SUBROUTINE IS EXACTLY THE SAME
3755C THAN I24NEXTTRIA2 EXCEPT THE CALL TO I24NEXTTRIA2
3756C
3757C*******************************************
3758C-----------------------------------------------
3759C M o d u l e s
3760C-----------------------------------------------
3761 USE tri7box
3762C-----------------------------------------------
3763C I m p l i c i t T y p e s
3764C-----------------------------------------------
3765#include "implicit_f.inc"
3766C-----------------------------------------------
3767C G l o b a l P a r a m e t e r s
3768C-----------------------------------------------
3769#include "comlock.inc"
3770C-----------------------------------------------
3771C D u m m y A r g u m e n t s
3772C-----------------------------------------------
3773 INTEGER ISTEP,SUBTRIA_N,SUBTRIA,IZLIM,LARGEP
3774 my_real
3775 . PENE,XXL(17),YYL(17),ZZL(17),XI ,YI ,ZI,
3776 . sx125,sy125,sz125,sx235,sy235,
3777 . sz235,sx345,sy345,sz345,sx415,
3778 . sy415,sz415,n1,n2,n3,la,lb,lc,gap2,eps
3779C-----------------------------------------------
3780C L o c a l V a r i a b l e s
3781C-----------------------------------------------
3782 INTEGER I, J, K,ITQ
3783 my_real
3784 . AAA,BBB,VOL,XA,YA,ZA,XB,YB,ZB,XC,YC,ZC,XD,YD,ZD,
3785 . XAB,YAB,ZAB,XBC,YBC,ZBC,XCA,YCA,ZCA,XBD,YBD,ZBD,
3786 . XIA,YIA,ZIA,XIB,YIB,ZIB,XIC,YIC,ZIC,SX ,SY, SZ,
3787 . S2 ,SAX,SAY,SAZ,SBX,SBY,SBZ,UNP,ZEROM
3788 INTEGER IT0(3,20),IC(10,20)
3789! +----+----------+----------+-------------+
3790! | IC |
3791! +----+----------+----------+-------------+
3792! | ST | noeuds T | noeuds TV| noeuds Quad |
3793! +----+----------+----------+-------------+ 11-------10
3794! | 1 | 5 1 2 | 14 3 4 | 3 4 1 2 | |\ 19 /|
3795! | 2 | 5 2 3 | 15 4 1 | 4 1 2 3 | | \ / |
3796! | 3 | 5 3 4 | 16 1 2 | 1 2 3 4 | | \ / |
3797! | 4 | 5 4 1 | 17 2 3 | 2 3 4 1 | | 16 |
3798! +----+----------+----------+-------------+ |15/ \11|
3799! | 5 | 14 2 1 | 5 6 7 | 6 7 2 1 | | / \ |
3800! | 6 | 15 3 2 | 5 8 9 | 8 9 3 2 | |/ 7 \|
3801! | 7 | 16 4 3 | 5 10 11 | 10 11 4 3 |12-------4-------3-------9
3802! | 8 | 17 1 4 | 5 12 13 | 12 13 1 4 | |\ 12 /|\ /|\ 14 /|
3803! +----+----------+----------+-------------+ | \ / | \ 3 / | \ / |
3804! | 9 | 14 1 6 |(13) 7 2 | 7 2 1 6 | | \ / | \ /2 |6 \ /18|
3805! | 10 | 15 2 8 |( 7) 9 3 | 9 3 2 8 | | 17 | 5 | 15 |
3806! | 11 | 16 3 10 |( 9)11 4 | 11 4 3 10 | |20/ \ 8| 4/ \ | / \ |
3807! | 12 | 17 4 12 |(11)13 1 | 13 1 4 12 | | / \ | / 1 \ | / \ |
3808! +----+----------+----------+-------------+ |/ 16 \|/ \|/ 10 \|
3809! | 13 | 14 7 2 |( 8) 1 6 | 1 6 7 2 |13-------1-------2-------8
3810! | 14 | 15 9 3 |(10) 2 8 | 2 8 9 3 | |\ 5 /|
3811! | 15 | 16 11 4 |(12) 3 10 | 3 10 11 4 | | \ / |
3812! | 16 | 17 13 1 |( 6) 4 12 | 4 12 13 1 | |9 \ /13|
3813! +----+----------+----------+-------------+ | 14 |
3814! | 17 | 14 6 7 | - 2 1 | 2 1 6 7 | | / \ |
3815! | 18 | 15 8 9 | - 3 2 | 3 2 8 9 | | / \ |
3816! | 19 | 16 10 11 | - 4 3 | 4 3 10 11 | |/ 17 \|
3817! | 20 | 17 12 13 | - 1 4 | 1 4 12 13 | 6-------7
3818! +----+----------+----------+-------------+
3819 DATA ic /
3820 1 5, 1, 2 , 14, 3, 4 , 3, 4, 1, 2,
3821 2 5, 2, 3 , 15, 4, 1 , 4, 1, 2, 3,
3822 3 5, 3, 4 , 16, 1, 2 , 1, 2, 3, 4,
3823 4 5, 4, 1 , 17, 2, 3 , 2, 3, 4, 1,
3824 5 14, 2, 1 , 5, 6, 7 , 6, 7, 2, 1,
3825 6 15, 3, 2 , 5, 8, 9 , 8, 9, 3, 2,
3826 7 16, 4, 3 , 5,10,11 , 10,11, 4, 3,
3827 8 17, 1, 4 , 5,12,13 , 12,13, 1, 4,
3828 9 14, 1, 6 , 0, 7, 2 , 7, 2, 1, 6,
3829 . 15, 2, 8 , 0, 9, 3 , 9, 3, 2, 8,
3830 1 16, 3,10 , 0,11, 4 , 11, 4, 3,10,
3831 2 17, 4,12 , 0,13, 1 , 13, 1, 4,12,
3832 3 14, 7, 2 , 0, 1, 6 , 1, 6, 7, 2,
3833 4 15, 9, 3 , 0, 2, 8 , 2, 8, 9, 3,
3834 5 16,11, 4 , 0, 3,10 , 3,10,11, 4,
3835 6 17,13, 1 , 0, 4,12 , 4,12,13, 1,
3836 7 14, 6, 7 , 0, 2, 1 , 2, 1, 6, 7,
3837 8 15, 8, 9 , 0, 3, 2 , 3, 2, 8, 9,
3838 9 16,10,11 , 0, 4, 3 , 4, 3,10,11,
3839 . 17,12,13 , 0, 1, 4 , 1, 4,12,13/
3840
3841! +------+-------------------------------------------+
3842! | | sous-triangle voisins |
3843! |sous +----------+----------+---------------------+
3844! |trian | n3/=n4 | n3=n4 |
3845! | +----------+----------+----------+----------+
3846! | | IT0 | 2,4,6,8=0| | 2,4,6,8=0|
3847! +------+----------+----------+----------+----------+
3848! | 1 | 5 2 4 | 5 2 4 | 5 6 8 | 5 6 8 |
3849! | 2 | 6 3 1 | 6 3 1 | - - - | - - - |
3850! | 3 | 7 4 2 | 7 4 2 | - - - | - - - |
3851! | 4 | 8 1 3 | 8 1 3 | - - - | - - - |
3852! +------+----------+----------+----------+----------+
3853! | 5 | 1 9 13 | 1 -1 -1 | 1 9 13 | 1 -1 -1 |
3854! | 6 | 2 10 14 | 2 -1 -1 | 1 10 14 | 1 -1 -1 |
3855! | 7 | 3 11 15 | 3 -1 -1 | - - - | - - - |
3856! | 8 | 4 12 16 | 4 -1 -1 | 1 12 16 | 1 -1 -1 |
3857! +------+----------+----------+----------+----------+
3858! | 9 | -1 17 5 | - - - | -1 17 5 | - - - |
3859! | 10 | -1 18 6 | - - - | -1 18 6 | - - - |
3860! | 11 | -1 19 7 | - - - | - - - | - - - |
3861! | 12 | -1 20 8 | - - - | -1 20 8 | - - - |
3862! +------+----------+----------+----------+----------+
3863! | 13 | -1 5 17 | - - - | -1 5 17 | - - - |
3864! | 14 | -1 6 18 | - - - | -1 6 18 | - - - |
3865! | 15 | -1 7 19 | - - - | - - - | - - - |
3866! | 16 | -1 8 20 | - - - | -1 8 20 | - - - |
3867! +------+----------+----------+----------+----------+
3868! | 17 | -1 9 13 | - - - | -1 9 13 | - - - |
3869! | 18 | -1 10 14 | - - - | -1 10 14 | - - - |
3870! | 19 | -1 11 15 | - - - | - - - | - - - |
3871! | 20 | -1 12 16 | - - - | -1 12 16 | - - - |
3872! +------+----------+----------+----------+----------+
3873 DATA it0 /
3874 1 5, 2, 4,
3875 2 6, 3, 1,
3876 3 7, 4, 2,
3877 4 8, 1, 3,
3878 5 1, 9,13,
3879 6 2,10,14,
3880 7 3,11,15,
3881 8 4,12,16,
3882 9 -1,17, 5,
3883 . -1,18, 6,
3884 1 -1,19, 7,
3885 2 -1,20, 8,
3886 3 -1, 5,17,
3887 4 -1, 6,18,
3888 5 -1, 7,19,
3889 6 -1, 8,20,
3890 7 -1, 9,13,
3891 8 -1,10,14,
3892 9 -1,11,15,
3893 . -1,12,16/
3894
3895 unp = one + em4
3896 zerom = zero - em4
3897c EPS = (TWO+HALF)/HUNDRED
3898
3899 itq = subtria
3900 k = ic(1,itq)
3901 xa = xxl(k)
3902 ya = yyl(k)
3903 za = zzl(k)
3904
3905 k = ic(2,itq)
3906 xb = xxl(k)
3907 yb = yyl(k)
3908 zb = zzl(k)
3909
3910 k = ic(3,itq)
3911 xc = xxl(k)
3912 yc = yyl(k)
3913 zc = zzl(k)
3914
3915 k = ic(4,itq)
3916 xd = xxl(k)
3917 yd = yyl(k)
3918 zd = zzl(k)
3919
3920 IF(itq ==1)THEN
3921 n1 = sx125
3922 n2 = sy125
3923 n3 = sz125
3924 ELSEIF(itq ==2)THEN
3925 n1 = sx235
3926 n2 = sy235
3927 n3 = sz235
3928 ELSEIF(itq ==3)THEN
3929 n1 = sx345
3930 n2 = sy345
3931 n3 = sz345
3932 ELSEIF(itq ==4)THEN
3933 n1 = sx415
3934 n2 = sy415
3935 n3 = sz415
3936 ENDIF
3937
3938 s2 = one/max(em30,sqrt(n1*n1 + n2*n2 + n3*n3))
3939 n1 = n1*s2
3940 n2 = n2*s2
3941 n3 = n3*s2
3942
3943 bbb = (xi-xa)*n1+(yi-ya)*n2+(zi-za)*n3
3944 pene = -min(zero,bbb)
3945
3946 xab = xb-xa
3947 yab = yb-ya
3948 zab = zb-za
3949 xbc = xc-xb
3950 ybc = yc-yb
3951 zbc = zc-zb
3952 xca = xa-xc
3953 yca = ya-yc
3954 zca = za-zc
3955
3956 xbd = xd-xb
3957 ybd = yd-yb
3958 zbd = zd-zb
3959
3960 xia = xa-xi
3961 yia = ya-yi
3962 zia = za-zi
3963 xib = xb-xi
3964 yib = yb-yi
3965 zib = zb-zi
3966 xic = xc-xi
3967 yic = yc-yi
3968 zic = zc-zi
3969 sx = - yab*zca + zab*yca
3970 sy = - zab*xca + xab*zca
3971 sz = - xab*yca + yab*xca
3972 s2 = sx*sx+sy*sy+sz*sz
3973 sax = yib*zic - zib*yic
3974 say = zib*xic - xib*zic
3975 saz = xib*yic - yib*xic
3976 la = (sx*sax+sy*say+sz*saz)/s2
3977 sbx = yic*zia - zic*yia
3978 sby = zic*xia - xic*zia
3979 sbz = xic*yia - yic*xia
3980 lb = (sx*sbx+sy*sby+sz*sbz)/s2
3981 lc = one - la - lb
3982
3983c------------------------------------------------------
3984c quadrangles
3985 aaa = zerom
3986 IF(lb<aaa)THEN
3987c outside side CA
3988 istep=3
3989 subtria_n=it0(2,itq)
3990 subtria =it0(2,itq)
3991 izlim=-1
3992c CALL I24NEXTTRIA2(
3993c 1 IZLIM,ISTEP,SUBTRIA_N,SUBTRIA,
3994c 2 PENE ,XXL ,YYL ,ZZL ,
3995c 3 SX125,SY125,SZ125,SX235,SY235,
3996c 4 SZ235,SX345,SY345,SZ345,SX415,
3997c 5 SY415,SZ415,XI ,YI ,ZI ,
3998c 6 N1 ,N2 ,N3 ,LA ,LB ,
3999c 7 LC )
4000 ELSEIF(lc<aaa)THEN
4001c outside side AB
4002 istep=3
4003 subtria_n=it0(3,itq)
4004 subtria =it0(3,itq)
4005 izlim=-1
4006c CALL I24NEXTTRIA2(
4007c 1 IZLIM,ISTEP,SUBTRIA_N,SUBTRIA,
4008c 2 PENE ,XXL ,YYL ,ZZL ,
4009c 3 SX125,SY125,SZ125,SX235,SY235,
4010c 4 SZ235,SX345,SY345,SZ345,SX415,
4011c 5 SY415,SZ415,XI ,YI ,ZI ,
4012c 6 N1 ,N2 ,N3 ,LA ,LB ,
4013c 7 LC )
4014c check if outside side BC
4015 ELSEIF(la<-eps)THEN
4016c hors zone d'extension de surface
4017 istep=3
4018 subtria_n=it0(1,itq)
4019 izlim=-1
4020 ELSEIF(la<eps)THEN
4021c zone limite interpolation des normales si convex
4022 vol = n1*xbd + n2*ybd + n3*zbd
4023 IF (gap2<eps.AND.(lb<eps.OR.lc<eps)) vol=-abs(vol)
4024 IF(largep == 1)THEN
4025c large penetration inside +- extension zone
4026 istep=3
4027 subtria_n=it0(1,itq)
4028 izlim = -1
4029 ELSEIF(vol < zero)THEN
4030c convex
4031 izlim=1
4032 ELSEIF(la<zerom)THEN
4033 istep=3
4034 subtria_n=it0(1,itq)
4035 izlim=-1
4036 ENDIF
4037 ENDIF
4038
4039 RETURN
4040 END
4041!||====================================================================
4042!|| s_in_m4 ../engine/source/interfaces/int24/i24dst3.F
4043!||--- called by ------------------------------------------------------
4044!|| i24dst3 ../engine/source/interfaces/int24/i24dst3.F
4045!||====================================================================
4046 SUBROUTINE s_in_m4(X1,Y1,Z1,X2,Y2,Z2,X3,Y3,Z3,X4,Y4,Z4,
4047 1 XI,YI,ZI,IER)
4048C-----------------------------------------------
4049C I m p l i c i t T y p e s
4050C-----------------------------------------------
4051#include "implicit_f.inc"
4052C-----------------------------------------------
4053C D u m m y A r g u m e n t s
4054C-----------------------------------------------
4055 INTEGER IER
4056C REAL
4057 my_real
4058 . x1,x2,x3,x4,
4059 . y1,y2,y3,y4,z1,z2,z3,z4,xi,yi,zi
4060C-----------------------------------------------
4061C L o c a l V a r i a b l e s
4062C-----------------------------------------------
4063C REAL
4064 my_real
4065 . x0, y0, z0, xl1, xl2, xl3, xl4, yy1, yy2, yy3, yy4,
4066 . zz1, zz2, zz3, zz4, xi1, xi2, xi3, xi4, yi1, yi2, yi3, yi4,
4067 . zi1, zi2, zi3, zi4, xn1, yn1, zn1, xn2, yn2, zn2, xn3, yn3,
4068 . zn3, xn4, yn4, zn4, an, area, a12, a23, a34, a41, b12, b23,
4069 . b34, b41, ab1, ab2, tp, tm, sp, sm, n1,n2,n3,alp,
4070 . nn(3), ssc, ttc
4071C------------------------------------------------
4072 ier = 0
4073 alp=em8
4074 x0 = fourth*(x1+x2+x3+x4)
4075 y0 = fourth*(y1+y2+y3+y4)
4076 z0 = fourth*(z1+z2+z3+z4)
4077C
4078 xl1 = x1-x0
4079 xl2 = x2-x0
4080 xl3 = x3-x0
4081 xl4 = x4-x0
4082 yy1 = y1-y0
4083 yy2 = y2-y0
4084 yy3 = y3-y0
4085 yy4 = y4-y0
4086 zz1 = z1-z0
4087 zz2 = z2-z0
4088 zz3 = z3-z0
4089 zz4 = z4-z0
4090C
4091 xi1 = x1-xi
4092 xi2 = x2-xi
4093 xi3 = x3-xi
4094 xi4 = x4-xi
4095 yi1 = y1-yi
4096 yi2 = y2-yi
4097 yi3 = y3-yi
4098 yi4 = y4-yi
4099 zi1 = z1-zi
4100 zi2 = z2-zi
4101 zi3 = z3-zi
4102 zi4 = z4-zi
4103C
4104 xn1 = yy1*zz2 - yy2*zz1
4105 yn1 = zz1*xl2 - zz2*xl1
4106 zn1 = xl1*yy2 - xl2*yy1
4107 n1=xn1
4108 n2=yn1
4109 n3=zn1
4110C
4111 xn2 = yy2*zz3 - yy3*zz2
4112 yn2 = zz2*xl3 - zz3*xl2
4113 zn2 = xl2*yy3 - xl3*yy2
4114 n1=n1+xn2
4115 n2=n2+yn2
4116 n3=n3+zn2
4117C
4118 xn3 = yy3*zz4 - yy4*zz3
4119 yn3 = zz3*xl4 - zz4*xl3
4120 zn3 = xl3*yy4 - xl4*yy3
4121 n1=n1+xn3
4122 n2=n2+yn3
4123 n3=n3+zn3
4124C
4125 xn4 = yy4*zz1 - yy1*zz4
4126 yn4 = zz4*xl1 - zz1*xl4
4127 zn4 = xl4*yy1 - xl1*yy4
4128 n1=n1+xn4
4129 n2=n2+yn4
4130 n3=n3+zn4
4131C
4132 an= max(em20,sqrt(n1*n1+n2*n2+n3*n3))
4133 n1=n1/an
4134 n2=n2/an
4135 n3=n3/an
4136 nn(1)=n1
4137 nn(2)=n2
4138 nn(3)=n3
4139 IF(an<=em19) THEN
4140 ier=-1
4141 RETURN
4142 ENDIF
4143 area=half*an
4144C
4145 a12=(n1*xn1+n2*yn1+n3*zn1)
4146 a23=(n1*xn2+n2*yn2+n3*zn2)
4147 a34=(n1*xn3+n2*yn3+n3*zn3)
4148 a41=(n1*xn4+n2*yn4+n3*zn4)
4149C
4150 xn1 = yi1*zi2 - yi2*zi1
4151 yn1 = zi1*xi2 - zi2*xi1
4152 zn1 = xi1*yi2 - xi2*yi1
4153 b12=(n1*xn1+n2*yn1+n3*zn1)
4154C
4155 xn2 = yi2*zi3 - yi3*zi2
4156 yn2 = zi2*xi3 - zi3*xi2
4157 zn2 = xi2*yi3 - xi3*yi2
4158 b23=(n1*xn2+n2*yn2+n3*zn2)
4159C
4160 xn3 = yi3*zi4 - yi4*zi3
4161 yn3 = zi3*xi4 - zi4*xi3
4162 zn3 = xi3*yi4 - xi4*yi3
4163 b34=(n1*xn3+n2*yn3+n3*zn3)
4164C
4165 xn4 = yi4*zi1 - yi1*zi4
4166 yn4 = zi4*xi1 - zi1*xi4
4167 zn4 = xi4*yi1 - xi1*yi4
4168 b41=(n1*xn4+n2*yn4+n3*zn4)
4169C
4170 ab1=a23*b41
4171 ab2=b23*a41
4172C
4173 IF(abs(ab1+ab2)/area>em10)THEN
4174 ssc=(ab1-ab2)/(ab1+ab2)
4175 ELSE
4176 ssc=zero
4177 ENDIF
4178 IF(abs(a34/area)>em10)THEN
4179 ab1=b12*a34
4180 ab2=b34*a12
4181 ttc=(ab1-ab2)/(ab1+ab2)
4182 ELSE
4183 ttc=b12/a12-one
4184 IF(b23<=zero.AND.b41<=zero)THEN
4185 IF(-b23/a12<=alp.AND.-b41/a12<=alp)ssc=zero
4186 ELSEIF(b23<=zero)THEN
4187 IF(-b23/a12<=alp) THEN
4188 ssc=one
4189 ELSE
4190 ssc=two
4191 END IF
4192 ELSEIF(b41<=zero)THEN
4193 IF(-b41/a12<=alp) THEN
4194 ssc=-one
4195 ELSE
4196 ssc=-two
4197 END IF
4198 ENDIF
4199 ENDIF
4200C
4201 IF(abs(ssc)>one+alp.OR.abs(ttc)>one+alp) THEN
4202 ier=1
4203 ENDIF
4204C
4205 RETURN
4206 END
4207!||====================================================================
4208!|| intersecv ../engine/source/interfaces/int24/i24dst3.F
4209!||--- called by ------------------------------------------------------
4210!|| i24dst3 ../engine/source/interfaces/int24/i24dst3.F
4211!||====================================================================
4212 SUBROUTINE intersecv(
4213 1 ISTEP ,SUBTRIA,IXX3 ,IXX4 ,
4214 1 INTERSECT,V12 ,V23 ,V34 ,V41 ,
4215 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
4216 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
4217 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
4218 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,ADT0 ,
4219 6 VXI ,VYI ,VZI ,XO1 ,XO2 ,
4220 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
4221 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
4222 9 ZO3 ,ZO4 ,ZO5 ,NX ,NY ,
4223 A NZ ,SX125 ,SX235 ,SX345 ,SX415 ,
4224 B SY125 ,SY235 ,SY345 ,SY415 ,SZ125 ,
4225 C SZ235 ,SZ345 ,SZ415 ,XI ,YI ,
4226 D ZI ,ZEROM ,UNP ,ZEROMT,UNPT )
4227C-----------------------------------------------
4228C I m p l i c i t T y p e s
4229C-----------------------------------------------
4230#include "implicit_f.inc"
4231C-----------------------------------------------
4232C D u m m y A r g u m e n t s
4233C-----------------------------------------------
4234 INTEGER ISTEP ,SUBTRIA,IXX3 ,IXX4 ,INTERSECT
4235C REAL
4236 my_real
4237 1 v12 ,v23 ,v34 ,v41 ,
4238 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
4239 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
4240 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
4241 5 zz4 ,xx5 ,yy5 ,zz5 ,adt0 ,
4242 6 vxi ,vyi ,vzi ,xo1 ,xo2 ,
4243 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
4244 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
4245 9 zo3 ,zo4 ,zo5 ,nx ,ny ,
4246 a nz ,sx125 ,sx235 ,sx345 ,sx415 ,
4247 b sy125 ,sy235 ,sy345 ,sy415 ,sz125 ,
4248 c sz235 ,sz345 ,sz415 ,zerom ,unp ,
4249 d zeromt,unpt ,xi,yi,zi
4250C-----------------------------------------------
4251C C o m m o n B l o c k s
4252C-----------------------------------------------
4253#include "com08_c.inc"
4254C-----------------------------------------------
4255C L o c a l V a r i a b l e s
4256C-----------------------------------------------
4257 my_real
4258 . x51, x52, x53, x54,
4259 . y51, y52, y53, y54,
4260 . z51, z52, z53, z54,
4261 . xi1, xi2, xi3, xi4, xi5,
4262 . yi1, yi2, yi3, yi4, yi5,
4263 . zi1, zi2, zi3, zi4, zi5,
4264 . xia, xib, xic,
4265 . yia, yib, yic,
4266 . zia, zib, zic,
4267 . xoi,yoi,zoi,xs,ys,zs,
4268 . xm1, xm2, xm3, xm4, xm5,
4269 . ym1, ym2, ym3, ym4, ym5,
4270 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
4271 my_real
4272 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
4273 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
4274 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,adt,s2
4275C------------------------------------------------
4276 IF(vo12 <= zero .and. v12 >= zero )THEN
4277 IF ((abs(sx125)+abs(sy125)+abs(sz125))>em20) THEN
4278 IF(abs(vo12) < em20)THEN
4279 aaa = one
4280 ELSEIF(abs(v12) < em20)THEN
4281 aaa = zero
4282 ELSE
4283 aaa = v12 / (v12-vo12)
4284 ENDIF
4285
4286 adt = aaa*dt1
4287
4288 xm1 = xx1 + aaa*(xo1-xx1)
4289 ym1 = yy1 + aaa*(yo1-yy1)
4290 zm1 = zz1 + aaa*(zo1-zz1)
4291
4292 xm2 = xx2 + aaa*(xo2-xx2)
4293 ym2 = yy2 + aaa*(yo2-yy2)
4294 zm2 = zz2 + aaa*(zo2-zz2)
4295
4296 xoi = xi - adt*vxi
4297 yoi = yi - adt*vyi
4298 zoi = zi - adt*vzi
4299
4300 xm5 = xx5 + aaa*(xo5-xx5)
4301 ym5 = yy5 + aaa*(yo5-yy5)
4302 zm5 = zz5 + aaa*(zo5-zz5)
4303
4304 x51 = xm1 - xm5
4305 y51 = ym1 - ym5
4306 z51 = zm1 - zm5
4307
4308 x52 = xm2 - xm5
4309 y52 = ym2 - ym5
4310 z52 = zm2 - zm5
4311
4312 xi1 = xm1 - xoi
4313 yi1 = ym1 - yoi
4314 zi1 = zm1 - zoi
4315
4316 xi2 = xm2 - xoi
4317 yi2 = ym2 - yoi
4318 zi2 = zm2 - zoi
4319
4320 xi5 = xm5 - xoi
4321 yi5 = ym5 - yoi
4322 zi5 = zm5 - zoi
4323
4324 sx0 = y51*z52 - z51*y52
4325 sy0 = z51*x52 - x51*z52
4326 sz0 = x51*y52 - y51*x52
4327
4328 sx1 = yi5*zi1 - zi5*yi1
4329 sy1 = zi5*xi1 - xi5*zi1
4330 sz1 = xi5*yi1 - yi5*xi1
4331
4332 sx5 = yi1*zi2 - zi1*yi2
4333 sy5 = zi1*xi2 - xi1*zi2
4334 sz5 = xi1*yi2 - yi1*xi2
4335
4336 sx2 = yi2*zi5 - zi2*yi5
4337 sy2 = zi2*xi5 - xi2*zi5
4338 sz2 = xi2*yi5 - yi2*xi5
4339
4340 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4341 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
4342 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
4343 la0 = one-lb0-lc0
4344
4345 IF(ixx3 /= ixx4)THEN
4346 IF(lb0 >= zerom .and. lb0 <= unp .and.
4347 . lc0 >= zerom .and. lc0 <= unp .and.
4348 . la0 >= zerom .and. la0 <= unp)THEN
4349 intersect = 1
4350 subtria= 1 ! subtriangle
4351 istep = 5
4352 adt0= adt
4353 nx = sx125
4354 ny = sy125
4355 nz = sz125
4356 ENDIF
4357C----------loose tolerance for tria
4358 ELSE
4359 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
4360 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
4361 . la0 >= zeromt .AND. la0 <= unpt)THEN
4362 intersect = 1
4363 subtria= 1 ! subtriangle
4364 istep = 5
4365 adt0= adt
4366 nx = sx125
4367 ny = sy125
4368 nz = sz125
4369 ENDIF
4370 END if! (IXX(I,3) /= IXX(I,4))THEN
4371 END if! ((ABS(SX125)+ABS(SY125)+ABS(SZ125))>EM20) THEN
4372 ENDIF
4373 IF(ixx3 /= ixx4)THEN
4374 IF(vo23 <= zero .and. v23 >= zero )THEN
4375 IF ((abs(sx235)+abs(sy235)+abs(sz235))>em20) THEN
4376 IF(abs(vo23) < em20)THEN
4377 aaa = one
4378 ELSEIF(abs(v23) < em20)THEN
4379 aaa = zero
4380 ELSE
4381 aaa = v23 / (v23-vo23)
4382 ENDIF
4383
4384 adt = aaa*dt1
4385
4386 xm2 = xx2 + aaa*(xo2-xx2)
4387 ym2 = yy2 + aaa*(yo2-yy2)
4388 zm2 = zz2 + aaa*(zo2-zz2)
4389
4390 xm3 = xx3 + aaa*(xo3-xx3)
4391 ym3 = yy3 + aaa*(yo3-yy3)
4392 zm3 = zz3 + aaa*(zo3-zz3)
4393
4394 xoi = xi - adt*vxi
4395 yoi = yi - adt*vyi
4396 zoi = zi - adt*vzi
4397
4398 xm5 = xx5 + aaa*(xo5-xx5)
4399 ym5 = yy5 + aaa*(yo5-yy5)
4400 zm5 = zz5 + aaa*(zo5-zz5)
4401
4402 x52 = xm2 - xm5
4403 y52 = ym2 - ym5
4404 z52 = zm2 - zm5
4405
4406 x53 = xm3 - xm5
4407 y53 = ym3 - ym5
4408 z53 = zm3 - zm5
4409
4410 xi2 = xm2 - xoi
4411 yi2 = ym2 - yoi
4412 zi2 = zm2 - zoi
4413
4414 xi3 = xm3 - xoi
4415 yi3 = ym3 - yoi
4416 zi3 = zm3 - zoi
4417
4418 xi5 = xm5 - xoi
4419 yi5 = ym5 - yoi
4420 zi5 = zm5 - zoi
4421
4422 sx0 = y52*z53 - z52*y53
4423 sy0 = z52*x53 - x52*z53
4424 sz0 = x52*y53 - y52*x53
4425
4426 sx2 = yi5*zi2 - zi5*yi2
4427 sy2 = zi5*xi2 - xi5*zi2
4428 sz2 = xi5*yi2 - yi5*xi2
4429
4430 sx5 = yi2*zi3 - zi2*yi3
4431 sy5 = zi2*xi3 - xi2*zi3
4432 sz5 = xi2*yi3 - yi2*xi3
4433
4434 sx3 = yi3*zi5 - zi3*yi5
4435 sy3 = zi3*xi5 - xi3*zi5
4436 sz3 = xi3*yi5 - yi3*xi5
4437
4438 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4439 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
4440 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
4441 la0 = one-lb0-lc0
4442
4443 IF(lb0 >= zerom .and. lb0 <= unp .and.
4444 . lc0 >= zerom .and. lc0 <= unp .and.
4445 . la0 >= zerom .and. la0 <= unp)THEN
4446 intersect = 1
4447
4448 IF(adt > adt0)THEN
4449 subtria= 2 ! subtriangle
4450 istep = 5
4451 adt0= adt
4452 nx = sx235
4453 ny = sy235
4454 nz = sz235
4455 ENDIF
4456 ENDIF
4457 END IF !((ABS(SX235)+ABS(SY235)+ABS(SZ235))>EM20) THEN
4458 ENDIF
4459 IF(vo34 <= zero .and. v34 >= zero )THEN
4460 IF ((abs(sx345)+abs(sy345)+abs(sz345))>em20) THEN
4461 IF(abs(vo34) < em20)THEN
4462 aaa = one
4463 ELSEIF(abs(v34) < em20)THEN
4464 aaa = zero
4465 ELSE
4466 aaa = v34 / (v34-vo34)
4467 ENDIF
4468
4469 adt = aaa*dt1
4470
4471 xm3 = xx3 + aaa*(xo3-xx3)
4472 ym3 = yy3 + aaa*(yo3-yy3)
4473 zm3 = zz3 + aaa*(zo3-zz3)
4474
4475 xm4 = xx4 + aaa*(xo4-xx4)
4476 ym4 = yy4 + aaa*(yo4-yy4)
4477 zm4 = zz4 + aaa*(zo4-zz4)
4478
4479 xoi = xi - adt*vxi
4480 yoi = yi - adt*vyi
4481 zoi = zi - adt*vzi
4482
4483 xm5 = xx5 + aaa*(xo5-xx5)
4484 ym5 = yy5 + aaa*(yo5-yy5)
4485 zm5 = zz5 + aaa*(zo5-zz5)
4486
4487 x53 = xm3 - xm5
4488 y53 = ym3 - ym5
4489 z53 = zm3 - zm5
4490
4491 x54 = xm4 - xm5
4492 y54 = ym4 - ym5
4493 z54 = zm4 - zm5
4494
4495 xi3 = xm3 - xoi
4496 yi3 = ym3 - yoi
4497 zi3 = zm3 - zoi
4498
4499 xi4 = xm4 - xoi
4500 yi4 = ym4 - yoi
4501 zi4 = zm4 - zoi
4502
4503 xi5 = xm5 - xoi
4504 yi5 = ym5 - yoi
4505 zi5 = zm5 - zoi
4506
4507 sx0 = y53*z54 - z53*y54
4508 sy0 = z53*x54 - x53*z54
4509 sz0 = x53*y54 - y53*x54
4510
4511 sx3 = yi5*zi3 - zi5*yi3
4512 sy3 = zi5*xi3 - xi5*zi3
4513 sz3 = xi5*yi3 - yi5*xi3
4514
4515 sx5 = yi3*zi4 - zi3*yi4
4516 sy5 = zi3*xi4 - xi3*zi4
4517 sz5 = xi3*yi4 - yi3*xi4
4518
4519 sx4 = yi4*zi5 - zi4*yi5
4520 sy4 = zi4*xi5 - xi4*zi5
4521 sz4 = xi4*yi5 - yi4*xi5
4522
4523 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4524 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
4525 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
4526 la0 = one-lb0-lc0
4527
4528 IF(lb0 >= zerom .and. lb0 <= unp .and.
4529 . lc0 >= zerom .and. lc0 <= unp .and.
4530 . la0 >= zerom .and. la0 <= unp)THEN
4531 intersect = 1
4532 IF(adt > adt0)THEN
4533 subtria= 3 ! subtriangle
4534 istep = 5
4535 adt0= adt
4536 nx = sx345
4537 ny = sy345
4538 nz = sz345
4539 ENDIF
4540 ENDIF
4541 END IF !((ABS(SX345)+ABS(SY345)+ABS(SZ345))>EM20) THEN
4542 ENDIF
4543
4544 IF(vo41 <= zero .and. v41 >= zero )THEN
4545 IF ((abs(sx415)+abs(sy415)+abs(sz415))>em20) THEN
4546 IF(abs(vo41) < em20)THEN
4547 aaa = one
4548 ELSEIF(abs(v41) < em20)THEN
4549 aaa = zero
4550 ELSE
4551 aaa = v41 / (v41-vo41)
4552 ENDIF
4553
4554 adt = aaa*dt1
4555
4556 xm4 = xx4 + aaa*(xo4-xx4)
4557 ym4 = yy4 + aaa*(yo4-yy4)
4558 zm4 = zz4 + aaa*(zo4-zz4)
4559
4560 xm1 = xx1 + aaa*(xo1-xx1)
4561 ym1 = yy1 + aaa*(yo1-yy1)
4562 zm1 = zz1 + aaa*(zo1-zz1)
4563
4564 xoi = xi - adt*vxi
4565 yoi = yi - adt*vyi
4566 zoi = zi - adt*vzi
4567
4568 xm5 = xx5 + aaa*(xo5-xx5)
4569 ym5 = yy5 + aaa*(yo5-yy5)
4570 zm5 = zz5 + aaa*(zo5-zz5)
4571
4572 x54 = xm4 - xm5
4573 y54 = ym4 - ym5
4574 z54 = zm4 - zm5
4575
4576 x51 = xm1 - xm5
4577 y51 = ym1 - ym5
4578 z51 = zm1 - zm5
4579
4580 xi4 = xm4 - xoi
4581 yi4 = ym4 - yoi
4582 zi4 = zm4 - zoi
4583
4584 xi1 = xm1 - xoi
4585 yi1 = ym1 - yoi
4586 zi1 = zm1 - zoi
4587
4588 xi5 = xm5 - xoi
4589 yi5 = ym5 - yoi
4590 zi5 = zm5 - zoi
4591
4592 sx0 = y54*z51 - z54*y51
4593 sy0 = z54*x51 - x54*z51
4594 sz0 = x54*y51 - y54*x51
4595
4596 sx4 = yi5*zi4 - zi5*yi4
4597 sy4 = zi5*xi4 - xi5*zi4
4598 sz4 = xi5*yi4 - yi5*xi4
4599
4600 sx5 = yi4*zi1 - zi4*yi1
4601 sy5 = zi4*xi1 - xi4*zi1
4602 sz5 = xi4*yi1 - yi4*xi1
4603
4604 sx1 = yi1*zi5 - zi1*yi5
4605 sy1 = zi1*xi5 - xi1*zi5
4606 sz1 = xi1*yi5 - yi1*xi5
4607
4608 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4609 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
4610 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
4611 la0 = one-lb0-lc0
4612
4613 IF(lb0 >= zerom .and. lb0 <= unp .and.
4614 . lc0 >= zerom .and. lc0 <= unp .and.
4615 . la0 >= zerom .and. la0 <= unp)THEN
4616 intersect = 1
4617 IF(adt > adt0)THEN
4618 subtria= 4 ! subtriangle
4619 istep = 5
4620 adt0= adt
4621 nx = sx415
4622 ny = sy415
4623 nz = sz415
4624 ENDIF
4625 ENDIF
4626 END IF !((ABS(SX415)+ABS(SY415)+ABS(SZ415))>EM20) THEN
4627 ENDIF
4628 ENDIF !IF(IXX3 /= IXX4)THEN
4629C
4630 RETURN
4631 END
4632!||====================================================================
4633!|| intersecv0 ../engine/source/interfaces/int24/i24dst3.F
4634!||--- called by ------------------------------------------------------
4635!|| i24dst3 ../engine/source/interfaces/int24/i24dst3.F
4636!||====================================================================
4637 SUBROUTINE intersecv0(
4638 1 ISTEP ,SUBTRIA,IXX3 ,IXX4 ,
4639 1 INTERSECT,V12 ,V23 ,V34 ,V41 ,
4640 2 VO12 ,VO23 ,VO34 ,VO41 ,XX1 ,
4641 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
4642 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
4643 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,ADT0 ,
4644 6 XOI ,YOI ,ZOI ,XO1 ,XO2 ,
4645 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
4646 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
4647 9 ZO3 ,ZO4 ,ZO5 ,NX ,NY ,
4648 A NZ ,SX125 ,SX235 ,SX345 ,SX415 ,
4649 B SY125 ,SY235 ,SY345 ,SY415 ,SZ125 ,
4650 C SZ235 ,SZ345 ,SZ415 ,XI ,YI ,
4651 D ZI ,ZEROMT ,UNPT )
4652C-----------------------------------------------
4653C I m p l i c i t T y p e s
4654C-----------------------------------------------
4655#include "implicit_f.inc"
4656C-----------------------------------------------
4657C D u m m y A r g u m e n t s
4658C-----------------------------------------------
4659 INTEGER ISTEP ,SUBTRIA,IXX3 ,IXX4 ,INTERSECT
4660C REAL
4661 my_real
4662 1 v12 ,v23 ,v34 ,v41 ,
4663 2 vo12 ,vo23 ,vo34 ,vo41 ,xx1 ,
4664 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
4665 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
4666 5 zz4 ,xx5 ,yy5 ,zz5 ,adt0 ,
4667 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
4668 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
4669 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
4670 9 zo3 ,zo4 ,zo5 ,nx ,ny ,
4671 a nz ,sx125 ,sx235 ,sx345 ,sx415 ,
4672 b sy125 ,sy235 ,sy345 ,sy415 ,sz125 ,
4673 c sz235 ,sz345 ,sz415 ,
4674 d zeromt,unpt ,xi,yi,zi
4675C-----------------------------------------------
4676C C o m m o n B l o c k s
4677C-----------------------------------------------
4678#include "com08_c.inc"
4679C-----------------------------------------------
4680C L o c a l V a r i a b l e s
4681C-----------------------------------------------
4682 my_real
4683 . x51, x52, x53, x54,
4684 . y51, y52, y53, y54,
4685 . z51, z52, z53, z54,
4686 . xi1, xi2, xi3, xi4, xi5,
4687 . yi1, yi2, yi3, yi4, yi5,
4688 . zi1, zi2, zi3, zi4, zi5,
4689 . xia, xib, xic,
4690 . yia, yib, yic,
4691 . zia, zib, zic,
4692 . xs,ys,zs,
4693 . xm1, xm2, xm3, xm4, xm5,
4694 . ym1, ym2, ym3, ym4, ym5,
4695 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
4696 my_real
4697 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
4698 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
4699 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,adt,s2
4700C------------------------------------------------
4701C---------------- cas V>0 V_old>0
4702 IF(vo12 > zero .and. v12 >= zero )THEN
4703 IF ((abs(sx125)+abs(sy125)+abs(sz125))>em20) THEN
4704 aaa = one
4705 adt = dt1
4706
4707 xm1 = xo1
4708 ym1 = yo1
4709 zm1 = zo1
4710
4711 xm2 = xo2
4712 ym2 = yo2
4713 zm2 = zo2
4714
4715c XOI = XI - ADT*VXI
4716c YOI = YI - ADT*VYI
4717c ZOI = ZI - ADT*VZI
4718
4719 xm5 = xo5
4720 ym5 = yo5
4721 zm5 = zo5
4722
4723 x51 = xm1 - xm5
4724 y51 = ym1 - ym5
4725 z51 = zm1 - zm5
4726
4727 x52 = xm2 - xm5
4728 y52 = ym2 - ym5
4729 z52 = zm2 - zm5
4730
4731 xi1 = xm1 - xoi
4732 yi1 = ym1 - yoi
4733 zi1 = zm1 - zoi
4734
4735 xi2 = xm2 - xoi
4736 yi2 = ym2 - yoi
4737 zi2 = zm2 - zoi
4738
4739 xi5 = xm5 - xoi
4740 yi5 = ym5 - yoi
4741 zi5 = zm5 - zoi
4742
4743 sx0 = y51*z52 - z51*y52
4744 sy0 = z51*x52 - x51*z52
4745 sz0 = x51*y52 - y51*x52
4746
4747 sx1 = yi5*zi1 - zi5*yi1
4748 sy1 = zi5*xi1 - xi5*zi1
4749 sz1 = xi5*yi1 - yi5*xi1
4750
4751 sx5 = yi1*zi2 - zi1*yi2
4752 sy5 = zi1*xi2 - xi1*zi2
4753 sz5 = xi1*yi2 - yi1*xi2
4754
4755 sx2 = yi2*zi5 - zi2*yi5
4756 sy2 = zi2*xi5 - xi2*zi5
4757 sz2 = xi2*yi5 - yi2*xi5
4758
4759 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4760 lb0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
4761 lc0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
4762 la0 = one-lb0-lc0
4763
4764 IF(lb0 >= zeromt .AND. lb0 <= unpt .AND.
4765 . lc0 >= zeromt .AND. lc0 <= unpt .AND.
4766 . la0 >= zeromt .AND. la0 <= unpt)THEN
4767 intersect = 1
4768 subtria= 1 ! subtriangle
4769 istep = 5
4770 adt0= adt
4771 nx = sx125
4772 ny = sy125
4773 nz = sz125
4774 ENDIF
4775 END IF !((ABS(SX125)+ABS(SY125)+ABS(SZ125))>EM20) THEN
4776 ENDIF
4777 IF(ixx3 /= ixx4)THEN
4778 IF(vo23 > zero .and. v23 >= zero )THEN
4779 IF ((abs(sx235)+abs(sy235)+abs(sz235))>em20) THEN
4780 aaa = one
4781 adt = dt1
4782
4783 xm2 = xo2
4784 ym2 = yo2
4785 zm2 = zo2
4786
4787 xm3 = xo3
4788 ym3 = yo3
4789 zm3 = zo3
4790
4791c XOI = XI - ADT*VXI
4792c YOI = YI - ADT*VYI
4793c ZOI = ZI - ADT*VZI
4794
4795 xm5 = xo5
4796 ym5 = yo5
4797 zm5 = zo5
4798
4799 x52 = xm2 - xm5
4800 y52 = ym2 - ym5
4801 z52 = zm2 - zm5
4802
4803 x53 = xm3 - xm5
4804 y53 = ym3 - ym5
4805 z53 = zm3 - zm5
4806
4807 xi2 = xm2 - xoi
4808 yi2 = ym2 - yoi
4809 zi2 = zm2 - zoi
4810
4811 xi3 = xm3 - xoi
4812 yi3 = ym3 - yoi
4813 zi3 = zm3 - zoi
4814
4815 xi5 = xm5 - xoi
4816 yi5 = ym5 - yoi
4817 zi5 = zm5 - zoi
4818
4819 sx0 = y52*z53 - z52*y53
4820 sy0 = z52*x53 - x52*z53
4821 sz0 = x52*y53 - y52*x53
4822
4823 sx2 = yi5*zi2 - zi5*yi2
4824 sy2 = zi5*xi2 - xi5*zi2
4825 sz2 = xi5*yi2 - yi5*xi2
4826
4827 sx5 = yi2*zi3 - zi2*yi3
4828 sy5 = zi2*xi3 - xi2*zi3
4829 sz5 = xi2*yi3 - yi2*xi3
4830
4831 sx3 = yi3*zi5 - zi3*yi5
4832 sy3 = zi3*xi5 - xi3*zi5
4833 sz3 = xi3*yi5 - yi3*xi5
4834
4835 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4836 lb0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
4837 lc0 = (sx0*sx2 + sy0*sy2 + sz0*sz2) * s2
4838 la0 = one-lb0-lc0
4839
4840 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
4841 . lc0 >= zeromt.and. lc0 <= unpt .and.
4842 . la0 >= zeromt.and. la0 <= unpt)THEN
4843 intersect = 1
4844
4845 IF(adt > adt0)THEN
4846 subtria= 2 ! subtriangle
4847 istep = 5
4848 adt0= adt
4849 nx = sx235
4850 ny = sy235
4851 nz = sz235
4852 ENDIF
4853 ENDIF
4854 END IF !((ABS(SX235)+ABS(SY235)+ABS(SZ235))>EM20) THEN
4855 END if!(VO23 > ZERO .and. V23 >= ZERO )THEN
4856 IF(vo34 > zero .and. v34 >= zero )THEN
4857 IF ((abs(sx345)+abs(sy345)+abs(sz345))>em20) THEN
4858 aaa = one
4859 adt = dt1
4860
4861 xm3 = xo3
4862 ym3 = yo3
4863 zm3 = zo3
4864
4865 xm4 = xo4
4866 ym4 = yo4
4867 zm4 = zo4
4868
4869c XOI = XI - ADT*VXI
4870c YOI = YI - ADT*VYI
4871c ZOI = ZI - ADT*VZI
4872
4873 xm5 = xo5
4874 ym5 = yo5
4875 zm5 = zo5
4876
4877 x53 = xm3 - xm5
4878 y53 = ym3 - ym5
4879 z53 = zm3 - zm5
4880
4881 x54 = xm4 - xm5
4882 y54 = ym4 - ym5
4883 z54 = zm4 - zm5
4884
4885 xi3 = xm3 - xoi
4886 yi3 = ym3 - yoi
4887 zi3 = zm3 - zoi
4888
4889 xi4 = xm4 - xoi
4890 yi4 = ym4 - yoi
4891 zi4 = zm4 - zoi
4892
4893 xi5 = xm5 - xoi
4894 yi5 = ym5 - yoi
4895 zi5 = zm5 - zoi
4896
4897 sx0 = y53*z54 - z53*y54
4898 sy0 = z53*x54 - x53*z54
4899 sz0 = x53*y54 - y53*x54
4900
4901 sx3 = yi5*zi3 - zi5*yi3
4902 sy3 = zi5*xi3 - xi5*zi3
4903 sz3 = xi5*yi3 - yi5*xi3
4904
4905 sx5 = yi3*zi4 - zi3*yi4
4906 sy5 = zi3*xi4 - xi3*zi4
4907 sz5 = xi3*yi4 - yi3*xi4
4908
4909 sx4 = yi4*zi5 - zi4*yi5
4910 sy4 = zi4*xi5 - xi4*zi5
4911 sz4 = xi4*yi5 - yi4*xi5
4912
4913 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4914 lb0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
4915 lc0 = (sx0*sx3 + sy0*sy3 + sz0*sz3) * s2
4916 la0 = one-lb0-lc0
4917
4918 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
4919 . lc0 >= zeromt.and. lc0 <= unpt .and.
4920 . la0 >= zeromt.and. la0 <= unpt)THEN
4921 intersect = 1
4922 IF(adt > adt0)THEN
4923 subtria= 3 ! subtriangle
4924 istep = 5
4925 adt0= adt
4926 nx = sx345
4927 ny = sy345
4928 nz = sz345
4929 ENDIF
4930 ENDIF
4931 END IF !((ABS(SX345)+ABS(SY345)+ABS(SZ345))>EM20) THEN
4932 ENDIF
4933
4934 IF(vo41 > zero .and. v41 >= zero )THEN
4935 IF ((abs(sx415)+abs(sy415)+abs(sz415))>em20) THEN
4936 aaa = one
4937 adt = dt1
4938
4939 xm4 = xo4
4940 ym4 = yo4
4941 zm4 = zo4
4942
4943 xm1 = xo1
4944 ym1 = yo1
4945 zm1 = zo1
4946
4947c XOI = XI - ADT*VXI
4948c YOI = YI - ADT*VYI
4949c ZOI = ZI - ADT*VZI
4950
4951 xm5 = xo5
4952 ym5 = yo5
4953 zm5 = zo5
4954
4955 x54 = xm4 - xm5
4956 y54 = ym4 - ym5
4957 z54 = zm4 - zm5
4958
4959 x51 = xm1 - xm5
4960 y51 = ym1 - ym5
4961 z51 = zm1 - zm5
4962
4963 xi4 = xm4 - xoi
4964 yi4 = ym4 - yoi
4965 zi4 = zm4 - zoi
4966
4967 xi1 = xm1 - xoi
4968 yi1 = ym1 - yoi
4969 zi1 = zm1 - zoi
4970
4971 xi5 = xm5 - xoi
4972 yi5 = ym5 - yoi
4973 zi5 = zm5 - zoi
4974
4975 sx0 = y54*z51 - z54*y51
4976 sy0 = z54*x51 - x54*z51
4977 sz0 = x54*y51 - y54*x51
4978
4979 sx4 = yi5*zi4 - zi5*yi4
4980 sy4 = zi5*xi4 - xi5*zi4
4981 sz4 = xi5*yi4 - yi5*xi4
4982
4983 sx5 = yi4*zi1 - zi4*yi1
4984 sy5 = zi4*xi1 - xi4*zi1
4985 sz5 = xi4*yi1 - yi4*xi1
4986
4987 sx1 = yi1*zi5 - zi1*yi5
4988 sy1 = zi1*xi5 - xi1*zi5
4989 sz1 = xi1*yi5 - yi1*xi5
4990
4991 s2 = one/max(em30,sx0*sx0 + sy0*sy0 + sz0*sz0)
4992 lb0 = (sx0*sx1 + sy0*sy1 + sz0*sz1) * s2
4993 lc0 = (sx0*sx4 + sy0*sy4 + sz0*sz4) * s2
4994 la0 = one-lb0-lc0
4995
4996 IF(lb0 >= zeromt.and. lb0 <= unpt .and.
4997 . lc0 >= zeromt.and. lc0 <= unpt .and.
4998 . la0 >= zeromt.and. la0 <= unpt)THEN
4999 intersect = 1
5000 IF(adt > adt0)THEN
5001 subtria= 4 ! subtriangle
5002 istep = 5
5003 adt0= adt
5004 nx = sx415
5005 ny = sy415
5006 nz = sz415
5007 ENDIF
5008 ENDIF
5009 END IF !((ABS(SX345)+ABS(SY345)+ABS(SZ345))>EM20) THEN
5010 ENDIF
5011 END if!(IXX3 /= IXX4)THEN
5012C
5013 RETURN
5014 END
5015!||====================================================================
5016!|| intersecsh ../engine/source/interfaces/int24/i24dst3.F
5017!||--- called by ------------------------------------------------------
5018!|| i24dst3 ../engine/source/interfaces/int24/i24dst3.F
5019!||--- calls -----------------------------------------------------
5020!|| n2edge3 ../engine/source/interfaces/int24/i24dst3.F
5021!|| n2edge3l ../engine/source/interfaces/int24/n2edge3l.F
5022!||====================================================================
5023 SUBROUTINE intersecsh(
5024 1 ISTEP ,SUBTRIA0,IXX3 ,IXX4 ,
5025 1 INTERSECT,XX1 ,
5026 3 YY1 ,ZZ1 ,XX2 ,YY2 ,ZZ2 ,
5027 4 XX3 ,YY3 ,ZZ3 ,XX4 ,YY4 ,
5028 5 ZZ4 ,XX5 ,YY5 ,ZZ5 ,ADT0 ,
5029 6 XOI ,YOI ,ZOI ,XO1 ,XO2 ,
5030 7 XO3 ,XO4 ,XO5 ,YO1 ,YO2 ,
5031 8 YO3 ,YO4 ,YO5 ,ZO1 ,ZO2 ,
5032 9 ZO3 ,ZO4 ,ZO5 ,NX ,NY ,
5033 A NZ ,SX125 ,SX235 ,SX345 ,SX415 ,
5034 B SY125 ,SY235 ,SY345 ,SY415 ,SZ125 ,
5035 C SZ235 ,SZ345 ,SZ415 ,XI ,YI ,
5036 D ZI ,SOX125 ,SOX235 ,SOX345 ,SOX415 ,
5037 B SOY125 ,SOY235 ,SOY345 ,SOY415 ,SOZ125 ,
5038 C SOZ235 ,SOZ345 ,SOZ415 ,MVOISIN )
5039C-----------------------------------------------
5040C I m p l i c i t T y p e s
5041C-----------------------------------------------
5042#include "implicit_f.inc"
5043C-----------------------------------------------
5044C D u m m y A r g u m e n t s
5045C-----------------------------------------------
5046 INTEGER ISTEP ,SUBTRIA0,IXX3 ,IXX4 ,INTERSECT,MVOISIN(4)
5047C REAL
5048 my_real
5049 1 xx1 ,
5050 3 yy1 ,zz1 ,xx2 ,yy2 ,zz2 ,
5051 4 xx3 ,yy3 ,zz3 ,xx4 ,yy4 ,
5052 5 zz4 ,xx5 ,yy5 ,zz5 ,adt0 ,
5053 6 xoi ,yoi ,zoi ,xo1 ,xo2 ,
5054 7 xo3 ,xo4 ,xo5 ,yo1 ,yo2 ,
5055 8 yo3 ,yo4 ,yo5 ,zo1 ,zo2 ,
5056 9 zo3 ,zo4 ,zo5 ,nx ,ny ,
5057 a nz ,sx125 ,sx235 ,sx345 ,sx415 ,
5058 b sy125 ,sy235 ,sy345 ,sy415 ,sz125 ,
5059 c sz235 ,sz345 ,sz415 ,
5060 d sox125 ,sox235 ,sox345 ,sox415 ,
5061 b soy125 ,soy235 ,soy345 ,soy415 ,soz125 ,
5062 c soz235 ,soz345 ,soz415 ,
5063 d xi,yi,zi
5064C-----------------------------------------------
5065C C o m m o n B l o c k s
5066C-----------------------------------------------
5067#include "com08_c.inc"
5068C-----------------------------------------------
5069C L o c a l V a r i a b l e s
5070C-----------------------------------------------
5071 INTEGER IEG,I,J,ipr,NEG,SUBTRIA
5072 my_real
5073 . x51, x52, x53, x54,
5074 . y51, y52, y53, y54,
5075 . z51, z52, z53, z54,
5076 . xi1, xi2, xi3, xi4, xi5,
5077 . yi1, yi2, yi3, yi4, yi5,
5078 . zi1, zi2, zi3, zi4, zi5,
5079 . xia, xib, xic,
5080 . yia, yib, yic,
5081 . zia, zib, zic,
5082 . xs,ys,zs,
5083 . xm1, xm2, xm3, xm4, xm5,
5084 . ym1, ym2, ym3, ym4, ym5,
5085 . zm1, zm2, zm3, zm4, zm5,la0,lb0,lc0
5086 my_real
5087 . sx,sx1,sx2,sx3,sx4,sx5,sx0,
5088 . sy,sy1,sy2,sy3,sy4,sy5,sy0,
5089 . sz,sz1,sz2,sz3,sz4,sz5,sz0,aaa,adt,s2,
5090 . dno(4),dmin,pij,pik,pij0
5091C------------------------------------------------
5092! 4-------3
5093! |\ /|
5094! | \ 3 / |
5095! | \ /2 |
5096! | 5 |
5097! | 4/ \ |
5098! | / 1 \ |
5099! |/ \|
5100! 1-------2
5101!C-----case secnd shell node/edge shell (12,23,34,41)
5102C--------edge selection: first look at smallest distance DNO,then final intersection
5103 neg = 0
5104 DO j=1,4
5105 IF (mvoisin(j)==0) neg=neg+1
5106 END DO
5107 IF (neg==0) RETURN
5108 dno(1) = (xo1-xoi)*(xo1-xoi)+(yo1-yoi)*(yo1-yoi)+(zo1-zoi)*(zo1-zoi)
5109 IF (mvoisin(1)/=0.AND.mvoisin(4)/=0) dno(1)=ep10
5110 dno(2) = (xo2-xoi)*(xo2-xoi)+(yo2-yoi)*(yo2-yoi)+(zo2-zoi)*(zo2-zoi)
5111 IF (mvoisin(1)/=0.AND.mvoisin(2)/=0) dno(2)=ep10
5112 dno(3) = (xo3-xoi)*(xo3-xoi)+(yo3-yoi)*(yo3-yoi)+(zo3-zoi)*(zo3-zoi)
5113 IF (mvoisin(3)/=0.AND.mvoisin(2)/=0) dno(3)=ep10
5114 IF(ixx3 /= ixx4)THEN
5115 dno(4) = (xo4-xoi)*(xo4-xoi)+(yo4-yoi)*(yo4-yoi)+(zo4-zoi)*(zo4-zoi)
5116 IF (mvoisin(3)/=0.AND.mvoisin(4)/=0) dno(4)=ep10
5117 ELSE
5118 dno(4) = dno(3)
5119 END IF
5120 dmin = ep20
5121 ieg =0
5122 DO j=1,4
5123 IF (dno(j)<dmin) THEN
5124 dmin = dno(j)
5125 ieg = j
5126 END IF
5127 END DO
5128C--------edge selection: smallest pen between IJ,IK
5129 pij = -ep10
5130 pik = zero
5131 SELECT CASE (ieg)
5132 CASE (1) ! 12,41
5133 IF (mvoisin(1)==0) THEN
5134 CALL n2edge3(xo1 ,yo1 ,zo1 ,
5135 2 xo2 ,yo2 ,zo2 ,
5136 3 sox125 ,soy125 ,soz125 ,
5137 4 xoi ,yoi ,zoi ,pij0 )
5138 IF (pij0>zero) THEN
5139 CALL n2edge3l(xx1 ,yy1 ,zz1 ,
5140 2 xx2 ,yy2 ,zz2 ,
5141 3 sx125 ,sy125 ,sz125 ,
5142 4 xi ,yi ,zi ,pij )
5143 IF (pij<zero) THEN
5144 intersect = 1
5145 subtria= 20 + 1
5146 END IF
5147 END IF
5148 END IF !(MVOISIN(1)==0) THEN
5149 IF (mvoisin(4)==0) THEN
5150 CALL n2edge3(xo4 ,yo4 ,zo4 ,
5151 2 xo1 ,yo1 ,zo1 ,
5152 3 sox415 ,soy415 ,soz415 ,
5153 4 xoi ,yoi ,zoi ,pij0 )
5154 IF (pij0>zero) THEN
5155 CALL n2edge3l(xx4 ,yy4 ,zz4 ,
5156 2 xx1 ,yy1 ,zz1 ,
5157 3 sx415 ,sy415 ,sz415 ,
5158 4 xi ,yi ,zi ,pik )
5159 IF (pik<zero.AND.pik>pij) THEN
5160 intersect = 1
5161 subtria= 20 + 4
5162 END IF
5163 END IF
5164 END IF
5165 CASE (2) ! 23,12
5166 IF (mvoisin(2)==0) THEN
5167 CALL n2edge3(xo2 ,yo2 ,zo2 ,
5168 2 xo3 ,yo3 ,zo3 ,
5169 3 sox235 ,soy235 ,soz235 ,
5170 4 xoi ,yoi ,zoi ,pij0 )
5171 IF (pij0>zero) THEN
5172 CALL n2edge3l(xx2 ,yy2 ,zz2 ,
5173 2 xx3 ,yy3 ,zz3 ,
5174 3 sx235 ,sy235 ,sz235 ,
5175 4 xi ,yi ,zi ,pij )
5176 IF (pij<zero) THEN
5177 intersect = 1
5178 subtria= 20 + 2
5179 END IF
5180 END IF
5181 END IF
5182 IF (mvoisin(1)==0) THEN
5183 CALL n2edge3(xo1 ,yo1 ,zo1 ,
5184 2 xo2 ,yo2 ,zo2 ,
5185 3 sox125 ,soy125 ,soz125 ,
5186 4 xoi ,yoi ,zoi ,pij0 )
5187 IF (pij0>zero) THEN
5188 CALL n2edge3l(xx1 ,yy1 ,zz1 ,
5189 2 xx2 ,yy2 ,zz2 ,
5190 3 sx125 ,sy125 ,sz125 ,
5191 4 xi ,yi ,zi ,pik )
5192 IF (pik<zero.AND.pik>pij) THEN
5193 intersect = 1
5194 subtria= 20 + 1
5195 END IF
5196 END IF
5197 END IF
5198 CASE (3) ! 34,23
5199 IF(ixx3 /= ixx4)THEN
5200 IF(mvoisin(3)==0)THEN
5201 CALL n2edge3(xo3 ,yo3 ,zo3 ,
5202 2 xo4 ,yo4 ,zo4 ,
5203 3 sox345 ,soy345 ,soz345 ,
5204 4 xoi ,yoi ,zoi ,pij0)
5205 IF (pij0>zero) THEN
5206 CALL n2edge3l(xx3 ,yy3 ,zz3 ,
5207 2 xx4 ,yy4 ,zz4 ,
5208 3 sx345 ,sy345 ,sz345 ,
5209 4 xi ,yi ,zi ,pij )
5210 IF (pij<zero) THEN
5211 intersect = 1
5212 subtria= 20 + 3
5213 END IF
5214 END IF
5215 END IF
5216 ELSE ! 3N : 41,23
5217 IF (mvoisin(4)==0) THEN
5218 CALL n2edge3(xo4 ,yo4 ,zo4 ,
5219 2 xo1 ,yo1 ,zo1 ,
5220 3 sox415 ,soy415 ,soz415 ,
5221 4 xoi ,yoi ,zoi ,pij0 )
5222 IF (pij0>zero) THEN
5223 CALL n2edge3l(xx4 ,yy4 ,zz4 ,
5224 2 xx1 ,yy1 ,zz1 ,
5225 3 sx415 ,sy415 ,sz415 ,
5226 4 xi ,yi ,zi ,pij )
5227 IF (pij<zero) THEN
5228 intersect = 1
5229 subtria= 20 + 4
5230 END IF
5231 END IF
5232 END IF
5233 END IF
5234C-------------- 23 for both 3N/4N
5235 IF(mvoisin(2)==0)THEN !23
5236 CALL n2edge3(xo2 ,yo2 ,zo2 ,
5237 2 xo3 ,yo3 ,zo3 ,
5238 3 sox235 ,soy235 ,soz235 ,
5239 4 xoi ,yoi ,zoi ,pij0 )
5240 IF (pij0>zero) THEN
5241 CALL n2edge3l(xx2 ,yy2 ,zz2 ,
5242 2 xx3 ,yy3 ,zz3 ,
5243 3 sx235 ,sy235 ,sz235 ,
5244 4 xi ,yi ,zi ,pik )
5245 IF (pik<zero.AND.pik>pij) THEN
5246 intersect = 1
5247 subtria= 20 + 2
5248 END IF
5249 END IF
5250 END IF
5251 CASE (4) ! 41,34
5252 IF(mvoisin(4)==0)THEN
5253 CALL n2edge3(xo4 ,yo4 ,zo4 ,
5254 2 xo1 ,yo1 ,zo1 ,
5255 3 sox415 ,soy415 ,soz415 ,
5256 4 xoi ,yoi ,zoi ,pij0 )
5257 IF (pij0>zero) THEN
5258 CALL n2edge3l(xx4 ,yy4 ,zz4 ,
5259 2 xx1 ,yy1 ,zz1 ,
5260 3 sx415 ,sy415 ,sz415 ,
5261 4 xi ,yi ,zi ,pij )
5262 IF (pij<zero) THEN
5263 intersect = 1
5264 subtria= 20 + 4
5265 END IF
5266 END IF
5267 END IF
5268 IF(ixx3 /= ixx4.AND.mvoisin(3)==0)THEN
5269 CALL n2edge3(xo3 ,yo3 ,zo3 ,
5270 2 xo4 ,yo4 ,zo4 ,
5271 3 sox345 ,soy345 ,soz345 ,
5272 4 xoi ,yoi ,zoi ,pij0)
5273 IF (pij0>zero) THEN
5274 CALL n2edge3l(xx3 ,yy3 ,zz3 ,
5275 2 xx4 ,yy4 ,zz4 ,
5276 3 sx345 ,sy345 ,sz345 ,
5277 4 xi ,yi ,zi ,pik )
5278 IF (pik<zero.AND.pik>pij) THEN
5279 intersect = 1
5280 subtria= 20 + 3
5281 END IF
5282 END IF
5283 END IF
5284 END SELECT
5285C------------out of seg
5286 IF (intersect == 1.AND.pik>zero) intersect = 0
5287 IF(intersect == 1) THEN
5288 subtria0 = subtria
5289 IF(ixx3 == ixx4) subtria0= 20 + 1
5290 istep = 6
5291 adt0 = dt1
5292 END IF
5293C
5294 RETURN
5295 END
5296!||====================================================================
5297!|| n2edge3 ../engine/source/interfaces/int24/i24dst3.F
5298!||--- called by ------------------------------------------------------
5299!|| intersecsh ../engine/source/interfaces/int24/i24dst3.F
5300!||====================================================================
5301 SUBROUTINE n2edge3(
5302 1 XXI ,YYI ,ZZI ,
5303 2 XXJ ,YYJ ,ZZJ ,
5304 3 NX ,NY ,NZ ,
5305 4 XI ,YI ,ZI ,BBB )
5306C-----------------------------------------------
5307C I m p l i c i t T y p e s
5308C-----------------------------------------------
5309#include "implicit_f.inc"
5310C-----------------------------------------------
5311C D u m m y A r g u m e n t s
5312C-----------------------------------------------
5313C REAL
5314 my_real
5315 1 xxi ,yyi ,zzi ,
5316 2 xxj ,yyj ,zzj ,
5317 3 nx ,ny ,nz ,
5318 4 xi ,yi ,zi ,bbb
5319C-----------------------------------------------
5320C L o c a l V a r i a b l e s
5321C-----------------------------------------------
5322 INTEGER I
5323 my_real
5324 . XMIJ(3),NIJ(3),DI(3),AA,NORM
5325C-----------------------------------------------
5326 xmij(1)=xxj-xxi
5327 xmij(2)=yyj-yyi
5328 xmij(3)=zzj-zzi
5329 nij(1)= ny*xmij(3) - nz*xmij(2)
5330 nij(2)= nz*xmij(1) - nx*xmij(3)
5331 nij(3)= nx*xmij(2) - ny*xmij(1)
5332 aa = nij(1)*nij(1)+nij(2)*nij(2)+nij(3)*nij(3)
5333 norm=max(em20,sqrt(aa))
5334 di(1)=half*(xxi+xxj)-xi
5335 di(2)=half*(yyi+yyj)-yi
5336 di(3)=half*(zzi+zzj)-zi
5337 bbb = (di(1)*nij(1)+di(2)*nij(2)+di(3)*nij(3))/norm
5338C
5339 RETURN
5340 END
my_real
#define my_real
Definition cppsort.cpp:32
if
if(complex_arithmetic) id
area
subroutine area(d1, x, x2, y, y2, eint, stif0)
Definition hm_read_prop32.F:567
intersecv0
subroutine intersecv0(istep, subtria, ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, adt0, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, nx, ny, nz, sx125, sx235, sx345, sx415, sy125, sy235, sy345, sy415, sz125, sz235, sz345, sz415, xi, yi, zi, zeromt, unpt)
Definition i24dst3.F:4652
i24nexttria2
subroutine i24nexttria2(izlim, istep, subtria_n, subtria, largep, pene, xxl, yyl, zzl, sx125, sy125, sz125, sx235, sy235, sz235, sx345, sy345, sz345, sx415, sy415, sz415, xi, yi, zi, n1, n2, n3, la, lb, lc, gap2, eps)
Definition i24dst3.F:3751
s_in_m4
subroutine s_in_m4(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4, xi, yi, zi, ier)
Definition i24dst3.F:4048
intersecv
subroutine intersecv(istep, subtria, ixx3, ixx4, intersect, v12, v23, v34, v41, vo12, vo23, vo34, vo41, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, adt0, vxi, vyi, vzi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, nx, ny, nz, sx125, sx235, sx345, sx415, sy125, sy235, sy345, sy415, sz125, sz235, sz345, sz415, xi, yi, zi, zerom, unp, zeromt, unpt)
Definition i24dst3.F:4227
intersecsh
subroutine intersecsh(istep, subtria0, ixx3, ixx4, intersect, xx1, yy1, zz1, xx2, yy2, zz2, xx3, yy3, zz3, xx4, yy4, zz4, xx5, yy5, zz5, adt0, xoi, yoi, zoi, xo1, xo2, xo3, xo4, xo5, yo1, yo2, yo3, yo4, yo5, zo1, zo2, zo3, zo4, zo5, nx, ny, nz, sx125, sx235, sx345, sx415, sy125, sy235, sy345, sy415, sz125, sz235, sz345, sz415, xi, yi, zi, sox125, sox235, sox345, sox415, soy125, soy235, soy345, soy415, soz125, soz235, soz345, soz415, mvoisin)
Definition i24dst3.F:5039
n2edge3
subroutine n2edge3(xxi, yyi, zzi, xxj, yyj, zzj, nx, ny, nz, xi, yi, zi, bbb)
Definition i24dst3.F:5306
i24nexttria
subroutine i24nexttria(izlim, istep, subtria_n, subtria, largep, pene, xxl, yyl, zzl, sx125, sy125, sz125, sx235, sy235, sz235, sx345, sy345, sz345, sx415, sy415, sz415, xi, yi, zi, n1, n2, n3, la, lb, lc, gap2, eps0)
Definition i24dst3.F:3453
i24dst3
subroutine i24dst3(jlt, cand_n, cand_e, cn_loc, ce_loc, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, xi, yi, zi, vx1, vx2, vx3, vx4, vxi, vy1, vy2, vy3, vy4, vyi, vz1, vz2, vz3, vz4, vzi, n1, n2, n3, h1, h2, h3, h4, nin, nsn, ix1, ix2, ix3, ix4, nsvg, stif, jlt_new, inacti, mseglo, gaps, gap_nm, kini, irect, irtlm, time_s, subtria, intth, nsms, pene, xx0, yy0, zz0, vx, vy, vz, ixx, mvoisin, pmax_gap, secnd_fr, gap_m, pene_old, stif_old, itriv, itab, cand_t, iedge, nft, penmin, eps0, nm1, nm2, nm3, intply, dgap_nm, icont_i, marge, nsne, ispt2, izero, iknon, penref)
Definition i24dst3.F:56
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
n2edge3l
subroutine n2edge3l(xxi, yyi, zzi, xxj, yyj, zzj, nx, ny, nz, xi, yi, zi, bbb)
Definition n2edge3l.F:33
tri7box::time_sfi
type(real_pointer), dimension(:), allocatable time_sfi
Definition tri7box.F:542
tri7box::irtlm_fi
type(int_pointer2), dimension(:), allocatable irtlm_fi
Definition tri7box.F:533
tri7box::pene_oldfi
type(real_pointer2), dimension(:), allocatable pene_oldfi
Definition tri7box.F:544
tri7box::icont_i_fi
type(int_pointer), dimension(:), allocatable icont_i_fi
Definition tri7box.F:532
intersect
static2 void intersect(char *uplo, char *diag, Int j, Int start, Int end, Int action, Int *ptrsizebuff, complex **pptrbuff, complex *ptrblock, Int m, Int n, MDESC *ma, Int ia, Int ja, Int templateheight0, Int templatewidth0, MDESC *mb, Int ib, Int jb, Int templateheight1, Int templatewidth1)
Definition pctrmr2.c:176