OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i16lagm.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!|| i16lagm ../engine/source/interfaces/int16/i16lagm.F
25!||--- called by ------------------------------------------------------
26!|| i16main ../engine/source/interfaces/int16/i16main.F
27!||--- calls -----------------------------------------------------
28!|| i16lll ../engine/source/interfaces/int16/i16lagm.F
29!||====================================================================
30 SUBROUTINE i16lagm(X ,V ,LLL ,JLL ,SLL ,
31 2 XLL ,CANDN ,CANDE ,I_STOK,IXS ,
32 3 IXS16 ,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,*),IXS16(8,*),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,17),LLT,NFT,LE,FIRST,LAST,
63 . I16
64 my_real
65 . XX(MVSIZ,17),YY(MVSIZ,17),ZZ(MVSIZ,17),
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*(17,51)) : 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 i16 = ie - numels8 - numels10 - numels20
130C-----------------------------------------------------------------------
131C test si shell 16
132C-----------------------------------------------------------------------
133 IF(i16.ge .1.AND.i16.le .numels16)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,17)=is
152c XX(LLT,17)=X(1,IS)+DT2*(V(1,IS)+DT12*A(1,IS))
153c YY(LLT,17)=X(2,IS)+DT2*(V(2,IS)+DT12*A(2,IS))
154c ZZ(LLT,17)=X(3,IS)+DT2*(V(3,IS)+DT12*A(3,IS))
155 xx(llt,17)=x(1,is)
156 yy(llt,17)=x(2,is)
157 zz(llt,17)=x(3,is)
158 DO k=1,8
159 iii(llt,k)=ixs(k+1,ie)
160 iii(llt,k+8)=ixs16(k,i16)
161 ENDDO
162 DO k=1,16
163 i = iii(llt,k)
164c XX(LLT,K)=X(1,I)+DT2*(V(1,I)+DT12*A(1,I))
165c YY(LLT,K)=X(2,I)+DT2*(V(2,I)+DT12*A(2,I))
166c ZZ(LLT,K)=X(3,I)+DT2*(V(3,I)+DT12*A(3,I))
167 xx(llt,k)=x(1,i)
168 yy(llt,k)=x(2,i)
169 zz(llt,k)=x(3,i)
170 ENDDO
171c
172c print *, "XX(1,1),XX(1,9)",XX(1,1),XX(1,9)
173c
174C-----------------------------------------------------------------------
175C calcul de [L] par paquet de mvsiz
176C-----------------------------------------------------------------------
177 IF(llt==mvsiz-1)THEN
178 CALL i16lll(
179 1 llt ,lll ,jll ,sll ,xll ,v ,
180 2 xx ,yy ,zz ,iii ,iadll ,
181 3 n_mul_mx ,a ,x ,itied ,nint ,nkmax ,
182 4 comntag )
183 nft=llt+1
184 llt = 0
185 ENDIF
186 ELSE
187c debug
188 k=0
189 ENDIF
190 ENDIF
191 ENDDO
192C-----------------------------------------------------------------------
193C calcul de [L] pour dernier paquet
194C-----------------------------------------------------------------------
195 IF(llt/=0) CALL i16lll(
196 1 llt ,lll ,jll ,sll ,xll ,v ,
197 2 xx ,yy ,zz ,iii ,iadll ,
198 3 n_mul_mx ,a ,x ,itied ,nint ,nkmax ,
199 4 comntag )
200C
201C-----------------------------------------------
202 RETURN
203 END
204!||====================================================================
205!|| i16lll ../engine/source/interfaces/int16/i16lagm.F
206!||--- called by ------------------------------------------------------
207!|| i16lagm ../engine/source/interfaces/int16/i16lagm.F
208!||--- calls -----------------------------------------------------
209!|| ancmsg ../engine/source/output/message/message.F
210!|| arret ../engine/source/system/arret.F
211!|| i16rst ../engine/source/interfaces/int16/i16lagm.F
212!||--- uses -----------------------------------------------------
213!|| message_mod ../engine/share/message_module/message_mod.F
214!||====================================================================
215 SUBROUTINE i16lll(LLT ,LLL ,JLL ,SLL ,XLL ,V ,
216 2 XX ,YY ,ZZ ,III ,IADLL ,
217 3 N_MUL_MX,A ,X ,ITIED,NINT ,NKMAX ,
218 4 COMNTAG )
219C-----------------------------------------------
220C M o d u l e s
221C-----------------------------------------------
222 USE message_mod
223C-----------------------------------------------
224C I m p l i c i t T y p e s
225C-----------------------------------------------
226#include "implicit_f.inc"
227#include "comlock.inc"
228C-----------------------------------------------
229C G l o b a l P a r a m e t e r s
230C-----------------------------------------------
231#include "mvsiz_p.inc"
232C-----------------------------------------------
233C C o m m o n B l o c k s
234C-----------------------------------------------
235#include "com08_c.inc"
236 COMMON /lagglob/n_mult
237 INTEGER N_MULT
238C-----------------------------------------------
239C D u m m y A r g u m e n t s
240C-----------------------------------------------
241 INTEGER LLT,N_MUL_MX,ITIED,NINT ,NKMAX
242 INTEGER LLL(*),JLL(*),SLL(*),COMNTAG(*),
243 . iii(mvsiz,17),iadll(*)
244C REAL
245 my_real
246 . xll(*),v(3,*),a(3,*)
247 my_real
248 . xx(mvsiz,17),yy(mvsiz,17),zz(mvsiz,17),x(3,*)
249C-----------------------------------------------
250C L o c a l V a r i a b l e s
251C-----------------------------------------------
252 INTEGER I,J,IK,NK,I1,I2,I3,I4,IAD,NN
253 my_real
254 . vx,vy,vz,vn,aa
255 my_real
256 . r(mvsiz),s(mvsiz),t(mvsiz),
257 . nx(mvsiz), ny(mvsiz), nz(mvsiz),
258 . ni(mvsiz,17)
259C-----------------------------------------------
260C calcul de r,s,t
261C-----------------------------------------------
262c
263c print *, "XX(1,1),XX(1,9)",XX(1,1),XX(1,9)
264c
265 CALL i16rst(llt ,r ,s ,t ,ni ,
266 2 nx ,ny ,nz ,xx ,yy ,zz )
267C-----------------------------------------------
268C calcul de [L]
269C-----------------------------------------------
270 IF(itied==0)THEN
271 DO i=1,llt
272C-----------------------------------------------
273C test si contact
274C-----------------------------------------------
275 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
276 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
277C
278 nk = 17
279 vx = zero
280 vy = zero
281 vz = zero
282 DO ik=1,nk
283 vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
284 vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
285 vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
286 ENDDO
287c
288c print *, "vx,vy,vz s-m",vx,vy,vz
289c print *, "nx,ny,nz ", NX(I),NY(I),NZ(I)
290c
291 vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
292C-----------------------------------------------
293C test si vitesse entrante en s
294C-----------------------------------------------
295 IF(s(i)*vn<=zero)THEN
296c
297c print *, "vitesse entrante",vn
298c print *, "s = ",S(I)
299c
300c AA = DT12/SQRT(NX(I)*NX(I)+NY(I)*NY(I)+NZ(I)*NZ(I))
301 aa = one/sqrt(nx(i)*nx(i)+ny(i)*ny(i)+nz(i)*nz(i))
302 nx(i) = nx(i)*aa
303 ny(i) = ny(i)*aa
304 nz(i) = nz(i)*aa
305#include "lockon.inc"
306 n_mult=n_mult+1
307 IF(n_mult>n_mul_mx)THEN
308#include "lockoff.inc"
309 CALL ancmsg(msgid=84,anmode=aninfo)
310 CALL arret(2)
311 ENDIF
312 iadll(n_mult+1)=iadll(n_mult) + 51
313 IF(iadll(n_mult+1)-1>nkmax)THEN
314#include "lockoff.inc"
315 CALL ancmsg(msgid=84,anmode=aninfo)
316 CALL arret(2)
317 ENDIF
318 iad = iadll(n_mult) - 1
319 DO ik=1,17
320 lll(iad+ik) = iii(i,ik)
321 jll(iad+ik) = 1
322 sll(iad+ik) = 0
323 xll(iad+ik) = nx(i)*ni(i,ik)
324 lll(iad+ik+17) = iii(i,ik)
325 jll(iad+ik+17) = 2
326 sll(iad+ik+17) = 0
327 xll(iad+ik+17) = ny(i)*ni(i,ik)
328 lll(iad+ik+34) = iii(i,ik)
329 jll(iad+ik+34) = 3
330 sll(iad+ik+34) = 0
331 xll(iad+ik+34) = nz(i)*ni(i,ik)
332 nn = lll(iad+ik)
333 comntag(nn) = comntag(nn) + 1
334 ENDDO
335 sll(iad+17) = nint
336 sll(iad+34) = nint
337 sll(iad+51) = nint
338#include "lockoff.inc"
339 ENDIF
340 ENDIF
341 ENDDO
342 ELSEIF(itied==1)THEN
343C-----------------------------------------------
344C ITIED = 1
345C-----------------------------------------------
346 DO i=1,llt
347C-----------------------------------------------
348C test si contact
349C-----------------------------------------------
350 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
351 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
352C
353 nk = 17
354 vx = zero
355 vy = zero
356 vz = zero
357 DO ik=1,nk
358 vx = vx - (v(1,iii(i,ik))+dt12*a(1,iii(i,ik)))*ni(i,ik)
359 vy = vy - (v(2,iii(i,ik))+dt12*a(2,iii(i,ik)))*ni(i,ik)
360 vz = vz - (v(3,iii(i,ik))+dt12*a(3,iii(i,ik)))*ni(i,ik)
361 ENDDO
362c
363c print *, "vx,vy,vz s-m",vx,vy,vz
364c print *, "nx,ny,nz ", NX(I),NY(I),NZ(I)
365c
366 vn = nx(i)*vx + ny(i)*vy + nz(i)*vz
367C-----------------------------------------------
368C test si vitesse entrante en s
369C-----------------------------------------------
370 IF(s(i)*vn<=zero)THEN
371c
372c print *, "vitesse entrante",vn
373c print *, "s = ",S(I)
374c
375#include "lockon.inc"
376 IF(n_mult+3>n_mul_mx)THEN
377#include "lockoff.inc"
378 CALL ancmsg(msgid=84,anmode=aninfo)
379 CALL arret(2)
380 ENDIF
381 IF(iadll(n_mult+1)-1+17*3>nkmax)THEN
382#include "lockoff.inc"
383 CALL ancmsg(msgid=84,anmode=aninfo)
384 CALL arret(2)
385 ENDIF
386C
387 n_mult=n_mult+1
388 iadll(n_mult+1)=iadll(n_mult) + 17
389 iad = iadll(n_mult) - 1
390 DO ik=1,17
391 lll(iad+ik) = iii(i,ik)
392 jll(iad+ik) = 1
393 sll(iad+ik) = 0
394 xll(iad+ik) = ni(i,ik)
395 nn = lll(iad+ik)
396 comntag(nn) = comntag(nn) + 1
397 ENDDO
398 sll(iad+17) = nint
399C
400 n_mult=n_mult+1
401 iadll(n_mult+1)=iadll(n_mult) + 17
402 iad = iadll(n_mult) - 1
403 DO ik=1,17
404 lll(iad+ik) = iii(i,ik)
405 jll(iad+ik) = 2
406 sll(iad+ik) = 0
407 xll(iad+ik) = ni(i,ik)
408 nn = lll(iad+ik)
409 comntag(nn) = comntag(nn) + 1
410 ENDDO
411 sll(iad+17) = nint
412C
413 n_mult=n_mult+1
414 iadll(n_mult+1)=iadll(n_mult) + 17
415 iad = iadll(n_mult) - 1
416 DO ik=1,17
417 lll(iad+ik) = iii(i,ik)
418 jll(iad+ik) = 3
419 sll(iad+ik) = 0
420 xll(iad+ik) = ni(i,ik)
421 nn = lll(iad+ik)
422 comntag(nn) = comntag(nn) + 1
423 ENDDO
424 sll(iad+17) = nint
425#include "lockoff.inc"
426 ENDIF
427 ENDIF
428 ENDDO
429 ELSE
430C-----------------------------------------------
431C ITIED = 2
432C-----------------------------------------------
433 DO i=1,llt
434C-----------------------------------------------
435C test si contact
436C-----------------------------------------------
437 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
438 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
439C
440 nk = 17
441C-----------------------------------------------
442c print *, "s = ",S(I)
443c
444#include "lockon.inc"
445 IF(n_mult+3>n_mul_mx)THEN
446#include "lockoff.inc"
447 CALL ancmsg(msgid=84,anmode=aninfo)
448 CALL arret(2)
449 ENDIF
450 IF(iadll(n_mult+1)-1+17*3>nkmax)THEN
451#include "lockoff.inc"
452 CALL ancmsg(msgid=84,anmode=aninfo)
453 CALL arret(2)
454 ENDIF
455 n_mult=n_mult+1
456 iadll(n_mult+1)=iadll(n_mult) + 17
457 iad = iadll(n_mult) - 1
458 DO ik=1,17
459 lll(iad+ik) = iii(i,ik)
460 jll(iad+ik) = 1
461 sll(iad+ik) = 0
462 xll(iad+ik) = ni(i,ik)
463 nn = lll(iad+ik)
464 comntag(nn) = comntag(nn) + 1
465 ENDDO
466 sll(iad+17) = nint
467C
468 n_mult=n_mult+1
469 iadll(n_mult+1)=iadll(n_mult) + 17
470 iad = iadll(n_mult) - 1
471 DO ik=1,17
472 lll(iad+ik) = iii(i,ik)
473 jll(iad+ik) = 2
474 sll(iad+ik) = 0
475 xll(iad+ik) = ni(i,ik)
476 nn = lll(iad+ik)
477 comntag(nn) = comntag(nn) + 1
478 ENDDO
479 sll(iad+17) = nint
480C
481 n_mult=n_mult+1
482 iadll(n_mult+1)=iadll(n_mult) + 17
483 iad = iadll(n_mult) - 1
484 DO ik=1,17
485 lll(iad+ik) = iii(i,ik)
486 jll(iad+ik) = 3
487 sll(iad+ik) = 0
488 xll(iad+ik) = ni(i,ik)
489 nn = lll(iad+ik)
490 comntag(nn) = comntag(nn) + 1
491 ENDDO
492 sll(iad+17) = nint
493C
494#include "lockoff.inc"
495 ENDIF
496 ENDDO
497 ENDIF
498c
499c print *, "r,s,t",r(1),s(1),t(1)
500C
501 RETURN
502 END
503C
504!||====================================================================
505!|| i16rst ../engine/source/interfaces/int16/i16lagm.F
506!||--- called by ------------------------------------------------------
507!|| i16lll ../engine/source/interfaces/int16/i16lagm.F
508!||--- calls -----------------------------------------------------
509!|| i16deri ../engine/source/interfaces/int16/i16lagm.F
510!|| i16ni ../engine/source/interfaces/int16/i16lagm.F
511!|| i16rstn ../engine/source/interfaces/int16/i16lagm.F
512!||====================================================================
513 SUBROUTINE i16rst(LLT,R ,S ,T ,NI ,
514 2 NX ,NY ,NZ ,XX ,YY ,ZZ )
515C-----------------------------------------------
516C I m p l i c i t T y p e s
517C-----------------------------------------------
518#include "implicit_f.inc"
519C-----------------------------------------------
520C G l o b a l P a r a m e t e r s
521C-----------------------------------------------
522#include "mvsiz_p.inc"
523C-----------------------------------------------
524C D u m m y A r g u m e n t s
525C-----------------------------------------------
526 INTEGER LLT
527C REAL
528 my_real
529 . XX(MVSIZ,17),YY(MVSIZ,17),ZZ(MVSIZ,17)
530 my_real
531 . R(MVSIZ),S(MVSIZ),T(MVSIZ),NI(MVSIZ,17) ,
532 . NX(MVSIZ),NY(MVSIZ),NZ(MVSIZ)
533C-----------------------------------------------
534C L o c a l V a r i a b l e s
535C-----------------------------------------------
536 INTEGER I,J,IK,NK,ITER,NITERMAX,JTER,NJTERMAX,CONV
537 my_real
538 . vx,vy,vz,vn
539 my_real
540 . drdx(mvsiz),drdy(mvsiz),drdz(mvsiz),
541 . dsdx(mvsiz),dsdy(mvsiz),dsdz(mvsiz),
542 . dtdx(mvsiz),dtdy(mvsiz),dtdz(mvsiz),
543 . dxdr(mvsiz),dydr(mvsiz),dzdr(mvsiz),
544 . dxdt(mvsiz),dydt(mvsiz),dzdt(mvsiz),
545 . rr(mvsiz),ss(mvsiz),tt(mvsiz)
546C-----------------------------------------------
547C
548C r=s=t=0
549C
550C +---> iter
551C |
552C | Ni(r,s,t) =
553C | dNi/dr =
554C | ... _
555C | \
556C | dx/dr = /_ (xi * dNi/dr)
557C | ...
558C |
559C | [dx/dr dy/dr dz/dr]
560C | [J] = |dx/ds dy/ds dz/ds|
561C | [dx/dt dy/dt dz/dt]
562C |
563C | +--> jter
564C | | _
565C | | \
566C | | x(r,s,t) = /_ (xi * Ni(r,s,t))
567C | | ...
568C | |
569C | | |r| |r| -1 |xs-x(r,s,t)|
570C | | {s} = {s} + [J] {ys-y(r,s,t)}
571C | | |t| |t| |zs-z(r,s,t)|
572C | |
573C | | Ni(r,s,t) =
574C +-+---
575C-----------------------------------------------
576 nitermax = 3
577 njtermax = 3
578 conv = 0
579C
580 DO i=1,llt
581 rr(i) = zero
582 ss(i) = zero
583 tt(i) = zero
584 ENDDO
585C-----------------------------------------------
586C calcul de r,s,t et Ni(r,s,t)
587C-----------------------------------------------
588 DO iter=1,nitermax
589c
590c print *, "iter",iter
591c
592C-----------------------------------------------
593C calcul de Ni(r,s,t); [J]; [J]-1
594C-----------------------------------------------
595 CALL i16deri(llt,rr ,ss ,tt ,ni ,
596 2 drdx ,drdy ,drdz ,dsdx ,dsdy ,dsdz ,
597 3 dtdx ,dtdy ,dtdz ,dxdr ,dydr ,dzdr ,
598 4 dxdt ,dydt ,dzdt ,xx ,yy ,zz )
599C
600 DO jter=1,njtermax
601c
602c print *, "jter",jter
603c
604C-----------------------------------------------
605C calcul de r,s,t new
606C-----------------------------------------------
607 CALL i16rstn(llt,rr,ss ,tt ,ni ,conv ,
608 2 drdx ,drdy ,drdz ,dsdx ,dsdy ,dsdz ,
609 3 dtdx ,dtdy ,dtdz ,xx ,yy ,zz ,
610 4 r ,s ,t )
611c
612c print *, "r,s,t",r(1),s(1),t(1)
613c print *, "rr,ss,tt",rr(1),ss(1),tt(1)
614c
615C-----------------------------------------------
616C calcul de Ni(-1<r<1 , -1<s<1 , -1<t<1)
617C-----------------------------------------------
618 CALL i16ni(llt,rr ,ss ,tt ,ni )
619C pb de parith on si conv fonction de mvsiz !!!!!!!
620C IF(CONV/=0)RETURN
621C
622 ENDDO
623 ENDDO
624C
625 DO i=1,llt
626 nx(i) = dydt(i)*dzdr(i) - dzdt(i)*dydr(i)
627 ny(i) = dzdt(i)*dxdr(i) - dxdt(i)*dzdr(i)
628 nz(i) = dxdt(i)*dydr(i) - dydt(i)*dxdr(i)
629 ENDDO
630C
631 RETURN
632 END
633C
634!||====================================================================
635!|| i16ni ../engine/source/interfaces/int16/i16lagm.F
636!||--- called by ------------------------------------------------------
637!|| i16rst ../engine/source/interfaces/int16/i16lagm.F
638!|| i17rst ../engine/source/interfaces/int17/i17lagm.F
639!||====================================================================
640 SUBROUTINE i16ni(LLT,RR ,SS ,TT ,NI )
641C-----------------------------------------------
642C I m p l i c i t T y p e s
643C-----------------------------------------------
644#include "implicit_f.inc"
645C-----------------------------------------------
646C G l o b a l P a r a m e t e r s
647C-----------------------------------------------
648#include "mvsiz_p.inc"
649C-----------------------------------------------
650C D u m m y A r g u m e n t s
651C-----------------------------------------------
652 INTEGER LLT
653 my_real
654 . RR(MVSIZ),SS(MVSIZ),TT(MVSIZ),NI(MVSIZ,17)
655C-----------------------------------------------
656C L o c a l V a r i a b l e s
657C-----------------------------------------------
658 INTEGER I
659 my_real
660 . u_m_r,u_p_r,u_m_s,u_p_s,u_m_t,u_p_t,
661 . ums_umt,ums_upt,ups_umt,ups_upt,
662 . umr_ums,umr_ups,upr_ums,upr_ups,
663 . umt_umr,umt_upr,upt_umr,upt_upr,
664 . a,r05,s05,t05
665C-----------------------------------------------------------------------
666C calcul de Ni
667C-----------------------------------------------------------------------
668 DO i=1,llt
669C
670 r05 = half*rr(i)
671 s05 = half*ss(i)
672 t05 = half*tt(i)
673C
674 u_m_r = half - r05
675 u_p_r = half + r05
676C
677 u_m_s = half - s05
678 u_p_s = half + s05
679C
680 u_m_t = half - t05
681 u_p_t = half + t05
682C
683 ums_umt = u_m_s * u_m_t
684 ums_upt = u_m_s * u_p_t
685 ups_umt = u_p_s * u_m_t
686 ups_upt = u_p_s * u_p_t
687C
688 umr_ums = u_m_r * u_m_s
689 umr_ups = u_m_r * u_p_s
690 upr_ums = u_p_r * u_m_s
691 upr_ups = u_p_r * u_p_s
692C
693 umt_umr = u_m_t * u_m_r
694 umt_upr = u_m_t * u_p_r
695 upt_umr = u_p_t * u_m_r
696 upt_upr = u_p_t * u_p_r
697C
698 a = -rr(i)-tt(i) - one
699 ni(i,1) = u_m_r * ums_umt * a
700 ni(i,5) = u_m_r * ups_umt * a
701 a = -rr(i)+tt(i) - one
702 ni(i,2) = u_m_r * ums_upt * a
703 ni(i,6) = u_m_r * ups_upt * a
704 a = rr(i)+tt(i) - one
705 ni(i,3) = u_p_r * ums_upt * a
706 ni(i,7) = u_p_r * ups_upt * a
707 a = rr(i)-tt(i) - one
708 ni(i,4) = u_p_r * ums_umt * a
709 ni(i,8) = u_p_r * ups_umt * a
710C------------------------------------
711 a = (one - rr(i)*rr(i))
712 ni(i,10) = a * ums_upt
713 ni(i,12) = a * ums_umt
714 ni(i,14) = a * ups_upt
715 ni(i,16) = a * ups_umt
716C------------------------------------
717 a = (one - tt(i)*tt(i))
718 ni(i,9) = a * umr_ums
719 ni(i,11) = a * upr_ums
720 ni(i,13) = a * umr_ups
721 ni(i,15) = a * upr_ups
722C------------------------------------
723 ni(i,17) = -one
724C------------------------------------
725 ENDDO
726C-----------------------------------------------
727 RETURN
728 END
729!||====================================================================
730!|| i16rstn ../engine/source/interfaces/int16/i16lagm.F
731!||--- called by ------------------------------------------------------
732!|| i16rst ../engine/source/interfaces/int16/i16lagm.F
733!|| i17rst ../engine/source/interfaces/int17/i17lagm.F
734!||====================================================================
735 SUBROUTINE i16rstn(LLT,RR,SS ,TT ,NI ,CONV ,
736 2 DRDX ,DRDY ,DRDZ ,DSDX ,DSDY ,DSDZ ,
737 3 DTDX ,DTDY ,DTDZ ,XX ,YY ,ZZ ,
738 4 R ,S ,T )
739C-----------------------------------------------
740C I m p l i c i t T y p e s
741C-----------------------------------------------
742c#include "implicit_f.inc"
743 implicit none
744C-----------------------------------------------
745C G l o b a l P a r a m e t e r s
746C-----------------------------------------------
747#include "mvsiz_p.inc"
748#include "constant.inc"
749C-----------------------------------------------
750C D u m m y A r g u m e n t s
751C-----------------------------------------------
752 INTEGER LLT,CONV
753 my_real
754 . r(mvsiz),s(mvsiz),t(mvsiz),ni(mvsiz,17) ,
755 . rr(mvsiz),ss(mvsiz),tt(mvsiz),
756 . xx(mvsiz,17) ,yy(mvsiz,17) ,zz(mvsiz,17) ,
757 . drdx(mvsiz),drdy(mvsiz),drdz(mvsiz),
758 . dsdx(mvsiz),dsdy(mvsiz),dsdz(mvsiz),
759 . dtdx(mvsiz),dtdy(mvsiz),dtdz(mvsiz)
760C-----------------------------------------------
761C L o c a l V a r i a b l e s
762C-----------------------------------------------
763 INTEGER I
764 my_real
765 . dx ,dy,dz,dr ,ds,dt,err
766C
767 err=zero
768C-----------------------------------------------
769 DO i=1,llt
770C
771 dx = xx(i,17)
772 + - ni(i, 1)*xx(i, 1) - ni(i, 2)*xx(i, 2) - ni(i, 3)*xx(i, 3)
773 + - ni(i, 4)*xx(i, 4) - ni(i, 5)*xx(i, 5) - ni(i, 6)*xx(i, 6)
774 + - ni(i, 7)*xx(i, 7) - ni(i, 8)*xx(i, 8) - ni(i, 9)*xx(i, 9)
775 + - ni(i,10)*xx(i,10) - ni(i,11)*xx(i,11) - ni(i,12)*xx(i,12)
776 + - ni(i,13)*xx(i,13) - ni(i,14)*xx(i,14) - ni(i,15)*xx(i,15)
777 + - ni(i,16)*xx(i,16)
778 dy = yy(i,17)
779 + - ni(i, 1)*yy(i, 1) - ni(i, 2)*yy(i, 2) - ni(i, 3)*yy(i, 3)
780 + - ni(i, 4)*yy(i, 4) - ni(i, 5)*yy(i, 5) - ni(i, 6)*yy(i, 6)
781 + - ni(i, 7)*yy(i, 7) - ni(i, 8)*yy(i, 8) - ni(i, 9)*yy(i, 9)
782 + - ni(i,10)*yy(i,10) - ni(i,11)*yy(i,11) - ni(i,12)*yy(i,12)
783 + - ni(i,13)*yy(i,13) - ni(i,14)*yy(i,14) - ni(i,15)*yy(i,15)
784 + - ni(i,16)*yy(i,16)
785 dz = zz(i,17)
786 + - ni(i, 1)*zz(i, 1) - ni(i, 2)*zz(i, 2) - ni(i, 3)*zz(i, 3)
787 + - ni(i, 4)*zz(i, 4) - ni(i, 5)*zz(i, 5) - ni(i, 6)*zz(i, 6)
788 + - ni(i, 7)*zz(i, 7) - ni(i, 8)*zz(i, 8) - ni(i, 9)*zz(i, 9)
789 + - ni(i,10)*zz(i,10) - ni(i,11)*zz(i,11) - ni(i,12)*zz(i,12)
790 + - ni(i,13)*zz(i,13) - ni(i,14)*zz(i,14) - ni(i,15)*zz(i,15)
791 + - ni(i,16)*zz(i,16)
792C
793 dr = drdx(i)*dx + drdy(i)*dy + drdz(i)*dz
794 ds = dsdx(i)*dx + dsdy(i)*dy + dsdz(i)*dz
795 dt = dtdx(i)*dx + dtdy(i)*dy + dtdz(i)*dz
796C
797c
798c print *, "DRDX(I),DRDY(I),DRDZ(I)",DRDX(I),DRDY(I),DRDZ(I)
799c print *, "DSDX(I),DSDY(I),DSDZ(I)",DSDX(I),DSDY(I),DSDZ(I)
800c print *, "DTDX(I),DTDY(I),DTDZ(I)",DTDX(I),DTDY(I),DTDZ(I)
801c print *, "Ni",ni(1,1),ni(1,2),ni(1,3),ni(1,4),ni(1,5),ni(1,9)
802c print *, "dx,dy,dz",dx ,dy ,dz
803c
804 rr(i) = rr(i) + dr
805 ss(i) = ss(i) + ds
806 tt(i) = tt(i) + dt
807C
808 r(i) = rr(i)
809 s(i) = ss(i)
810 t(i) = tt(i)
811C
812 IF(r(i)>=-one.AND.s(i)>=-one.AND.t(i)>=-one.AND.
813 . r(i)<= one.AND.s(i)<= one.AND.t(i)<= one)THEN
814 err = max(err,abs(dr),abs(ds),abs(dt))
815 ELSE
816 rr(i) = max(min(rr(i),one),-one)
817 ss(i) = max(min(ss(i),one),-one)
818 tt(i) = max(min(tt(i),one),-one)
819 ENDIF
820c
821c print *, "3r,s,t",r(1),s(1),t(1)
822c print *, "3rr,ss,tt",rr(1),ss(1),tt(1)
823c print *, "dr,ds,dt",dr ,ds ,dt
824c print *, "r,s,t",r(1),s(1),t(1)
825c print *, "ERR",ERR
826c
827C
828 ENDDO
829C
830 IF(err<em4) conv = 1
831C-----------------------------------------------
832 RETURN
833 END
834!||====================================================================
835!|| i16deri ../engine/source/interfaces/int16/i16lagm.f
836!||--- called by ------------------------------------------------------
837!|| i16rst ../engine/source/interfaces/int16/i16lagm.F
838!|| i17rst ../engine/source/interfaces/int17/i17lagm.F
839!||====================================================================
840 SUBROUTINE i16deri(LLT,RR,SS ,TT ,NI ,
841 2 DRDX ,DRDY ,DRDZ ,DSDX ,DSDY ,DSDZ ,
842 3 DTDX ,DTDY ,DTDZ ,DXDR ,DYDR ,DZDR ,
843 4 DXDT ,DYDT ,DZDT ,XX ,YY ,ZZ )
844C-----------------------------------------------
845C I m p l i c i t T y p e s
846C-----------------------------------------------
847#include "implicit_f.inc"
848C-----------------------------------------------
849C G l o b a l P a r a m e t e r s
850C-----------------------------------------------
851#include "mvsiz_p.inc"
852C-----------------------------------------------
853C D u m m y A r g u m e n t s
854C-----------------------------------------------
855 INTEGER LLT
856C REAL
857 my_real
858 . DXDR(MVSIZ), DYDR(MVSIZ), DZDR(MVSIZ),
859 . DXDT(MVSIZ), DYDT(MVSIZ), DZDT(MVSIZ),
860 . DRDX(MVSIZ), DSDX(MVSIZ), DTDX(MVSIZ),
861 . DRDY(MVSIZ), DSDY(MVSIZ), DTDY(MVSIZ),
862 . DRDZ(MVSIZ), DSDZ(MVSIZ), DTDZ(MVSIZ),
863 . XX(MVSIZ,17) ,YY(MVSIZ,17),ZZ(MVSIZ,17),
864 . NI(MVSIZ,17) ,RR(MVSIZ) ,SS(MVSIZ) ,TT(MVSIZ)
865C-----------------------------------------------
866C L o c a l V a r i a b l e s
867C-----------------------------------------------
868 INTEGER I,N
869 my_real
870 . dxds(mvsiz), dyds(mvsiz), dzds(mvsiz),
871 . dnidr(16),dnids(16),dnidt(16),
872 . d, aa, bb, det(mvsiz),r9 ,r13 ,s9 ,s10 ,s11 ,s12 ,t10 ,t14
873 my_real
874 . u_m_r,u_p_r,u_m_s,u_p_s,u_m_t,u_p_t,
875 . ums_umt,ums_upt,ups_umt,ups_upt,
876 . umr_ums,umr_ups,upr_ums,upr_ups,
877 . umt_umr,umt_upr,upt_umr,upt_upr,
878 . a,r05,s05,t05
879C-----------------------------------------------
880C/*
881C
882C
883C
884C ^ S _ T
885C | /|
886C | /
887C | /
888C ( 7)===========|===(14)==============( 6)
889C //| | / //|
890C // | | / //||
891C // | | / // ||
892C // | (Is) | + // ||
893C // | | / // ||
894C (15) | + / (13) ||
895C // | | / // ||
896C // ( 3)-----------|---(10)-------//-----( 2)
897C // / |/ // //
898C R <-------//- -+ - - - - - - - - -# - - - - - -//- -+ //
899C // / / // //
900C ( 8)===============(16)==============( 5) //
901C || / / || //
902C || (11) / || ( 9)
903C || / / || //
904C || / + || //
905C || / || //
906C || / ||//
907C ||/ ||/
908C ( 4)===============(12)==============( 1)
909C
910C
911C*/
912C-----------------------------------------------
913C
914C-----------------------------------------------
915C _
916C \
917C x(r,s,t) = /_ (xi * Ni(r,s,t))
918C _
919C \
920C y(r,s,t) = /_ (yi * Ni(r,s,t))
921C _
922C \
923C z(r,s,t) = /_ zi * Ni(r,s,t))
924C
925C _
926C \
927C dx/dr = /_ (xi * dNi/dr)
928C ...
929C
930C [dx/dr dy/dr dz/dr]
931C [J] = |dx/ds dy/ds dz/ds|
932C [dx/dt dy/dt dz/dt]
933C
934C |r| |r| -1 |xs-x|
935C {s} = {s} + [J] {ys-y}
936C |t| |t| |zs-z|
937C
938C-----------------------------------------------------------------------
939C Ni; dNi/dr; dNi/ds; dNi/dt
940C-----------------------------------------------------------------------
941 DO i=1,llt
942 r05 = half*rr(i)
943 s05 = half*ss(i)
944 t05 = half*tt(i)
945C
946 u_m_r = half - r05
947 u_p_r = half + r05
948C
949 u_m_s = half - s05
950 u_p_s = half + s05
951C
952 u_m_t = half - t05
953 u_p_t = half + t05
954C
955 ums_umt = u_m_s * u_m_t
956 ums_upt = u_m_s * u_p_t
957 ups_umt = u_p_s * u_m_t
958 ups_upt = u_p_s * u_p_t
959C
960 umr_ums = u_m_r * u_m_s
961 umr_ups = u_m_r * u_p_s
962 upr_ums = u_p_r * u_m_s
963 upr_ups = u_p_r * u_p_s
964C
965 umt_umr = u_m_t * u_m_r
966 umt_upr = u_m_t * u_p_r
967 upt_umr = u_p_t * u_m_r
968 upt_upr = u_p_t * u_p_r
969C
970 a = -rr(i)-tt(i)-one
971 ni(i,1) = u_m_r * ums_umt * a
972 ni(i,5) = u_m_r * ups_umt * a
973 a = -rr(i)+tt(i)-one
974 ni(i,2) = u_m_r * ums_upt * a
975 ni(i,6) = u_m_r * ups_upt * a
976 a = rr(i)+tt(i)-one
977 ni(i,3) = u_p_r * ums_upt * a
978 ni(i,7) = u_p_r * ups_upt * a
979 a = rr(i)-tt(i)-one
980 ni(i,4) = u_p_r * ums_umt * a
981 ni(i,8) = u_p_r * ups_umt * a
982C
983 a = -t05 - rr(i)
984 dnidr(1) = -ums_umt * a
985 dnidr(5) = -ups_umt * a
986 a = t05 - rr(i)
987 dnidr(2) = -ums_upt * a
988 dnidr(6) = -ups_upt * a
989 a = t05 + rr(i)
990 dnidr(3) = ums_upt * a
991 dnidr(7) = ups_upt * a
992 a = -t05 + rr(i)
993 dnidr(4) = ums_umt * a
994 dnidr(8) = ups_umt * a
995C
996 a = -r05 - u_p_t
997 dnids(1) = -umt_umr * a
998 dnids(5) = umt_umr * a
999 a = -r05 - u_m_t
1000 dnids(2) = -upt_umr * a
1001 dnids(6) = upt_umr * a
1002 a = r05 - u_m_t
1003 dnids(3) = -upt_upr * a
1004 dnids(7) = upt_upr * a
1005 a = r05 - u_p_t
1006 dnids(4) = -umt_upr * a
1007 dnids(8) = umt_upr * a
1008C
1009 a = -r05 - tt(i)
1010 dnidt(1) = -umr_ums * a
1011 dnidt(5) = -umr_ups * a
1012 a = -r05 + tt(i)
1013 dnidt(2) = umr_ums * a
1014 dnidt(6) = umr_ups * a
1015 a = +r05 + tt(i)
1016 dnidt(3) = upr_ums * a
1017 dnidt(7) = upr_ups * a
1018 a = +r05 - tt(i)
1019 dnidt(4) = -upr_ums * a
1020 dnidt(8) = -upr_ups * a
1021C------------------------------------
1022 a = (one-rr(i)*rr(i))
1023 ni(i,10) = a * ums_upt
1024 ni(i,12) = a * ums_umt
1025 ni(i,14) = a * ups_upt
1026 ni(i,16) = a * ups_umt
1027C
1028 a = half*a
1029 dnidt(10) = a * u_m_s
1030 dnidt(14) = a * u_p_s
1031 dnids(10) = -a * u_p_t
1032 dnids(12) = -a * u_m_t
1033C
1034 a = -two*rr(i)
1035 dnidr(10) = a * ums_upt
1036 dnidr(12) = a * ums_umt
1037 dnidr(14) = a * ups_upt
1038 dnidr(16) = a * ups_umt
1039C------------------------------------
1040 a = (one-tt(i)*tt(i))
1041 ni(i,9) = a * umr_ums
1042 ni(i,11) = a * upr_ums
1043 ni(i,13) = a * umr_ups
1044 ni(i,15) = a * upr_ups
1045C
1046 ni(i,17) = -one
1047C
1048 a = half*a
1049 dnidr(9) = -a * u_m_s
1050 dnidr(13) = -a * u_p_s
1051C
1052 dnids(9) = -a * u_m_r
1053 dnids(11) = -a * u_p_r
1054C
1055 a = -two*tt(i)
1056 dnidt(9) = a * umr_ums
1057 dnidt(11) = a * upr_ums
1058 dnidt(13) = a * umr_ups
1059 dnidt(15) = a * upr_ups
1060c
1061c print *, "DNIDr(1),DNIDr(9)",DNIDr(1),DNIDr(9)
1062c print *, "DNIDs(1),DNIDs(9)",DNIDs(1),DNIDs(9)
1063c print *, "DNIDT(1),DNIDT(9)",DNIDT(1),DNIDT(9)
1064c print *, "XX(I,1),XX(I,9)",XX(I,1),XX(I,9)
1065c
1066C-----------------------------------------------------------------------
1067C dx/dr; dx/ds; dx/dt
1068C-----------------------------------------------------------------------
1069 dxdr(i) = dnidr(1)*xx(i,1) + dnidr(2)*xx(i,2) + dnidr(3)*xx(i,3)
1070 + + dnidr(4)*xx(i,4) + dnidr(5)*xx(i,5) + dnidr(6)*xx(i,6)
1071 + + dnidr(7)*xx(i,7) + dnidr(8)*xx(i,8)
1072 + + dnidr(9)*(xx(i,9) - xx(i,11)) + dnidr(10)*xx(i,10)
1073 + + dnidr(12)*xx(i,12) + dnidr(13)*(xx(i,13) - xx(i,15))
1074 + + dnidr(14)*xx(i,14) + dnidr(16)*xx(i,16)
1075C
1076 dxds(i) = dnids(1)*xx(i,1) + dnids(2)*xx(i,2) + dnids(3)*xx(i,3)
1077 + + dnids(4)*xx(i,4) + dnids(5)*xx(i,5) + dnids(6)*xx(i,6)
1078 + + dnids(7)*xx(i,7) + dnids(8)*xx(i,8)
1079 + + dnids(9)* (xx(i,9) - xx(i,13))
1080 + + dnids(10)*(xx(i,10) - xx(i,14))
1081 + + dnids(11)*(xx(i,11) - xx(i,15))
1082 + + dnids(12)*(xx(i,12) - xx(i,16))
1083C
1084 dxdt(i) = dnidt(1)*xx(i,1) + dnidt(2)*xx(i,2) + dnidt(3)*xx(i,3)
1085 + + dnidt(4)*xx(i,4) + dnidt(5)*xx(i,5) + dnidt(6)*xx(i,6)
1086 + + dnidt(7)*xx(i,7) + dnidt(8)*xx(i,8)
1087 + + dnidt(9)*xx(i,9) + dnidt(10)*(xx(i,10) - xx(i,12))
1088 + + dnidt(11)*xx(i,11) + dnidt(13)*xx(i,13)
1089 + + dnidt(14)*(xx(i,14) - xx(i,16)) + dnidt(15)*xx(i,15)
1090C-----------------------------------------------------------------------
1091C dy/dr; dy/ds; dy/dt
1092C-----------------------------------------------------------------------
1093 dydr(i) = dnidr(1)*yy(i,1) + dnidr(2)*yy(i,2) + dnidr(3)*yy(i,3)
1094 + + dnidr(4)*yy(i,4) + dnidr(5)*yy(i,5) + dnidr(6)*yy(i,6)
1095 + + dnidr(7)*yy(i,7) + dnidr(8)*yy(i,8)
1096 + + dnidr(9)*(yy(i,9) - yy(i,11)) + dnidr(10)*yy(i,10)
1097 + + dnidr(12)*yy(i,12) + dnidr(13)*(yy(i,13) - yy(i,15))
1098 + + dnidr(14)*yy(i,14) + dnidr(16)*yy(i,16)
1099C
1100 dyds(i) = dnids(1)*yy(i,1) + dnids(2)*yy(i,2) + dnids(3)*yy(i,3)
1101 + + dnids(4)*yy(i,4) + dnids(5)*yy(i,5) + dnids(6)*yy(i,6)
1102 + + dnids(7)*yy(i,7) + dnids(8)*yy(i,8)
1103 + + dnids(9)* (yy(i,9) - yy(i,13))
1104 + + dnids(10)*(yy(i,10) - yy(i,14))
1105 + + dnids(11)*(yy(i,11) - yy(i,15))
1106 + + dnids(12)*(yy(i,12) - yy(i,16))
1107C
1108 dydt(i) = dnidt(1)*yy(i,1) + dnidt(2)*yy(i,2) + dnidt(3)*yy(i,3)
1109 + + dnidt(4)*yy(i,4) + dnidt(5)*yy(i,5) + dnidt(6)*yy(i,6)
1110 + + dnidt(7)*yy(i,7) + dnidt(8)*yy(i,8)
1111 + + dnidt(9)*yy(i,9) + dnidt(10)*(yy(i,10) - yy(i,12))
1112 + + dnidt(11)*yy(i,11) + dnidt(13)*yy(i,13)
1113 + + dnidt(14)*(yy(i,14) - yy(i,16)) + dnidt(15)*yy(i,15)
1114C-----------------------------------------------------------------------
1115C dz/dr; dz/ds; dz/dt
1116C-----------------------------------------------------------------------
1117 dzdr(i) = dnidr(1)*zz(i,1) + dnidr(2)*zz(i,2) + dnidr(3)*zz(i,3)
1118 + + dnidr(4)*zz(i,4) + dnidr(5)*zz(i,5) + dnidr(6)*zz(i,6)
1119 + + dnidr(7)*zz(i,7) + dnidr(8)*zz(i,8)
1120 + + dnidr(9)*(zz(i,9) - zz(i,11)) + dnidr(10)*zz(i,10)
1121 + + dnidr(12)*zz(i,12) + dnidr(13)*(zz(i,13) - zz(i,15))
1122 + + dnidr(14)*zz(i,14) + dnidr(16)*zz(i,16)
1123C
1124 dzds(i) = dnids(1)*zz(i,1) + dnids(2)*zz(i,2) + dnids(3)*zz(i,3)
1125 + + dnids(4)*zz(i,4) + dnids(5)*zz(i,5) + dnids(6)*zz(i,6)
1126 + + dnids(7)*zz(i,7) + dnids(8)*zz(i,8)
1127 + + dnids(9)* (zz(i,9) - zz(i,13))
1128 + + dnids(10)*(zz(i,10) - zz(i,14))
1129 + + dnids(11)*(zz(i,11) - zz(i,15))
1130 + + dnids(12)*(zz(i,12) - zz(i,16))
1131C
1132 dzdt(i) = dnidt(1)*zz(i,1) + dnidt(2)*zz(i,2) + dnidt(3)*zz(i,3)
1133 + + dnidt(4)*zz(i,4) + dnidt(5)*zz(i,5) + dnidt(6)*zz(i,6)
1134 + + dnidt(7)*zz(i,7) + dnidt(8)*zz(i,8)
1135 + + dnidt(9)*zz(i,9) + dnidt(10)*(zz(i,10) - zz(i,12))
1136 + + dnidt(11)*zz(i,11) + dnidt(13)*zz(i,13)
1137 + + dnidt(14)*(zz(i,14) - zz(i,16)) + dnidt(15)*zz(i,15)
1138C-----------------------------------------------------------------------
1139C -1
1140C [J] Inversion du jacobien
1141C-----------------------------------------------------------------------
1142 drdx(i)=dyds(i)*dzdt(i)-dzds(i)*dydt(i)
1143 drdy(i)=dzds(i)*dxdt(i)-dxds(i)*dzdt(i)
1144 drdz(i)=dxds(i)*dydt(i)-dyds(i)*dxdt(i)
1145C
1146 dsdz(i)=dxdt(i)*dydr(i)-dydt(i)*dxdr(i)
1147 dsdy(i)=dzdt(i)*dxdr(i)-dxdt(i)*dzdr(i)
1148 dsdx(i)=dydt(i)*dzdr(i)-dzdt(i)*dydr(i)
1149C
1150 dtdx(i)=dydr(i)*dzds(i)-dzdr(i)*dyds(i)
1151 dtdy(i)=dzdr(i)*dxds(i)-dxdr(i)*dzds(i)
1152 dtdz(i)=dxdr(i)*dyds(i)-dydr(i)*dxds(i)
1153C
1154 det(i) = dxdr(i) * drdx(i)
1155 . + dydr(i) * drdy(i)
1156 . + dzdr(i) * drdz(i)
1157C
1158c
1159c print *, "Det",DET(I)
1160c print *, "DXDR(I),DYDR(I),DZDR(I)",DXDR(I),DYDR(I),DZDR(I)
1161c print *, "DXDs(I),DYDs(I),DZDs(I)",DXDs(I),DYDs(I),DZDs(I)
1162c print *, "DXDt(I),DYDt(I),DZDt(I)",DXDt(I),DYDt(I),DZDt(I)
1163c
1164 ENDDO
1165C
1166 DO i=1,llt
1167C-----------------------------------------------------------------------
1168C -1
1169C [J] Inversion du jacobien suite
1170C-----------------------------------------------------------------------
1171 d = one/det(i)
1172 drdx(i)=d*drdx(i)
1173 dsdx(i)=d*dsdx(i)
1174 dtdx(i)=d*dtdx(i)
1175C
1176 drdy(i)=d*drdy(i)
1177 dsdy(i)=d*dsdy(i)
1178 dtdy(i)=d*dtdy(i)
1179C
1180 drdz(i)=d*drdz(i)
1181 dsdz(i)=d*dsdz(i)
1182 dtdz(i)=d*dtdz(i)
1183C
1184c
1185c print *, "DRDX(I),DRDY(I),DRDZ(I)",DRDX(I),DRDY(I),DRDZ(I)
1186c print *, "DSDX(I),DSDY(I),DSDZ(I)",DSDX(I),DSDY(I),DSDZ(I)
1187c print *, "DTDX(I),DTDY(I),DTDZ(I)",DTDX(I),DTDY(I),DTDZ(I)
1188c
1189 ENDDO
1190C-----------------------------------------------------------------------
1191 RETURN
1192 END
my_real
#define my_real
Definition cppsort.cpp:32
i16rst
subroutine i16rst(llt, r, s, t, ni, nx, ny, nz, xx, yy, zz)
Definition i16lagm.F:515
i16lagm
subroutine i16lagm(x, v, lll, jll, sll, xll, candn, cande, i_stok, ixs, ixs16, iadll, eminx, nsv, nelem, n_mul_mx, itask, a, itied, nint, nkmax, comntag)
Definition i16lagm.F:35
i16lll
subroutine i16lll(llt, lll, jll, sll, xll, v, xx, yy, zz, iii, iadll, n_mul_mx, a, x, itied, nint, nkmax, comntag)
Definition i16lagm.F:219
i16ni
subroutine i16ni(llt, rr, ss, tt, ni)
Definition i16lagm.F:641
i16rstn
subroutine i16rstn(llt, rr, ss, tt, ni, conv, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, xx, yy, zz, r, s, t)
Definition i16lagm.F:739
i16deri
subroutine i16deri(llt, rr, ss, tt, ni, drdx, drdy, drdz, dsdx, dsdy, dsdz, dtdx, dtdy, dtdz, dxdr, dydr, dzdr, dxdt, dydt, dzdt, xx, yy, zz)
Definition i16lagm.F:844
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