OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i20lagm.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!|| i20lagm ../engine/source/interfaces/int16/i20lagm.F
25!||--- called by ------------------------------------------------------
26!|| i16main ../engine/source/interfaces/int16/i16main.F
27!||--- calls -----------------------------------------------------
28!|| i20lll ../engine/source/interfaces/int16/i20lagm.F
29!||====================================================================
30 SUBROUTINE i20lagm(X ,V ,LLL ,JLL ,SLL ,
31 2 XLL ,CANDN ,CANDE ,I_STOK,IXS ,
32 3 IXS20 ,IADLL ,EMINX ,NSV ,NELEM ,
33 4 N_MUL_MX,ITASK ,A ,ITIED ,
34 5 NINT ,NKMAX ,COMNTAG)
35C-----------------------------------------------
36C I m p l i c i t T y p e s
37C-----------------------------------------------
38#include "implicit_f.inc"
39C-----------------------------------------------
40C G l o b a l P a r a m e t e r s
41C-----------------------------------------------
42#include "mvsiz_p.inc"
43C-----------------------------------------------
44C C o m m o n B l o c k s
45C-----------------------------------------------
46#include "task_c.inc"
47#include "com04_c.inc"
48#include "com08_c.inc"
49C-----------------------------------------------
50C D u m m y A r g u m e n t s
51C-----------------------------------------------
52 INTEGER I_STOK,N_MUL_MX,ITASK,ITIED,NINT,NKMAX ,
53 . LLL(*),JLL(*),SLL(*),CANDN(*),CANDE(*),COMNTAG(*),
54 . IXS(NIXS,*),IXS20(12,*),IADLL(*),NSV(*) ,NELEM(*)
55C REAL
56 my_real
57 . x(3,*),v(3,*),xll(*),
58 . eminx(6,*),a(3,*)
59C-----------------------------------------------
60C L o c a l V a r i a b l e s
61C-----------------------------------------------
62 INTEGER I,J,K,IK,IE,IS,IC,NK,III(MVSIZ,21),LLT,NFT,LE,FIRST,LAST,
63 . I20
64 my_real
65 . XX(MVSIZ,21),YY(MVSIZ,21),ZZ(MVSIZ,21),
66 . aa,xmin,ymin,zmin,xmax,ymax,zmax,dist
67C-----------------------------------------------
68C
69C
70C | M | Lt| | a | M ao
71C |---+---| | = |
72C | L | 0 | | la | bo
73C
74C [M] a + [L]t la = [M] ao
75C [L] a = bo
76C
77C a = -[M]-1[L]t la + ao
78C [L][M]-1[L]t la = [L] ao - bo
79C
80C on pose:
81C [H] = [L][M]-1[L]t
82C b = [L] ao - bo
83C
84C [H] la = b
85C
86C a = ao - [M]-1[L]t la
87C-----------------------------------------------
88C
89C la : LAMBDA(NC)
90C ao : A(NUMNOD)
91C L : XLL(NK,NC)
92C M : MAS(NUMNOD)
93C [L][M]-1[L]t la : HLA(NC)
94C [L] ao - b : B(NC)
95C [M]-1[L]t la : LTLA(NUMNOD)
96C
97C NC : nombre de contact
98C NK : nombre de noeud pour un contact (8+1,16+1,8+8,16+16)
99C
100C IC : numero du contact (1,NC)
101C IK : numero de noeud local a un contact (1,NK)
102C I : numero global du noeud (1,NUMNOD)
103C
104C IADLL(NC) : IAD = IADLL(IC)
105C LLL(NC*(21,63)) : I = LLL(IAD+1,2...IADNEXT-1)
106C-----------------------------------------------
107C evaluation de b:
108C
109C Vs = Somme(Ni Vi)
110C Vs_ + dt As = Somme(Ni Vi_) + Somme(dt Ni Ai)
111C Somme(dt Ni Ai) - dt As = Vs_ -Somme(Ni Vi_)
112C [L] = dt {N1,N2,..,N15,-1}
113C bo = [L] a = -[L]/dt v_
114C b = [L] ao - bo
115C b = [L] ao + [L]/dt v_ = [L] (v_ + ao dt)/dt
116C-----------------------------------------------
117C b = [L] vo+/dt + vout
118C-----------------------------------------------
119C-----------------------------------------------------------------------
120C boucle sur les candidats au contact
121C-----------------------------------------------------------------------
122 first = 1 + i_stok * itask / nthread
123 last = i_stok*(itask+1) / nthread
124 llt = 0
125 nft=llt+1
126 DO ic=first,last
127 le = cande(ic)
128 ie = nelem(le)
129 i20 = ie - numels8 - numels10
130C-----------------------------------------------------------------------
131C test si brick 20
132C-----------------------------------------------------------------------
133 IF(i20.ge .1.AND.i20.le .numels20)THEN
134 is = nsv(candn(ic))
135 dist = -ep30
136 dist = max(eminx(1,le)-x(1,is)-dt2*(v(1,is)+dt12*a(1,is)),
137 . x(1,is)+dt2*(v(1,is)+dt12*a(1,is))-eminx(4,le),dist)
138 dist = max(eminx(2,le)-x(2,is)-dt2*(v(2,is)+dt12*a(2,is)),
139 . x(2,is)+dt2*(v(2,is)+dt12*a(2,is))-eminx(5,le),dist)
140 dist = max(eminx(3,le)-x(3,is)-dt2*(v(3,is)+dt12*a(3,is)),
141 . x(3,is)+dt2*(v(3,is)+dt12*a(3,is))-eminx(6,le),dist)
142c IF (DIST<0.) CANDN(I) = -CANDN(I)
143C-----------------------------------------------------------------------
144C test si dans la boite
145C-----------------------------------------------------------------------
146 IF(dist.le .zero)THEN
147c
148c print *, "dans la boite",XMIN,XMAX,YMIN,YMAX,ZMIN,ZMAX
149c
150 llt = llt+1
151 iii(llt,21)=is
152 xx(llt,21)=x(1,is)
153 yy(llt,21)=x(2,is)
154 zz(llt,21)=x(3,is)
155 DO k=1,8
156 iii(llt,k)=ixs(k+1,ie)
157 ENDDO
158 DO k=1,12
159 iii(llt,k+8)=ixs20(k,i20)
160 ENDDO
161 DO k=1,20
162 i = iii(llt,k)
163 IF(i/=0)THEN
164 xx(llt,k)=x(1,i)
165 yy(llt,k)=x(2,i)
166 zz(llt,k)=x(3,i)
167 ELSE
168 IF(k==9)THEN
169 xx(llt,k)=0.500*(x(1,iii(llt,1))+x(1,iii(llt,2)))
170 yy(llt,k)=0.500*(x(2,iii(llt,1))+x(2,iii(llt,2)))
171 zz(llt,k)=0.500*(x(3,iii(llt,1))+x(3,iii(llt,2)))
172 ENDIF
173 IF(k==10)THEN
174 xx(llt,k)=0.500*(x(1,iii(llt,2))+x(1,iii(llt,3)))
175 yy(llt,k)=0.500*(x(2,iii(llt,2))+x(2,iii(llt,3)))
176 zz(llt,k)=0.500*(x(3,iii(llt,2))+x(3,iii(llt,3)))
177 ENDIF
178 IF(k==11)THEN
179 xx(llt,k)=0.500*(x(1,iii(llt,3))+x(1,iii(llt,4)))
180 yy(llt,k)=0.500*(x(2,iii(llt,3))+x(2,iii(llt,4)))
181 zz(llt,k)=0.500*(x(3,iii(llt,3))+x(3,iii(llt,4)))
182 ENDIF
183 IF(k==12)THEN
184 xx(llt,k)=0.500*(x(1,iii(llt,4))+x(1,iii(llt,1)))
185 yy(llt,k)=0.500*(x(2,iii(llt,4))+x(2,iii(llt,1)))
186 zz(llt,k)=0.500*(x(3,iii(llt,4))+x(3,iii(llt,1)))
187 ENDIF
188 IF(k==13)THEN
189 xx(llt,k)=0.500*(x(1,iii(llt,1))+x(1,iii(llt,5)))
190 yy(llt,k)=0.500*(x(2,iii(llt,1))+x(2,iii(llt,5)))
191 zz(llt,k)=0.500*(x(3,iii(llt,1))+x(3,iii(llt,5)))
192 ENDIF
193 IF(k==14)THEN
194 xx(llt,k)=0.500*(x(1,iii(llt,2))+x(1,iii(llt,6)))
195 yy(llt,k)=0.500*(x(2,iii(llt,2))+x(2,iii(llt,6)))
196 zz(llt,k)=0.500*(x(3,iii(llt,2))+x(3,iii(llt,6)))
197 ENDIF
198 IF(k==15)THEN
199 xx(llt,k)=0.500*(x(1,iii(llt,3))+x(1,iii(llt,7)))
200 yy(llt,k)=0.500*(x(2,iii(llt,3))+x(2,iii(llt,7)))
201 zz(llt,k)=0.500*(x(3,iii(llt,3))+x(3,iii(llt,7)))
202 ENDIF
203 IF(k==16)THEN
204 xx(llt,k)=0.500*(x(1,iii(llt,4))+x(1,iii(llt,8)))
205 yy(llt,k)=0.500*(x(2,iii(llt,4))+x(2,iii(llt,8)))
206 zz(llt,k)=0.500*(x(3,iii(llt,4))+x(3,iii(llt,8)))
207 ENDIF
208 IF(k==17)THEN
209 xx(llt,k)=0.500*(x(1,iii(llt,5))+x(1,iii(llt,6)))
210 yy(llt,k)=0.500*(x(2,iii(llt,5))+x(2,iii(llt,6)))
211 zz(llt,k)=0.500*(x(3,iii(llt,5))+x(3,iii(llt,6)))
212 ENDIF
213 IF(k==18)THEN
214 xx(llt,k)=0.500*(x(1,iii(llt,6))+x(1,iii(llt,7)))
215 yy(llt,k)=0.500*(x(2,iii(llt,6))+x(2,iii(llt,7)))
216 zz(llt,k)=0.500*(x(3,iii(llt,6))+x(3,iii(llt,7)))
217 ENDIF
218 IF(k==19)THEN
219 xx(llt,k)=0.500*(x(1,iii(llt,7))+x(1,iii(llt,8)))
220 yy(llt,k)=0.500*(x(2,iii(llt,7))+x(2,iii(llt,8)))
221 zz(llt,k)=0.500*(x(3,iii(llt,7))+x(3,iii(llt,8)))
222 ENDIF
223 IF(k==20)THEN
224 xx(llt,k)=0.500*(x(1,iii(llt,5))+x(1,iii(llt,8)))
225 yy(llt,k)=0.500*(x(2,iii(llt,5))+x(2,iii(llt,8)))
226 zz(llt,k)=0.500*(x(3,iii(llt,5))+x(3,iii(llt,8)))
227 ENDIF
228 ENDIF
229 ENDDO
230C-----------------------------------------------------------------------
231C calcul de [L] par paquet de mvsiz
232C-----------------------------------------------------------------------
233 IF(llt==mvsiz-1)THEN
234 CALL i20lll(
235 1 llt ,lll ,jll ,sll ,xll ,v ,
236 2 xx ,yy ,zz ,iii ,iadll ,
237 3 n_mul_mx ,a ,x ,itied ,nint ,nkmax ,
238 4 comntag )
239 nft=llt+1
240 llt = 0
241 ENDIF
242 ELSE
243c debug
244 k=0
245 ENDIF
246 ENDIF
247 ENDDO
248C-----------------------------------------------------------------------
249C calcul de [L] pour dernier paquet
250C-----------------------------------------------------------------------
251 IF(llt/=0) CALL i20lll(
252 1 llt ,lll ,jll ,sll ,xll ,v ,
253 2 xx ,yy ,zz ,iii ,iadll ,
254 3 n_mul_mx ,a ,x ,itied ,nint ,nkmax ,
255 4 comntag )
256C
257C-----------------------------------------------
258 RETURN
259 END
260!||====================================================================
261!|| i20lll ../engine/source/interfaces/int16/i20lagm.F
262!||--- called by ------------------------------------------------------
263!|| i20lagm ../engine/source/interfaces/int16/i20lagm.F
264!||--- calls -----------------------------------------------------
265!|| ancmsg ../engine/source/output/message/message.F
266!|| arret ../engine/source/system/arret.F
267!|| i20rst ../engine/source/interfaces/int16/i20lagm.F
268!||--- uses -----------------------------------------------------
269!|| message_mod ../engine/share/message_module/message_mod.F
270!||====================================================================
271 SUBROUTINE i20lll(LLT ,LLL ,JLL ,SLL ,XLL ,V ,
272 2 XX ,YY ,ZZ ,III ,IADLL ,
273 3 N_MUL_MX,A ,X ,ITIED,NINT ,NKMAX ,
274 4 COMNTAG )
275C-----------------------------------------------
276C M o d u l e s
277C-----------------------------------------------
278 USE message_mod
279C-----------------------------------------------
280C I m p l i c i t T y p e s
281C-----------------------------------------------
282#include "implicit_f.inc"
283#include "comlock.inc"
284C-----------------------------------------------
285C G l o b a l P a r a m e t e r s
286C-----------------------------------------------
287#include "mvsiz_p.inc"
288C-----------------------------------------------
289C C o m m o n B l o c k s
290C-----------------------------------------------
291#include "com08_c.inc"
292 COMMON /lagglob/n_mult
293 INTEGER N_MULT
294C-----------------------------------------------
295C D u m m y A r g u m e n t s
296C-----------------------------------------------
297 INTEGER LLT,N_MUL_MX,ITIED,NINT ,NKMAX
298 INTEGER LLL(*),JLL(*),SLL(*),COMNTAG(*),
299 . iii(mvsiz,21),iadll(*)
300C REAL
301 my_real
302 . xll(*),v(3,*),a(3,*)
303 my_real
304 . xx(mvsiz,21),yy(mvsiz,21),zz(mvsiz,21),x(3,*)
305C-----------------------------------------------
306C L o c a l V a r i a b l e s
307C-----------------------------------------------
308 INTEGER I,J,IK,NK,I1,I2,I3,I4,IAD,NN,J1,J2,IIK,
309 . IPERM1(20),IPERM2(20),COMPTIK
310 my_real
311 . vx,vy,vz,vn,aa
312 my_real
313 . r(mvsiz),s(mvsiz),t(mvsiz),
314 . nsx(mvsiz), nsy(mvsiz), nsz(mvsiz),
315 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
316 . ni(mvsiz,21)
317 DATA iperm1 /0,0,0,0,0,0,0,0,1,2,3,4,5,6,7,8,
318 . 5,6,7,8/
319 DATA iperm2 /0,0,0,0,0,0,0,0,2,3,4,1,1,2,3,4,
320 . 6,7,8,5/
321C-----------------------------------------------
322C calcul de r,s,t
323C-----------------------------------------------
324c
325c print *, "XX(1,1),XX(1,9)",XX(1,1),XX(1,9)
326c
327 CALL i20rst(llt ,r ,s ,t ,ni ,
328 2 nsx ,nsy ,nsz ,nx ,ny ,nz ,
329 3 xx ,yy ,zz )
330C-----------------------------------------------
331C calcul de [L]
332C-----------------------------------------------
333 IF(itied==0)THEN
334 DO i=1,llt
335C-----------------------------------------------
336C test si contact
337C-----------------------------------------------
338 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
339 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
340C
341 nk = 21
342 vx = zero
343 vy = zero
344 vz = zero
345 DO ik=1,8
346 vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
347 vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
348 vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
349 ENDDO
350 DO ik=9,21
351 IF(iii(i,ik)/=0)THEN
352 vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
353 vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
354 vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
355 ELSE
356 vx=vx
357 vy=vy
358 vz=vz
359 ENDIF
360 ENDDO
361c
362c print *, "vx,vy,vz s-m",vx,vy,vz
363c print *, "nx,ny,nz ", NX(I),NY(I),NZ(I)
364c
365 vn = nsx(i)*vx + nsy(i)*vy + nsz(i)*vz
366C-----------------------------------------------
367C test si vitesse entrante en s
368C-----------------------------------------------
369 IF(s(i)*vn<=0.0)THEN
370c
371c print *, "vitesse entrante",vn
372c
373c AA = DT12/SQRT(NX(I)*NX(I)+NY(I)*NY(I)+NZ(I)*NZ(I))
374 aa = one/sqrt(nsx(i)*nsx(i)+nsy(i)*nsy(i)+nsz(i)*nsz(i))
375 nsx(i) = nsx(i)*aa
376 nsy(i) = nsy(i)*aa
377 nsz(i) = nsz(i)*aa
378#include "lockon.inc"
379 n_mult=n_mult+1
380 IF(n_mult>n_mul_mx)THEN
381#include "lockoff.inc"
382 CALL ancmsg(msgid=84,anmode=aninfo)
383 CALL arret(2)
384 ENDIF
385 iad = iadll(n_mult) - 1
386 iik=0
387 comptik=0
388 DO ik=1,21
389 IF(iii(i,ik)/=0)THEN
390 comptik=comptik+1
391 ENDIF
392 ENDDO
393 DO ik=1,21
394 IF(iii(i,ik)/=0)THEN
395 iik=iik+1
396 lll(iad+iik) = iii(i,ik)
397 jll(iad+iik) = 1
398 sll(iad+iik) = 0
399 xll(iad+iik) = nsx(i)*ni(i,ik)
400 lll(iad+comptik+iik) = iii(i,ik)
401 jll(iad+comptik+iik) = 2
402 sll(iad+comptik+iik) = 0
403 xll(iad+comptik+iik) = nsy(i)*ni(i,ik)
404 lll(iad+(2*comptik)+iik) = iii(i,ik)
405 jll(iad+(2*comptik)+iik) = 3
406 sll(iad+(2*comptik)+iik) = 0
407 xll(iad+(2*comptik)+iik) = nsz(i)*ni(i,ik)
408 nn = lll(iad+iik)
409 comntag(nn) = comntag(nn) + 1
410 ELSE
411 j1=iperm1(ik)
412 j2=iperm2(ik)
413 xll(iad+j1)=xll(iad+j1)+0.5*(nsx(i)*ni(i,ik))
414 xll(iad+j2)=xll(iad+j2)+0.5*(nsx(i)*ni(i,ik))
415 xll(iad+comptik+j1)=xll(iad+comptik+j1)+
416 . 0.5*(nsy(i)*ni(i,ik))
417 xll(iad+comptik+j2)=xll(iad+comptik+j2)+
418 . 0.5*(nsy(i)*ni(i,ik))
419 xll(iad+(2*comptik)+j1)=xll(iad+(2*comptik)+j1)+
420 . 0.5*(nsz(i)*ni(i,ik))
421 xll(iad+(2*comptik)+j2)=xll(iad+(2*comptik)+j2)+
422 . 0.5*(nsz(i)*ni(i,ik))
423 ENDIF
424 ENDDO
425 iadll(n_mult+1)=iadll(n_mult)+(3*comptik)
426 IF(iadll(n_mult+1)-1>nkmax)THEN
427#include "lockoff.inc"
428 CALL ancmsg(msgid=84,anmode=aninfo)
429 CALL arret(2)
430 ENDIF
431 sll(iad+comptik) = nint
432 sll(iad+(2*comptik)) = nint
433 sll(iad+(3*comptik)) = nint
434#include "lockoff.inc"
435 ENDIF
436 ENDIF
437 ENDDO
438 ELSEIF(itied==1)THEN
439C-----------------------------------------------
440C ITIED = 1
441C-----------------------------------------------
442 DO i=1,llt
443C-----------------------------------------------
444C test si contact
445C-----------------------------------------------
446 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
447 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
448C
449 nk = 21
450 vx = zero
451 vy = zero
452 vz = zero
453 DO ik=1,8
454 vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
455 vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
456 vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
457 ENDDO
458 DO ik=9,21
459 IF(iii(i,ik)/=0)THEN
460 vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
461 vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
462 vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
463 ELSE
464 vx=vx
465 vy=vy
466 vz=vz
467 ENDIF
468 ENDDO
469c
470c print *, "vx,vy,vz s-m",vx,vy,vz
471c print *, "nx,ny,nz ", NX(I),NY(I),NZ(I)
472c
473 vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
474C-----------------------------------------------
475C test si vitesse entrante en r,s ou t
476C-----------------------------------------------
477 IF(vn<=zero)THEN
478c
479c print *, "vitesse entrante",vn
480c
481#include "lockon.inc"
482 IF(n_mult+3>n_mul_mx)THEN
483#include "lockoff.inc"
484 CALL ancmsg(msgid=84,anmode=aninfo)
485 CALL arret(2)
486 ENDIF
487 IF(iadll(n_mult+1)-1+21*3>nkmax)THEN
488#include "lockoff.inc"
489 CALL ancmsg(msgid=84,anmode=aninfo)
490 CALL arret(2)
491 ENDIF
492C
493 n_mult=n_mult+1
494 iad = iadll(n_mult) - 1
495 iik=0
496 DO ik=1,21
497 IF(iii(i,ik)/=0)THEN
498 iik=iik+1
499 lll(iad+iik) = iii(i,ik)
500 jll(iad+iik) = 1
501 sll(iad+iik) = 0
502 xll(iad+iik) = ni(i,ik)
503 nn = lll(iad+iik)
504 comntag(nn) = comntag(nn) + 1
505 ELSE
506 j1=iperm1(ik)
507 j2=iperm2(ik)
508 xll(iad+j1)=xll(iad+j1)+0.5*ni(i,ik)
509 xll(iad+j2)=xll(iad+j2)+0.5*ni(i,ik)
510 ENDIF
511 ENDDO
512 sll(iad+iik) = nint
513 iadll(n_mult+1)=iadll(n_mult) + iik
514C
515 n_mult=n_mult+1
516 iad = iadll(n_mult) - 1
517 iik=0
518 DO ik=1,21
519 IF(iii(i,ik)/=0)THEN
520 iik=iik+1
521 lll(iad+iik) = iii(i,ik)
522 jll(iad+iik) = 2
523 sll(iad+iik) = 0
524 xll(iad+iik) = ni(i,ik)
525 nn = lll(iad+iik)
526 comntag(nn) = comntag(nn) + 1
527 ELSE
528 j1=iperm1(ik)
529 j2=iperm2(ik)
530 xll(iad+j1)=xll(iad+j1)+0.5*ni(i,ik)
531 xll(iad+j2)=xll(iad+j2)+0.5*ni(i,ik)
532 ENDIF
533 ENDDO
534 iadll(n_mult+1)=iadll(n_mult) + iik
535 sll(iad+iik) = nint
536C
537 n_mult=n_mult+1
538 iad = iadll(n_mult) - 1
539 iik=0
540 DO ik=1,21
541 IF(iii(i,ik)/=0)THEN
542 iik=iik+1
543 lll(iad+iik) = iii(i,ik)
544 jll(iad+iik) = 3
545 sll(iad+iik) = 0
546 xll(iad+iik) = ni(i,ik)
547 nn = lll(iad+iik)
548 comntag(nn) = comntag(nn) + 1
549 ELSE
550 j1=iperm1(ik)
551 j2=iperm2(ik)
552 xll(iad+j1)=xll(iad+j1)+0.5*ni(i,ik)
553 xll(iad+j2)=xll(iad+j2)+0.5*ni(i,ik)
554 ENDIF
555 ENDDO
556 iadll(n_mult+1)=iadll(n_mult) + iik
557 sll(iad+iik) = nint
558#include "lockoff.inc"
559 ENDIF
560 ENDIF
561 ENDDO
562 ELSE
563C-----------------------------------------------
564C ITIED = 2
565C-----------------------------------------------
566 DO i=1,llt
567C-----------------------------------------------
568C test si contact
569C-----------------------------------------------
570 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
571 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
572C
573 nk = 21
574C-----------------------------------------------
575c
576#include "lockon.inc"
577 IF(n_mult+3>n_mul_mx)THEN
578#include "lockoff.inc"
579 CALL ancmsg(msgid=84,anmode=aninfo)
580 CALL arret(2)
581 ENDIF
582 IF(iadll(n_mult+1)-1+21*3>nkmax)THEN
583#include "lockoff.inc"
584 CALL ancmsg(msgid=84,anmode=aninfo)
585 CALL arret(2)
586 ENDIF
587 n_mult=n_mult+1
588 iad = iadll(n_mult) - 1
589 iik=0
590 DO ik=1,21
591 IF(iii(i,ik)/=0)THEN
592 iik=iik+1
593 lll(iad+iik) = iii(i,ik)
594 jll(iad+iik) = 1
595 sll(iad+iik) = 0
596 xll(iad+iik) = ni(i,ik)
597 nn = lll(iad+iik)
598 comntag(nn) = comntag(nn) + 1
599 ELSE
600 j1=iperm1(ik)
601 j2=iperm2(ik)
602 xll(iad+j1)=xll(iad+j1)+0.5*ni(i,ik)
603 xll(iad+j2)=xll(iad+j2)+0.5*ni(i,ik)
604 ENDIF
605 ENDDO
606 iadll(n_mult+1)=iadll(n_mult) + iik
607 sll(iad+iik) = nint
608C
609 n_mult=n_mult+1
610 iad = iadll(n_mult) - 1
611 iik=0
612 DO ik=1,21
613 IF(iii(i,ik)/=0)THEN
614 iik=iik+1
615 lll(iad+iik) = iii(i,ik)
616 jll(iad+iik) = 2
617 sll(iad+iik) = 0
618 xll(iad+iik) = ni(i,ik)
619 nn = lll(iad+iik)
620 comntag(nn) = comntag(nn) + 1
621 ELSE
622 j1=iperm1(ik)
623 j2=iperm2(ik)
624 xll(iad+j1)=xll(iad+j1)+0.5*ni(i,ik)
625 xll(iad+j2)=xll(iad+j2)+0.5*ni(i,ik)
626 ENDIF
627 ENDDO
628 iadll(n_mult+1)=iadll(n_mult) + iik
629 sll(iad+iik) = nint
630C
631 n_mult=n_mult+1
632 iad = iadll(n_mult) - 1
633 iik=0
634 DO ik=1,21
635 IF(iii(i,ik)/=0)THEN
636 iik=iik+1
637 lll(iad+iik) = iii(i,ik)
638 jll(iad+iik) = 3
639 sll(iad+iik) = 0
640 xll(iad+iik) = ni(i,ik)
641 nn = lll(iad+iik)
642 comntag(nn) = comntag(nn) + 1
643 ELSE
644 j1=iperm1(ik)
645 j2=iperm2(ik)
646 xll(iad+j1)=xll(iad+j1)+0.5*ni(i,ik)
647 xll(iad+j2)=xll(iad+j2)+0.5*ni(i,ik)
648 ENDIF
649 ENDDO
650 iadll(n_mult+1)=iadll(n_mult) + iik
651 sll(iad+iik) = nint
652C
653#include "lockoff.inc"
654 ENDIF
655 ENDDO
656 ENDIF
657c
658c print *, "r,s,t",r(1),s(1),t(1)
659C
660 RETURN
661 END
662C
663!||====================================================================
664!|| i20rst ../engine/source/interfaces/int16/i20lagm.F
665!||--- called by ------------------------------------------------------
666!|| i20lll ../engine/source/interfaces/int16/i20lagm.F
667!|| i21lll ../engine/source/interfaces/int17/i21lagm.F
668!||--- calls -----------------------------------------------------
669!|| i20deri ../engine/source/interfaces/int16/i20lagm.F
670!|| i20ni ../engine/source/interfaces/int16/i20lagm.F
671!|| i20rstn ../engine/source/interfaces/int16/i20lagm.F
672!||====================================================================
673 SUBROUTINE i20rst(LLT,R ,S ,T ,NI ,
674 2 NSX ,NSY ,NSZ ,NX ,NY ,NZ ,
675 3 XX ,YY ,ZZ )
676C-----------------------------------------------
677C I m p l i c i t T y p e s
678C-----------------------------------------------
679#include "implicit_f.inc"
680C-----------------------------------------------
681C G l o b a l P a r a m e t e r s
682C-----------------------------------------------
683#include "mvsiz_p.inc"
684C-----------------------------------------------
685C D u m m y A r g u m e n t s
686C-----------------------------------------------
687 INTEGER LLT
688C REAL
689 my_real
690 . XX(MVSIZ,21),YY(MVSIZ,21),ZZ(MVSIZ,21)
691 my_real
692 . R(MVSIZ),S(MVSIZ),T(MVSIZ),NI(MVSIZ,21) ,
693 . NSX(MVSIZ),NSY(MVSIZ),NSZ(MVSIZ),
694 . NX(MVSIZ),NY(MVSIZ),NZ(MVSIZ)
695C-----------------------------------------------
696C L o c a l V a r i a b l e s
697C-----------------------------------------------
698 INTEGER I,J,IK,NK,ITER,NITERMAX,JTER,NJTERMAX,CONV
699 my_real
700 . vx,vy,vz,vn
701 my_real
702 . sn, rn, tn,
703 . drdx(mvsiz),drdy(mvsiz),drdz(mvsiz),
704 . dsdx(mvsiz),dsdy(mvsiz),dsdz(mvsiz),
705 . dtdx(mvsiz),dtdy(mvsiz),dtdz(mvsiz),
706 . dxdr(mvsiz),dydr(mvsiz),dzdr(mvsiz),
707 . dxds(mvsiz),dyds(mvsiz),dzds(mvsiz),
708 . dxdt(mvsiz),dydt(mvsiz),dzdt(mvsiz),
709 . rr(mvsiz),ss(mvsiz),tt(mvsiz)
710C-----------------------------------------------
711C
712C r=s=t=0
713C
714C +---> iter
715C |
716C | Ni(r,s,t) =
717C | dNi/dr =
718C | ... _
719C | \
720C | dx/dr = /_ (xi * dNi/dr)
721C | ...
722C |
723C | [dx/dr dy/dr dz/dr]
724C | [J] = |dx/ds dy/ds dz/ds|
725C | [dx/dt dy/dt dz/dt]
726C |
727C | +--> jter
728C | | _
729C | | \
730C | | x(r,s,t) = /_ (xi * Ni(r,s,t))
731C | | ...
732C | |
733C | | |r| |r| -1 |xs-x(r,s,t)|
734C | | {s} = {s} + [J] {ys-y(r,s,t)}
735C | | |t| |t| |zs-z(r,s,t)|
736C | |
737C | | Ni(r,s,t) =
738C +-+---
739C-----------------------------------------------
740 nitermax = 3
741 njtermax = 3
742 conv = 0
743C
744 DO i=1,llt
745 rr(i) = zero
746 ss(i) = zero
747 tt(i) = zero
748 ENDDO
749C-----------------------------------------------
750C calcul de r,s,t et Ni(r,s,t)
751C-----------------------------------------------
752 DO iter=1,nitermax
753c
754c print *, "iter",iter
755c
756C-----------------------------------------------
757C calcul de Ni(r,s,t); [J]; [J]-1
758C-----------------------------------------------
759 CALL i20deri(llt,rr ,ss ,tt ,ni ,
760 2 drdx ,drdy ,drdz ,dsdx ,dsdy ,dsdz ,
761 3 dtdx ,dtdy ,dtdz ,dxdr ,dydr ,dzdr ,
762 4 dxds ,dyds ,dzds ,dxdt ,dydt ,dzdt ,
763 5 xx ,yy ,zz )
764C
765 DO jter=1,njtermax
766c
767c print *, "jter",jter
768c
769C-----------------------------------------------
770C calcul de r,s,t new
771C-----------------------------------------------
772 CALL i20rstn(llt,rr,ss ,tt ,ni ,conv ,
773 2 drdx ,drdy ,drdz ,dsdx ,dsdy ,dsdz ,
774 3 dtdx ,dtdy ,dtdz ,xx ,yy ,zz ,
775 4 r ,s ,t )
776c
777c print *, "r,s,t",r(1),s(1),t(1)
778c print *, "rr,ss,tt",rr(1),ss(1),tt(1)
779c
780C-----------------------------------------------
781C calcul de Ni(-1<r<1 , -1<s<1 , -1<t<1)
782C-----------------------------------------------
783 CALL i20ni(llt,rr ,ss ,tt ,ni )
784C pb de parith on si conv fonction de mvsiz !!!!!!!
785C IF(CONV/=0)RETURN
786C
787 ENDDO
788 ENDDO
789C
790 DO i=1,llt
791 nsx(i) = dydt(i)*dzdr(i) - dzdt(i)*dydr(i)
792 nsy(i) = dzdt(i)*dxdr(i) - dxdt(i)*dzdr(i)
793 nsz(i) = dxdt(i)*dydr(i) - dydt(i)*dxdr(i)
794C
795 sn = ss(i) * ss(i)
796 sn = sn * sn
797 sn = sn * sn
798 sn = sn * sn
799 sn = sn * sn
800 sn = sn * ss(i)
801 nx(i) = sn * nsx(i)
802 ny(i) = sn * nsy(i)
803 nz(i) = sn * nsz(i)
804C
805 rn = rr(i) * rr(i)
806 rn = rn * rn
807 rn = rn * rn
808 rn = rn * rn
809 rn = rn * rn
810 rn = rn * rr(i)
811 nx(i) = nx(i) + rn * (dyds(i)*dzdt(i) - dzds(i)*dydt(i))
812 ny(i) = ny(i) + rn * (dzds(i)*dxdt(i) - dxds(i)*dzdt(i))
813 nz(i) = nz(i) + rn * (dxds(i)*dydt(i) - dyds(i)*dxdt(i))
814C
815 tn = tt(i) * tt(i)
816 tn = tn * tn
817 tn = tn * tn
818 tn = tn * tn
819 tn = tn * tn
820 tn = tn * tt(i)
821 nx(i) = nx(i) + tn * (dydr(i)*dzds(i) - dzdr(i)*dyds(i))
822 ny(i) = ny(i) + tn * (dzdr(i)*dxds(i) - dxdr(i)*dzds(i))
823 nz(i) = nz(i) + tn * (dxdr(i)*dyds(i) - dydr(i)*dxds(i))
824C
825 ENDDO
826C
827 RETURN
828 END
829C
830!||====================================================================
831!|| i20ni ../engine/source/interfaces/int16/i20lagm.F
832!||--- called by ------------------------------------------------------
833!|| i20rst ../engine/source/interfaces/int16/i20lagm.F
834!||====================================================================
835 SUBROUTINE i20ni(LLT,RR ,SS ,TT ,NI )
836C-----------------------------------------------
837C I m p l i c i t T y p e s
838C-----------------------------------------------
839#include "implicit_f.inc"
840C-----------------------------------------------
841C G l o b a l P a r a m e t e r s
842C-----------------------------------------------
843#include "mvsiz_p.inc"
844C-----------------------------------------------
845C D u m m y A r g u m e n t s
846C-----------------------------------------------
847 INTEGER LLT
848 my_real
849 . RR(MVSIZ),SS(MVSIZ),TT(MVSIZ),NI(MVSIZ,21)
850C-----------------------------------------------
851C L o c a l V a r i a b l e s
852C-----------------------------------------------
853 INTEGER I
854 my_real
855 . U_M_R,U_P_R,U_M_S,U_P_S,U_M_T,U_P_T,
856 . ums_umt,ums_upt,ups_umt,ups_upt,
857 . umr_ums,umr_ups,upr_ums,upr_ups,
858 . umt_umr,umt_upr,upt_umr,upt_upr,
859 . a,r05,s05,t05
860C-----------------------------------------------------------------------
861C calcul de Ni
862C-----------------------------------------------------------------------
863 DO i=1,llt
864C
865 r05 = half*rr(i)
866 s05 = half*ss(i)
867 t05 = half*tt(i)
868C
869 u_m_r = half - r05
870 u_p_r = half + r05
871C
872 u_m_s = half - s05
873 u_p_s = half + s05
874C
875 u_m_t = half - t05
876 u_p_t = half + t05
877C
878 ums_umt = u_m_s * u_m_t
879 ums_upt = u_m_s * u_p_t
880 ups_umt = u_p_s * u_m_t
881 ups_upt = u_p_s * u_p_t
882C
883 umr_ums = u_m_r * u_m_s
884 umr_ups = u_m_r * u_p_s
885 upr_ums = u_p_r * u_m_s
886 upr_ups = u_p_r * u_p_s
887C
888 umt_umr = u_m_t * u_m_r
889 umt_upr = u_m_t * u_p_r
890 upt_umr = u_p_t * u_m_r
891 upt_upr = u_p_t * u_p_r
892C
893 ni(i,1) = u_m_r * ums_umt * (-rr(i)-ss(i)-tt(i)-two)
894 ni(i,2) = u_m_r * ums_upt * (-rr(i)-ss(i)+tt(i)-two)
895 ni(i,3) = u_p_r * ums_upt * ( rr(i)-ss(i)+tt(i)-two)
896 ni(i,4) = u_p_r * ums_umt * ( rr(i)-ss(i)-tt(i)-two)
897 ni(i,5) = u_m_r * ups_umt * (-rr(i)+ss(i)-tt(i)-two)
898 ni(i,6) = u_m_r * ups_upt * (-rr(i)+ss(i)+tt(i)-two)
899 ni(i,7) = u_p_r * ups_upt * ( rr(i)+ss(i)+tt(i)-two)
900 ni(i,8) = u_p_r * ups_umt * ( rr(i)+ss(i)-tt(i)-two)
901C------------------------------------
902 a = (one-rr(i)*rr(i))
903 ni(i,10) = a * ums_upt
904 ni(i,12) = a * ums_umt
905 ni(i,18) = a * ups_upt
906 ni(i,20) = a * ups_umt
907C------------------------------------
908 a = (one-ss(i)*ss(i))
909 ni(i,13) = a * umt_umr
910 ni(i,14) = a * upt_umr
911 ni(i,15) = a * upt_upr
912 ni(i,16) = a * umt_upr
913C------------------------------------
914 a = (one-tt(i)*tt(i))
915 ni(i,9) = a * umr_ums
916 ni(i,11) = a * upr_ums
917 ni(i,17) = a * umr_ups
918 ni(i,19) = a * upr_ups
919C------------------------------------
920 ni(i,21) = -one
921C------------------------------------
922 ENDDO
923C-----------------------------------------------
924 RETURN
925 END
926!||====================================================================
927!|| i20rstn ../engine/source/interfaces/int16/i20lagm.F
928!||--- called by ------------------------------------------------------
929!|| i20rst ../engine/source/interfaces/int16/i20lagm.F
930!||====================================================================
931 SUBROUTINE i20rstn(LLT,RR,SS ,TT ,NI ,CONV ,
932 2 DRDX ,DRDY ,DRDZ ,DSDX ,DSDY ,DSDZ ,
933 3 DTDX ,DTDY ,DTDZ ,XX ,YY ,ZZ ,
934 4 R ,S ,T )
935C-----------------------------------------------
936C I m p l i c i t T y p e s
937C-----------------------------------------------
938c#include "implicit_f.inc"
939 implicit none
940C-----------------------------------------------
941C G l o b a l P a r a m e t e r s
942C-----------------------------------------------
943#include "mvsiz_p.inc"
944#include "constant.inc"
945C-----------------------------------------------
946C D u m m y A r g u m e n t s
947C-----------------------------------------------
948 INTEGER LLT,CONV
949 my_real
950 . R(MVSIZ),S(MVSIZ),T(MVSIZ),NI(MVSIZ,21) ,
951 . rr(mvsiz),ss(mvsiz),tt(mvsiz),
952 . xx(mvsiz,21) ,yy(mvsiz,21) ,zz(mvsiz,21) ,
953 . drdx(mvsiz),drdy(mvsiz),drdz(mvsiz),
954 . dsdx(mvsiz),dsdy(mvsiz),dsdz(mvsiz),
955 . dtdx(mvsiz),dtdy(mvsiz),dtdz(mvsiz)
956C-----------------------------------------------
957C L o c a l V a r i a b l e s
958C-----------------------------------------------
959 INTEGER I
960 my_real
961 . dx ,dy,dz,dr ,ds,dt,err
962C
963 err=zero
964C-----------------------------------------------
965 DO i=1,llt
966C
967 dx = xx(i,21)
968 + - ni(i, 1)*xx(i, 1) - ni(i, 2)*xx(i, 2) - ni(i, 3)*xx(i, 3)
969 + - ni(i, 4)*xx(i, 4) - ni(i, 5)*xx(i, 5) - ni(i, 6)*xx(i, 6)
970 + - ni(i, 7)*xx(i, 7) - ni(i, 8)*xx(i, 8) - ni(i, 9)*xx(i, 9)
971 + - ni(i,10)*xx(i,10) - ni(i,11)*xx(i,11) - ni(i,12)*xx(i,12)
972 + - ni(i,13)*xx(i,13) - ni(i,14)*xx(i,14) - ni(i,15)*xx(i,15)
973 + - ni(i,16)*xx(i,16) - ni(i,17)*xx(i,17) - ni(i,18)*xx(i,18)
974 + - ni(i,19)*xx(i,19) - ni(i,20)*xx(i,20)
975 dy = yy(i,21)
976 + - ni(i, 1)*yy(i, 1) - ni(i, 2)*yy(i, 2) - ni(i, 3)*yy(i, 3)
977 + - ni(i, 4)*yy(i, 4) - ni(i, 5)*yy(i, 5) - ni(i, 6)*yy(i, 6)
978 + - ni(i, 7)*yy(i, 7) - ni(i, 8)*yy(i, 8) - ni(i, 9)*yy(i, 9)
979 + - ni(i,10)*yy(i,10) - ni(i,11)*yy(i,11) - ni(i,12)*yy(i,12)
980 + - ni(i,13)*yy(i,13) - ni(i,14)*yy(i,14) - ni(i,15)*yy(i,15)
981 + - ni(i,16)*yy(i,16) - ni(i,17)*yy(i,17) - ni(i,18)*yy(i,18)
982 + - ni(i,19)*yy(i,19) - ni(i,20)*yy(i,20)
983 dz = zz(i,21)
984 + - ni(i, 1)*zz(i, 1) - ni(i, 2)*zz(i, 2) - ni(i, 3)*zz(i, 3)
985 + - ni(i, 4)*zz(i, 4) - ni(i, 5)*zz(i, 5) - ni(i, 6)*zz(i, 6)
986 + - ni(i, 7)*zz(i, 7) - ni(i, 8)*zz(i, 8) - ni(i, 9)*zz(i, 9)
987 + - ni(i,10)*zz(i,10) - ni(i,11)*zz(i,11) - ni(i,12)*zz(i,12)
988 + - ni(i,13)*zz(i,13) - ni(i,14)*zz(i,14) - ni(i,15)*zz(i,15)
989 + - ni(i,16)*zz(i,16) - ni(i,17)*zz(i,17) - ni(i,18)*zz(i,18)
990 + - ni(i,19)*zz(i,19) - ni(i,20)*zz(i,20)
991C
992 dr = drdx(i)*dx + drdy(i)*dy + drdz(i)*dz
993 ds = dsdx(i)*dx + dsdy(i)*dy + dsdz(i)*dz
994 dt = dtdx(i)*dx + dtdy(i)*dy + dtdz(i)*dz
995C
996c
997c print *, "DRDX(I),DRDY(I),DRDZ(I)",DRDX(I),DRDY(I),DRDZ(I)
998c print *, "DSDX(I),DSDY(I),DSDZ(I)",DSDX(I),DSDY(I),DSDZ(I)
999c print *, "DTDX(I),DTDY(I),DTDZ(I)",DTDX(I),DTDY(I),DTDZ(I)
1000c print *, "Ni",ni(1,1),ni(1,2),ni(1,3),ni(1,4),ni(1,5),ni(1,9)
1001c print *, "dx,dy,dz",dx ,dy ,dz
1002c
1003 rr(i) = rr(i) + dr
1004 ss(i) = ss(i) + ds
1005 tt(i) = tt(i) + dt
1006C
1007 r(i) = rr(i)
1008 s(i) = ss(i)
1009 t(i) = tt(i)
1010C
1011 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
1012 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
1013 err = max(err,abs(dr),abs(ds),abs(dt))
1014 ELSE
1015 rr(i) = max(min(rr(i),one),-one)
1016 ss(i) = max(min(ss(i),one),-one)
1017 tt(i) = max(min(tt(i),one),-one)
1018 ENDIF
1019c
1020c print *, "3r,s,t",r(1),s(1),t(1)
1021c print *, "3rr,ss,tt",rr(1),ss(1),tt(1)
1022c print *, "dr,ds,dt",dr ,ds ,dt
1023c print *, "r,s,t",r(1),s(1),t(1)
1024c print *, "ERR",ERR
1025c
1026C
1027 ENDDO
1028C
1029 IF(err<em4) conv = 1
1030C-----------------------------------------------
1031 RETURN
1032 END
1033!||====================================================================
1034!|| i20deri ../engine/source/interfaces/int16/i20lagm.F
1035!||--- called by ------------------------------------------------------
1036!|| i20rst ../engine/source/interfaces/int16/i20lagm.F
1037!||====================================================================
1038 SUBROUTINE i20deri(LLT,RR,SS ,TT ,NI ,
1039 2 DRDX ,DRDY ,DRDZ ,DSDX ,DSDY ,DSDZ ,
1040 3 DTDX ,DTDY ,DTDZ ,DXDR ,DYDR ,DZDR ,
1041 4 DXDS ,DYDS ,DZDS ,DXDT ,DYDT ,DZDT ,
1042 5 XX ,YY ,ZZ )
1043C-----------------------------------------------
1044C I m p l i c i t T y p e s
1045C-----------------------------------------------
1046#include "implicit_f.inc"
1047C-----------------------------------------------
1048C G l o b a l P a r a m e t e r s
1049C-----------------------------------------------
1050#include "mvsiz_p.inc"
1051C-----------------------------------------------
1052C D u m m y A r g u m e n t s
1053C-----------------------------------------------
1054 INTEGER LLT
1055C REAL
1056 my_real
1057 . DXDR(MVSIZ), DYDR(MVSIZ), DZDR(MVSIZ),
1058 . DXDS(MVSIZ), DYDS(MVSIZ), DZDS(MVSIZ),
1059 . DXDT(MVSIZ), DYDT(MVSIZ), DZDT(MVSIZ),
1060 . DRDX(MVSIZ), DSDX(MVSIZ), DTDX(MVSIZ),
1061 . DRDY(MVSIZ), DSDY(MVSIZ), DTDY(MVSIZ),
1062 . DRDZ(MVSIZ), DSDZ(MVSIZ), DTDZ(MVSIZ),
1063 . XX(MVSIZ,21) ,YY(MVSIZ,21),ZZ(MVSIZ,21),
1064 . NI(MVSIZ,21) ,RR(MVSIZ) ,SS(MVSIZ) ,TT(MVSIZ)
1065C-----------------------------------------------
1066C L o c a l V a r i a b l e s
1067C-----------------------------------------------
1068 INTEGER I,N
1069 my_real
1070 . dnidr(20),dnids(20),dnidt(20),
1071 . d, aa, bb, det(mvsiz),r9 ,r13 ,s9 ,s10 ,s11 ,s12 ,t10 ,t14
1072 my_real
1073 . u_m_r,u_p_r,u_m_s,u_p_s,u_m_t,u_p_t,
1074 . ums_umt,ums_upt,ups_umt,ups_upt,
1075 . umr_ums,umr_ups,upr_ums,upr_ups,
1076 . umt_umr,umt_upr,upt_umr,upt_upr,
1077 . a,r05,s05,t05
1078C-----------------------------------------------
1079C
1080C
1081C ^ S _ T
1082C | /|
1083C | /
1084C | /
1085C ( 7)==========|===(18)===============( 6)
1086C //| | / //|
1087C // | | / //||
1088C // | | // ||
1089C // | | / // ||
1090C // | | / // ||
1091C (19) | *- - - + - * - - - - -* (17) ||
1092C // | /| | /| / /| // ||
1093C // | / // ||
1094C // (15)/ | |/ | + / // (14)
1095C // |#- - - - - # - -/- - -# // ||
1096C // | * - - /|- -*- - -/ -//* ||
1097C ( 8)===============(20)==============( 5) ||
1098C || / || / | / | / | / ||| | ||
1099C || @- | - - - @ - - - - -@ || ||
1100C R <-----||- - -+ -|- -# - - -| - # - - -|- -#|| - -+ ||
1101C || | * - - /| - *- - -/-|| * ||
1102C || | || / | / | |||/ ||
1103C || ( 3)-------/--|---(10)----||---------( 2)
1104C || @ /- / - - @ - -/ - - @ || //
1105C || |/ #- - -/| - # - - -|- -#|| //
1106C (16) / / + /| /(13) //
1107C || /| | | || //
1108C || / / / | / || //
1109C || (11) @- - - - - @ - + - - -@ || ( 9)
1110C || / || //
1111C || / || //
1112C || / || //
1113C || / ||//
1114C ||/ ||/
1115C ( 4)===============(12)==============( 1)
1116C
1117C
1118C
1119C
1120C
1121C*/
1122C-----------------------------------------------
1123C _
1124C \
1125C x(r,s,t) = /_ (xi * Ni(r,s,t))
1126C _
1127C \
1128C y(r,s,t) = /_ (yi * Ni(r,s,t))
1129C _
1130C \
1131C z(r,s,t) = /_ zi * Ni(r,s,t))
1132C
1133C _
1134C \
1135C dx/dr = /_ (xi * dNi/dr)
1136C ...
1137C
1138C [dx/dr dy/dr dz/dr]
1139C [J] = |dx/ds dy/ds dz/ds|
1140C [dx/dt dy/dt dz/dt]
1141C
1142C |r| |r| -1 |xs-x|
1143C {s} = {s} + [J] {ys-y}
1144C |t| |t| |zs-z|
1145C
1146C-----------------------------------------------------------------------
1147C Ni; dNi/dr; dNi/ds; dNi/dt
1148C-----------------------------------------------------------------------
1149 DO i=1,llt
1150 r05 = half*rr(i)
1151 s05 = half*ss(i)
1152 t05 = half*tt(i)
1153C
1154 u_m_r = half - r05
1155 u_p_r = half + r05
1156C
1157 u_m_s = half - s05
1158 u_p_s = half + s05
1159C
1160 u_m_t = half - t05
1161 u_p_t = half + t05
1162C
1163 ums_umt = u_m_s * u_m_t
1164 ums_upt = u_m_s * u_p_t
1165 ups_umt = u_p_s * u_m_t
1166 ups_upt = u_p_s * u_p_t
1167C
1168 umr_ums = u_m_r * u_m_s
1169 umr_ups = u_m_r * u_p_s
1170 upr_ums = u_p_r * u_m_s
1171 upr_ups = u_p_r * u_p_s
1172C
1173 umt_umr = u_m_t * u_m_r
1174 umt_upr = u_m_t * u_p_r
1175 upt_umr = u_p_t * u_m_r
1176 upt_upr = u_p_t * u_p_r
1177C
1178C
1179C
1180C
1181 ni(i,1) = u_m_r * ums_umt * (-rr(i)-ss(i)-tt(i)-two)
1182 ni(i,2) = u_m_r * ums_upt * (-rr(i)-ss(i)+tt(i)-two)
1183 ni(i,3) = u_p_r * ums_upt * ( rr(i)-ss(i)+tt(i)-two)
1184 ni(i,4) = u_p_r * ums_umt * ( rr(i)-ss(i)-tt(i)-two)
1185 ni(i,5) = u_m_r * ups_umt * (-rr(i)+ss(i)-tt(i)-two)
1186 ni(i,6) = u_m_r * ups_upt * (-rr(i)+ss(i)+tt(i)-two)
1187 ni(i,7) = u_p_r * ups_upt * ( rr(i)+ss(i)+tt(i)-two)
1188 ni(i,8) = u_p_r * ups_umt * ( rr(i)+ss(i)-tt(i)-two)
1189C
1190 dnidr(1) = -ums_umt * (u_m_s + u_m_t - rr(i) -three_half)
1191 dnidr(2) = -ums_upt * (u_m_s + u_p_t - rr(i) -three_half)
1192 dnidr(3) = ums_upt * (u_m_s + u_p_t + rr(i) -three_half)
1193 dnidr(4) = ums_umt * (u_m_s + u_m_t + rr(i) -three_half)
1194 dnidr(5) = -ups_umt * (u_p_s + u_m_t - rr(i) -three_half)
1195 dnidr(6) = -ups_upt * (u_p_s + u_p_t - rr(i) -three_half)
1196 dnidr(7) = ups_upt * (u_p_s + u_p_t + rr(i) -three_half)
1197 dnidr(8) = ups_umt * (u_p_s + u_m_t + rr(i) -three_half)
1198C
1199 dnids(1) = -umt_umr * (u_m_r + u_m_t - ss(i) -three_half)
1200 dnids(2) = -upt_umr * (u_m_r + u_p_t - ss(i) -three_half)
1201 dnids(3) = -upt_upr * (u_p_r + u_p_t - ss(i) -three_half)
1202 dnids(4) = -umt_upr * (u_p_r + u_m_t - ss(i) -three_half)
1203 dnids(5) = umt_umr * (u_m_r + u_m_t + ss(i) -three_half)
1204 dnids(6) = upt_umr * (u_m_r + u_p_t + ss(i) -three_half)
1205 dnids(7) = upt_upr * (u_p_r + u_p_t + ss(i) -three_half)
1206 dnids(8) = umt_upr * (u_p_r + u_m_t + ss(i) -three_half)
1207C
1208 dnidt(1) = -umr_ums * (u_m_r + u_m_s - tt(i) -three_half)
1209 dnidt(2) = umr_ums * (u_m_r + u_m_s + tt(i) -three_half)
1210 dnidt(3) = upr_ums * (u_p_r + u_m_s + tt(i) -three_half)
1211 dnidt(4) = -upr_ums * (u_p_r + u_m_s - tt(i) -three_half)
1212 dnidt(5) = -umr_ups * (u_m_r + u_p_s - tt(i) -three_half)
1213 dnidt(6) = umr_ups * (u_m_r + u_p_s + tt(i) -three_half)
1214 dnidt(7) = upr_ups * (u_p_r + u_p_s + tt(i) -three_half)
1215 dnidt(8) = -upr_ups * (u_p_r + u_p_s - tt(i) -three_half)
1216C------------------------------------
1217 a = (one-rr(i)*rr(i))
1218 ni(i,10) = a * ums_upt
1219 ni(i,12) = a * ums_umt
1220 ni(i,18) = a * ups_upt
1221 ni(i,20) = a * ups_umt
1222C
1223 a = half*a
1224 dnidt(10) = a * u_m_s
1225 dnidt(18) = a * u_p_s
1226 dnidt(12) = -dnidt(10)
1227 dnidt(20) = -dnidt(18)
1228C
1229 dnids(18) = a * u_p_t
1230 dnids(20) = a * u_m_t
1231 dnids(10) = -dnids(18)
1232 dnids(12) = -dnids(20)
1233C
1234 a = -two*rr(i)
1235 dnidr(10) = a * ums_upt
1236 dnidr(12) = a * ums_umt
1237 dnidr(18) = a * ups_upt
1238 dnidr(20) = a * ups_umt
1239C------------------------------------
1240 a = (one-ss(i)*ss(i))
1241 ni(i,13) = a * umt_umr
1242 ni(i,14) = a * upt_umr
1243 ni(i,15) = a * upt_upr
1244 ni(i,16) = a * umt_upr
1245C
1246 a = half*a
1247 dnidr(15) = a * u_p_t
1248 dnidr(16) = a * u_m_t
1249 dnidr(13) = -dnidr(16)
1250 dnidr(14) = -dnidr(15)
1251C
1252 dnidt(14) = a * u_m_r
1253 dnidt(15) = a * u_p_r
1254 dnidt(13) = -dnidt(14)
1255 dnidt(16) = -dnidt(15)
1256C
1257 a = -two*ss(i)
1258 dnids(13) = a * umt_umr
1259 dnids(14) = a * upt_umr
1260 dnids(15) = a * upt_upr
1261 dnids(16) = a * umt_upr
1262C------------------------------------
1263 a = (one-tt(i)*tt(i))
1264 ni(i,9) = a * umr_ums
1265 ni(i,11) = a * upr_ums
1266 ni(i,17) = a * umr_ups
1267 ni(i,19) = a * upr_ups
1268C
1269 ni(i,21) = -one
1270C
1271 a = half*a
1272 dnidr(11) = a * u_m_s
1273 dnidr(19) = a * u_p_s
1274 dnidr(9) = -dnidr(11)
1275 dnidr(17) = -dnidr(19)
1276C
1277 dnids(17) = a * u_m_r
1278 dnids(19) = a * u_p_r
1279 dnids(9) = -dnids(17)
1280 dnids(11) = -dnids(19)
1281C
1282 a = -two*tt(i)
1283 dnidt(9) = a * umr_ums
1284 dnidt(11) = a * upr_ums
1285 dnidt(17) = a * umr_ups
1286 dnidt(19) = a * upr_ups
1287C
1288c
1289c print *, "DNIDr(1),DNIDr(9)",DNIDr(1),DNIDr(9)
1290c print *, "DNIDs(1),DNIDs(9)",DNIDs(1),DNIDs(9)
1291c print *, "DNIDT(1),DNIDT(9)",DNIDT(1),DNIDT(9)
1292c print *, "XX(I,1),XX(I,9)",XX(I,1),XX(I,9)
1293c
1294C-----------------------------------------------------------------------
1295C dx/dr; dx/ds; dx/dt
1296C-----------------------------------------------------------------------
1297 dxdr(i) = dnidr(1)*xx(i,1) + dnidr(2)*xx(i,2) + dnidr(3)*xx(i,3)
1298 + + dnidr(4)*xx(i,4) + dnidr(5)*xx(i,5) + dnidr(6)*xx(i,6)
1299 + + dnidr(7)*xx(i,7) + dnidr(8)*xx(i,8)
1300 + + dnidr(9)*xx(i,9) + dnidr(10)*xx(i,10) + dnidr(11)*xx(i,11)
1301 + + dnidr(12)*xx(i,12) + dnidr(13)*xx(i,13) + dnidr(14)*xx(i,14)
1302 + + dnidr(15)*xx(i,15) + dnidr(16)*xx(i,16) + dnidr(17)*xx(i,17)
1303 + + dnidr(18)*xx(i,18) + dnidr(19)*xx(i,19) + dnidr(20)*xx(i,20)
1304C
1305 dxds(i) = dnids(1)*xx(i,1) + dnids(2)*xx(i,2) + dnids(3)*xx(i,3)
1306 + + dnids(4)*xx(i,4) + dnids(5)*xx(i,5) + dnids(6)*xx(i,6)
1307 + + dnids(7)*xx(i,7) + dnids(8)*xx(i,8)
1308 + + dnids(9)*xx(i,9) + dnids(10)*xx(i,10) + dnids(11)*xx(i,11)
1309 + + dnids(12)*xx(i,12) + dnids(13)*xx(i,13) + dnids(14)*xx(i,14)
1310 + + dnids(15)*xx(i,15) + dnids(16)*xx(i,16) + dnids(17)*xx(i,17)
1311 + + dnids(18)*xx(i,18) + dnids(19)*xx(i,19) + dnids(20)*xx(i,20)
1312C
1313 dxdt(i) = dnidt(1)*xx(i,1) + dnidt(2)*xx(i,2) + dnidt(3)*xx(i,3)
1314 + + dnidt(4)*xx(i,4) + dnidt(5)*xx(i,5) + dnidt(6)*xx(i,6)
1315 + + dnidt(7)*xx(i,7) + dnidt(8)*xx(i,8)
1316 + + dnidt(9)*xx(i,9) + dnidt(10)*xx(i,10) + dnidt(11)*xx(i,11)
1317 + + dnidt(12)*xx(i,12) + dnidt(13)*xx(i,13) + dnidt(14)*xx(i,14)
1318 + + dnidt(15)*xx(i,15) + dnidt(16)*xx(i,16) + dnidt(17)*xx(i,17)
1319 + + dnidt(18)*xx(i,18) + dnidt(19)*xx(i,19) + dnidt(20)*xx(i,20)
1320C-----------------------------------------------------------------------
1321C dy/dr; dy/ds; dy/dt
1322C-----------------------------------------------------------------------
1323 dydr(i) = dnidr(1)*yy(i,1) + dnidr(2)*yy(i,2) + dnidr(3)*yy(i,3)
1324 + + dnidr(4)*yy(i,4) + dnidr(5)*yy(i,5) + dnidr(6)*yy(i,6)
1325 + + dnidr(7)*yy(i,7) + dnidr(8)*yy(i,8)
1326 + + dnidr(9)*yy(i,9) + dnidr(10)*yy(i,10) + dnidr(11)*yy(i,11)
1327 + + dnidr(12)*yy(i,12) + dnidr(13)*yy(i,13) + dnidr(14)*yy(i,14)
1328 + + dnidr(15)*yy(i,15) + dnidr(16)*yy(i,16) + dnidr(17)*yy(i,17)
1329 + + dnidr(18)*yy(i,18) + dnidr(19)*yy(i,19) + dnidr(20)*yy(i,20)
1330C
1331 dyds(i) = dnids(1)*yy(i,1) + dnids(2)*yy(i,2) + dnids(3)*yy(i,3)
1332 + + dnids(4)*yy(i,4) + dnids(5)*yy(i,5) + dnids(6)*yy(i,6)
1333 + + dnids(7)*yy(i,7) + dnids(8)*yy(i,8)
1334 + + dnids(9)*yy(i,9) + dnids(10)*yy(i,10) + dnids(11)*yy(i,11)
1335 + + dnids(12)*yy(i,12) + dnids(13)*yy(i,13) + dnids(14)*yy(i,14)
1336 + + dnids(15)*yy(i,15) + dnids(16)*yy(i,16) + dnids(17)*yy(i,17)
1337 + + dnids(18)*yy(i,18) + dnids(19)*yy(i,19) + dnids(20)*yy(i,20)
1338C
1339 dydt(i) = dnidt(1)*yy(i,1) + dnidt(2)*yy(i,2) + dnidt(3)*yy(i,3)
1340 + + dnidt(4)*yy(i,4) + dnidt(5)*yy(i,5) + dnidt(6)*yy(i,6)
1341 + + dnidt(7)*yy(i,7) + dnidt(8)*yy(i,8)
1342 + + dnidt(9)*yy(i,9) + dnidt(10)*yy(i,10) + dnidt(11)*yy(i,11)
1343 + + dnidt(12)*yy(i,12) + dnidt(13)*yy(i,13) + dnidt(14)*yy(i,14)
1344 + + dnidt(15)*yy(i,15) + dnidt(16)*yy(i,16) + dnidt(17)*yy(i,17)
1345 + + dnidt(18)*yy(i,18) + dnidt(19)*yy(i,19) + dnidt(20)*yy(i,20)
1346C-----------------------------------------------------------------------
1347C dz/dr; dz/ds; dz/dt
1348C-----------------------------------------------------------------------
1349 dzdr(i) = dnidr(1)*zz(i,1) + dnidr(2)*zz(i,2) + dnidr(3)*zz(i,3)
1350 + + dnidr(4)*zz(i,4) + dnidr(5)*zz(i,5) + dnidr(6)*zz(i,6)
1351 + + dnidr(7)*zz(i,7) + dnidr(8)*zz(i,8)
1352 + + dnidr(9)*zz(i,9) + dnidr(10)*zz(i,10) + dnidr(11)*zz(i,11)
1353 + + dnidr(12)*zz(i,12) + dnidr(13)*zz(i,13) + dnidr(14)*zz(i,14)
1354 + + dnidr(15)*zz(i,15) + dnidr(16)*zz(i,16) + dnidr(17)*zz(i,17)
1355 + + dnidr(18)*zz(i,18) + dnidr(19)*zz(i,19) + dnidr(20)*zz(i,20)
1356C
1357 dzds(i) = dnids(1)*zz(i,1) + dnids(2)*zz(i,2) + dnids(3)*zz(i,3)
1358 + + dnids(4)*zz(i,4) + dnids(5)*zz(i,5) + dnids(6)*zz(i,6)
1359 + + dnids(7)*zz(i,7) + dnids(8)*zz(i,8)
1360 + + dnids(9)*zz(i,9) + dnids(10)*zz(i,10) + dnids(11)*zz(i,11)
1361 + + dnids(12)*zz(i,12) + dnids(13)*zz(i,13) + dnids(14)*zz(i,14)
1362 + + dnids(15)*zz(i,15) + dnids(16)*zz(i,16) + dnids(17)*zz(i,17)
1363 + + dnids(18)*zz(i,18) + dnids(19)*zz(i,19) + dnids(20)*zz(i,20)
1364C
1365 dzdt(i) = dnidt(1)*zz(i,1) + dnidt(2)*zz(i,2) + dnidt(3)*zz(i,3)
1366 + + dnidt(4)*zz(i,4) + dnidt(5)*zz(i,5) + dnidt(6)*zz(i,6)
1367 + + dnidt(7)*zz(i,7) + dnidt(8)*zz(i,8)
1368 + + dnidt(9)*zz(i,9) + dnidt(10)*zz(i,10) + dnidt(11)*zz(i,11)
1369 + + dnidt(12)*zz(i,12) + dnidt(13)*zz(i,13) + dnidt(14)*zz(i,14)
1370 + + dnidt(15)*zz(i,15) + dnidt(16)*zz(i,16) + dnidt(17)*zz(i,17)
1371 + + dnidt(18)*zz(i,18) + dnidt(19)*zz(i,19) + dnidt(20)*zz(i,20)
1372C-----------------------------------------------------------------------
1373C -1
1374C [J] Inversion du jacobien
1375C-----------------------------------------------------------------------
1376 drdx(i)=dyds(i)*dzdt(i)-dzds(i)*dydt(i)
1377 drdy(i)=dzds(i)*dxdt(i)-dxds(i)*dzdt(i)
1378 drdz(i)=dxds(i)*dydt(i)-dyds(i)*dxdt(i)
1379C
1380 dsdz(i)=dxdt(i)*dydr(i)-dydt(i)*dxdr(i)
1381 dsdy(i)=dzdt(i)*dxdr(i)-dxdt(i)*dzdr(i)
1382 dsdx(i)=dydt(i)*dzdr(i)-dzdt(i)*dydr(i)
1383C
1384 dtdx(i)=dydr(i)*dzds(i)-dzdr(i)*dyds(i)
1385 dtdy(i)=dzdr(i)*dxds(i)-dxdr(i)*dzds(i)
1386 dtdz(i)=dxdr(i)*dyds(i)-dydr(i)*dxds(i)
1387C
1388 det(i) = dxdr(i) * drdx(i)
1389 . + dydr(i) * drdy(i)
1390 . + dzdr(i) * drdz(i)
1391C
1392c
1393c print *, "Det",DET(I)
1394c print *, "DXDR(I),DYDR(I),DZDR(I)",DXDR(I),DYDR(I),DZDR(I)
1395c print *, "DXDs(I),DYDs(I),DZDs(I)",DXDs(I),DYDs(I),DZDs(I)
1396c print *, "DXDt(I),DYDt(I),DZDt(I)",DXDt(I),DYDt(I),DZDt(I)
1397c
1398 ENDDO
1399C
1400 DO i=1,llt
1401C-----------------------------------------------------------------------
1402C -1
1403C [J] Inversion du jacobien suite
1404C-----------------------------------------------------------------------
1405 d = one/det(i)
1406 drdx(i)=d*drdx(i)
1407 dsdx(i)=d*dsdx(i)
1408 dtdx(i)=d*dtdx(i)
1409C
1410 drdy(i)=d*drdy(i)
1411 dsdy(i)=d*dsdy(i)
1412 dtdy(i)=d*dtdy(i)
1413C
1414 drdz(i)=d*drdz(i)
1415 dsdz(i)=d*dsdz(i)
1416 dtdz(i)=d*dtdz(i)
1417C
1418c
1419c print *, "DRDX(I),DRDY(I),DRDZ(I)",DRDX(I),DRDY(I),DRDZ(I)
1420c print *, "DSDX(I),DSDY(I),DSDZ(I)",DSDX(I),DSDY(I),DSDZ(I)
1421c print *, "DTDX(I),DTDY(I),DTDZ(I)",DTDX(I),DTDY(I),DTDZ(I)
1422c
1423 ENDDO
1424C-----------------------------------------------------------------------
1425 RETURN
1426 END
my_real
#define my_real
Definition cppsort.cpp:32
i20deri
subroutine i20deri(llt, rr, ss, tt, ni, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, dxdr, dydr, dzdr, dxds, dyds, dzds, dxdt, dydt, dzdt, xx, yy, zz)
Definition i20lagm.F:1043
i20lagm
subroutine i20lagm(x, v, lll, jll, sll, xll, candn, cande, i_stok, ixs, ixs20, iadll, eminx, nsv, nelem, n_mul_mx, itask, a, itied, nint, nkmax, comntag)
Definition i20lagm.F:35
i20rstn
subroutine i20rstn(llt, rr, ss, tt, ni, conv, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, xx, yy, zz, r, s, t)
Definition i20lagm.F:935
i20ni
subroutine i20ni(llt, rr, ss, tt, ni)
Definition i20lagm.F:836
i20rst
subroutine i20rst(llt, r, s, t, ni, nsx, nsy, nsz, nx, ny, nz, xx, yy, zz)
Definition i20lagm.F:676
i20lll
subroutine i20lll(llt, lll, jll, sll, xll, v, xx, yy, zz, iii, iadll, n_mul_mx, a, x, itied, nint, nkmax, comntag)
Definition i20lagm.F:275
ymax
subroutine ymax(idn, fac, npc, pld, stiffmin, stiffmax, stiffini, stiffavg)
Definition law100_upd.F:272
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
message_mod
Definition message_mod.F:1249
ancmsg
subroutine ancmsg(msgid, msgtype, anmode, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10, i11, i12, i13, i14, i15, i16, i17, i18, i19, i20, r1, r2, r3, r4, r5, r6, r7, r8, r9, c1, c2, c3, c4, c5, c6, c7, c8, c9, prmode)
Definition message.F:889
arret
subroutine arret(nn)
Definition arret.F:87