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