OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
i16crit.F File Reference
#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "com04_c.inc"
#include "com08_c.inc"
#include "task_c.inc"
#include "lockon.inc"
#include "lockoff.inc"

Go to the source code of this file.

Functions/Subroutines

subroutine i16crit (x, nsv, nelem, nsn, eminx, nme, itask, xsav, ixs, ixs16, ixs20, ixs10, v, a, xmsrg, xslvg)
subroutine i16box (lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs16, size, xmsr, index, xsav)
subroutine i20box (lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs20, size, xmsr, index, xsav)
subroutine i10box (lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs10, size, xmsr, index, xsav)
subroutine i8box (lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, size, xmsr, index, xsav)

Function/Subroutine Documentation

◆ i10box()

subroutine i10box ( integer lft,
integer llt,
integer, dimension(*) nelem,
eminx,
integer nmef,
integer nmel,
x,
v,
a,
integer, dimension(nixs,*) ixs,
integer, dimension(6,*) ixs10,
size,
xmsr,
integer, dimension(*) index,
xsav )

Definition at line 865 of file i16crit.F.

868 use element_mod , only : nixs
869C-----------------------------------------------
870C I m p l i c i t T y p e s
871C-----------------------------------------------
872#include "implicit_f.inc"
873C-----------------------------------------------
874C C o m m o n B l o c k s
875C-----------------------------------------------
876#include "com04_c.inc"
877#include "com08_c.inc"
878C-----------------------------------------------
879C D u m m y A r g u m e n t s
880C-----------------------------------------------
881 INTEGER LFT ,LLT,NMEF,NMEL,
882 . IXS(NIXS,*),IXS10(6,*),NELEM(*),INDEX(*)
883C REAL
884 my_real
885 . x(3,*),v(3,*),a(3,*),eminx(6,*),SIZE,xmsr(*),xsav(3,*)
886C-----------------------------------------------
887C L o c a l V a r i a b l e s
888C-----------------------------------------------
889 INTEGER I,J,L,NE,IDIR,N10
890 my_real
891 .
892 . x1,x2,x3,x4,x5,x6,x7,x8,
893 . x9,x10,xc,xx,xn
894C------------------------------------
895C calculation of element bounds
896C------------------------------------
897C-----------------------------------------------------------------------
898C Face 1 2 3 4 or 5 6 7 8
899C-----------------------------------------------------------------------
900 DO idir=1,3
901C-----------------------------------------------------------------------
902C x y or z
903C-----------------------------------------------------------------------
904 DO l=lft,llt
905 i = index(l)
906 ne = nelem(i)
907 n10= ne - numels8
908C-----------------------------------------------------------------------
909 j = ixs(2,ne)
910 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
911 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
912 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
913 j = ixs(4,ne)
914 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
915 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
916 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
917 j = ixs(6,ne)
918 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
919 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
920 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
921 j = ixs(7,ne)
922 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
923 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
924 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
925C
926 j = ixs10(1,n10)
927 x5 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
928 xmsr(idir) =max(xmsr(idir) ,x5-xsav(idir,j))
929 xmsr(idir+3)=min(xmsr(idir+3),x5-xsav(idir,j))
930 j = ixs10(2,n10)
931 x6 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
932 xmsr(idir) =max(xmsr(idir) ,x6-xsav(idir,j))
933 xmsr(idir+3)=min(xmsr(idir+3),x6-xsav(idir,j))
934 j = ixs10(3,n10)
935 x7 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
936 xmsr(idir) =max(xmsr(idir) ,x7-xsav(idir,j))
937 xmsr(idir+3)=min(xmsr(idir+3),x7-xsav(idir,j))
938 j = ixs10(4,n10)
939 x8 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
940 xmsr(idir) =max(xmsr(idir) ,x8-xsav(idir,j))
941 xmsr(idir+3)=min(xmsr(idir+3),x8-xsav(idir,j))
942 j = ixs10(5,n10)
943 x9 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
944 xmsr(idir) =max(xmsr(idir) ,x9-xsav(idir,j))
945 xmsr(idir+3)=min(xmsr(idir+3),x9-xsav(idir,j))
946 j = ixs10(6,n10)
947 x10 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
948 xmsr(idir) =max(xmsr(idir) ,x10-xsav(idir,j))
949 xmsr(idir+3)=min(xmsr(idir+3),x10-xsav(idir,j))
950C-----------------------------------------------------------------------
951 xx=max(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,x9,x10)
952 xn=min(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 ,x9,x10)
953 eminx(idir,i) = min( eminx(idir,i) , xn )
954 eminx(idir+3,i) = max( eminx(idir+3,i), xx )
955C-----------------------------------------------------------------------
956C Face 1 2 3 4
957C-----------------------------------------------------------------------
958 xc = (two*(x5+x6+x7) - (x1+x2+x3))* third
959 eminx(idir,i) = min( eminx(idir,i) , xc )
960 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
961C
962 xc = (two*(x5+x8+x9) - (x1+x2+x4))*third
963 eminx(idir,i) = min( eminx(idir,i) , xc )
964 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
965C
966 xc = (two*(x6+x9+x10) - (x2+x3+x4)) * third
967 eminx(idir,i) = min( eminx(idir,i) , xc )
968 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
969C
970 xc = (two*(x7+x8+x10) - (x3+x1+x4)) * third
971 eminx(idir,i) = min( eminx(idir,i) , xc )
972 eminx(idir+3,i) = max( eminx(idir+3,i), xc )
973C
974C-----------------------------------------------------------------------
975C
976 ENDDO
977 ENDDO
978C--------------------------------------------------------------
979C
980 RETURN
#define my_real
Definition cppsort.cpp:32
#define min(a, b)
Definition macros.h:20
#define max(a, b)
Definition macros.h:21

◆ i16box()

subroutine i16box ( integer lft,
integer llt,
integer, dimension(*) nelem,
eminx,
integer nmef,
integer nmel,
x,
v,
a,
integer, dimension(nixs,*) ixs,
integer, dimension(8,*) ixs16,
size,
xmsr,
integer, dimension(*) index,
xsav )

Definition at line 227 of file i16crit.F.

230 use element_mod , only : nixs
231C-----------------------------------------------
232C I m p l i c i t T y p e s
233C-----------------------------------------------
234#include "implicit_f.inc"
235C-----------------------------------------------
236C C o m m o n B l o c k s
237C-----------------------------------------------
238#include "com04_c.inc"
239#include "com08_c.inc"
240C-----------------------------------------------
241C D u m m y A r g u m e n t s
242C-----------------------------------------------
243 INTEGER LFT ,LLT,NMEF,NMEL,
244 . IXS(NIXS,*),IXS16(8,*),NELEM(*),INDEX(*)
245C REAL
246 my_real
247 . x(3,*),v(3,*),a(3,*),eminx(6,*),SIZE,xmsr(*),xsav(3,*)
248C-----------------------------------------------
249C L o c a l V a r i a b l e s
250C-----------------------------------------------
251 INTEGER I,J,L,NE,IFACE,IDIR,IPERM(8,2),N16
252 my_real
253 . an,ax,bn,bx,cn,cx,dx,dn,d4,d8,x1,x2,x3,x4,
254 . x9,x10,x11,x12,xc
255 DATA iperm / 2, 3, 4, 5, 1, 2, 3, 4,
256 2 6, 7, 8, 9, 5, 6, 7, 8/
257C
258C-----------------------------------------------
259C/*
260C
261C ( 7)8==============(14)6=============( 6)7
262C //| //|
263C // | //||
264C // | // ||
265C // | // ||
266C // | // ||
267C (15)7 | (13)5 ||
268C // | // ||
269C // ( 3)4--------------(10)2------//-----( 2)3
270C // / // //
271C // / // //
272C // / // //
273C ( 8)9==============(16)8=============( 5)6 //
274C || / || //
275C || (11)3 (C) || ( 9)1
276C || / || //
277C || / || //
278C || / || //
279C || / ||//
280C ||/ ||/
281C ( 4)5==============(12)4=============( 1)2
282C
283C*/
284C---------------------------------------------------------
285C MONODIM
286C---------------------------------------------------------
287C
288C x
289C / \
290C / (2)
291C (3)
292C /
293C /
294C (1)
295C
296C N1 = -0.5 (1-r)r dN1/dr = r - 0.5
297C N2 = 0.5 (1+r)r dN2/dr = r + 0.5
298C N3 = (1-r^2) dN3/dr = - 2 r
299C
300C x = N1 x1 + N2 x2 + N3 x3
301C x = -0.5 (1-r)r x1 + 0.5 (1+r)r x2 + (1-r^2) x3
302C 2 x = (x1 + x2 - 2 x3) r^2 + (x2-x1) r + 2 x3
303C
304C 0) Point XMAX
305C
306C dx/dr = (r - 0.5) x1 + (r + 0.5) x2 - 2 r x3 = 0
307C r = 0.5 (x1 - x2) / (x1 + x2 - 2 x3)
308C
309C 2 x (x1 + x2 - 2 x3) = (x1 + x2 - 2 x3)^2 r^2
310C + (x2-x1)(x1 + x2 - 2 x3) r
311C + 2 (x1 + x2 - 2 x3)x3
312C
313C 2 x (x1 + x2 - 2 x3) = - 0.25 (x1 - x2)^2
314C + 2 (x1 + x2 - 2 x3)x3
315C
316C x = x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
317C
318C------------------------------------------------------------
319C solution 0 => x < x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
320C si x3 -> (x1 + x2)/ 2 x -> infini
321C------------------------------------------------------------
322C
323C 1) search for the xmax point between 1 and 2
324C
325C si r > 1
326C => x = x2
327C
328C si r < -1
329C => x = x1
330C
331C si -1 < r < 1
332C => -1/4 < 0.125 (x1 - x2) / (x1 + x2 - 2 x3) < 1/4
333C
334C (x - x3)/(x1 - x2) = - 0.125 (x1 - x2) / (x1 + x2 - 2 x3)
335C -1/4 < (x - x3)/ < 1/4
336C
337C si x2 > x1 x3 -1/4 (x2 - x1) < x < x3 + 1/4 (x2 - x1)
338C si x2 < x1 x3 -1/4 (x2 - x1) > x > x3 + 1/4 (x2 - x1)
339C
340C => x3 - 1/4 |x2 - x1| < x < x3 + 1/4 |x2 - x1|
341C
342C------------------------------------------------------------
343C solution 1 => x < max (x1 , x2 , x3 + 1/4 |x2 - x1|)
344C------------------------------------------------------------
345C
346C 2) search for the most unfavorable position of x3
347C
348C x = x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
349C dx/dx3 = 1 - 0.25 (x1 - x2)^2 / (x1 + x2 - 2 x3)^2 = 0
350C (2 x3 - x1 - x2 )^2 = 0.25 (x2 - x1)^2
351C if x2 > x1 and x3 > (x1 + x2)/2 or
352C si x2 < x1 et x3 < (x1 + x2)/2
353C (2 x3 - x1 - x2 ) = 0.5 (x2 - x1)
354C x3 = (x1 + x2)/2 + (x2 - x1)/4
355C if x2 > x1 and x3 < (x1 + x2)/2 or
356C si x2 < x1 et x3 > (x1 + x2)/2
357C no solution
358C
359C si x3 < (x1 + x2)/2 + (x2 - x1)/4
360C => x < max (x1 , x2 ) et
361C => x> min (x1, x2) => Verified by solution 1 and 2
362C
363C si x3 = (x1 + x2)/2 + (x2 - x1)/4
364C x = x3 + 0.125 (x1 - x2)^2 / (2 x3 - x1 - x2 )
365C x = (x1 + x2)/2 + (x2 - x1)/4
366C + 0.125 (x1 - x2)^2 / (2 ((x1 + x2)/2 + (x2 - x1)/4) - x1 - x2)
367C x = (x1 + x2)/2 + (x2 - x1)/4
368C + 0.25 (x2 - x1)^2 / (x2 - x1)
369C si x2 > x1 => x = x2
370C si x2 < x1 => x = x1
371C
372C et x < x3 + 1/4 |x2 - x1|
373C
374C si x3 = x2
375C x = x3 + 0.125 (x1 - x2)^2 / (2 x3 - x1 - x2 )
376C x = x3 + 0.125 (x1 - x2)^2 / (x2 - x1)
377C x = x3 +- 1/8 |x1 - x2| = x2 +- 1/8 |x1 - x2|
378C
379C si (x1 + x2)/2 + (x2 - x1)/4 < x3 < max(x1,x2)
380C x < max(x1,x2) + 1/8 |x1 - x2|
381C et x < x3 + 1/4 |x1 - x2|
382C
383C si x3 > max(x1,x2)
384C x < x3 + 1/8 |x1 - x2|
385C------------------------------------------------------------
386C solution 2 => x < max (x1 , x2 , x3) + 1/8 |x2 - x1|
387C------------------------------------------------------------
388C solution 1 : x < max (x1 , x2 , x3 + 1/4 |x2 - x1|)
389C------------------------------------------------------------
390C => x <min (solution 1, solution 2)
391C------------------------------------------------------------
392C
393C 3) Exact solution (0 terminal solution by solution 1)
394C
395C solution 0 :
396C x = x3 - 0.125 (x1 - x2)^2 / (x1 + x2 - 2 x3)
397C solution 1 :
398C x < max (x1 , x2 , x3 + 1/4 |x2 - x1|)
399C------------------------------------------------------------
400C Solution 3 => x <min (solution 1, solution 0)
401C------------------------------------------------------------
402C s = (x1+x2)/2 d = |x2-x1|
403C-------------------------------------------------------------------------
404C x3 x min x max
405C-------------------------------------------------------------------------
406C -inf < x3 < s - d/4 x3 + d^2 / 16(x3-s) max(x1,x2)
407C s - d/4 < x3 < s + d/4 min(x1,x2) max(x1,x2)
408C s + d/4 < x3 < +inf min(x1,x2) x3 + d^2 / 16(x3-s)
409C-------------------------------------------------------------------------
410C---------------------------------------------------------
411C BIDIM ( pas resolu )
412C---------------------------------------------------------
413C
414C-----------------------------------------------
415C i ri si ti Ni
416C--------------------------------------------------------------------
417C 1 -1 -1 -1 1/4(1-r)(1-t)(-r-t-1)
418C 2 -1 -1 +1 1/4(1-r)(1+t)(-r+t-1)
419C 3 +1 -1 +1 1/4(1+r)(1+t)(+r+t-1)
420C 4 +1 -1 -1 1/4(1+r)(1-t)(+r-t-1)
421C 9 -1 -1 0 1/2(1-t^2)(1-r)
422C 10 0 -1 +1 1/2(1-r^2) (1+t)
423C 11 +1 -1 0 1/2(1-t^2)(1+r)
424C 12 0 -1 -1 1/2(1-r^2) (1-t)
425C
426C x = N1 x1 + N2 x2 + N3 x3 + N4 x4
427C + N9 x9 + N10 x10 + N11 x11 + N12 x12
428C
429C 0) Point XMAX
430C
431C dx/dr = -1/4(1-t)(-2r-t) x1
432C -1/4(1+t)(-2r+t) x2
433C +1/4(1+t)(+2r+t) x3
434C +1/4(1-t)(+2r-t) x4
435C -1/2(1-t^2) x9
436C -r(1+t) x10
437C +1/2(1-t^2) x11
438C -r(1-t) x12 = 0
439C dx/dt = -1/4(1-r)(-2t-r) x1
440C +1/4(1-r)(+2t-r) x2
441C +1/4(1+r)(+2t+r) x3
442C +1/4(1+r)(-2t+r) x4
443C -t(1-r) x9
444C +1/2(1-r^2) x10
445C -t(1+r) x11
446C -1/2(1-r^2) x12 = 0
447C------------------------------------
448C calculation of element bounds
449C------------------------------------
450 DO iface=1,2
451C-----------------------------------------------------------------------
452C Face 1 2 3 4 or 5 6 7 8
453C-----------------------------------------------------------------------
454 DO idir=1,3
455C-----------------------------------------------------------------------
456C x y or z
457C-----------------------------------------------------------------------
458 DO l=lft,llt
459 i = index(l)
460 ne = nelem(i)
461 n16= ne - numels8 - numels10 - numels20
462C
463 j = ixs(iperm(1,iface),ne)
464 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
465 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
466 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
467 j = ixs(iperm(2,iface),ne)
468 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
469 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
470 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
471 j = ixs(iperm(3,iface),ne)
472 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
473 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
474 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
475 j = ixs(iperm(4,iface),ne)
476 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
477 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
478 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
479 j = ixs16(iperm(5,iface),n16)
480 x9 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
481 xmsr(idir) =max(xmsr(idir) ,x9-xsav(idir,j))
482 xmsr(idir+3)=min(xmsr(idir+3),x9-xsav(idir,j))
483 j = ixs16(iperm(6,iface),n16)
484 x10= x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
485 xmsr(idir) =max(xmsr(idir) ,x10-xsav(idir,j))
486 xmsr(idir+3)=min(xmsr(idir+3),x10-xsav(idir,j))
487 j = ixs16(iperm(7,iface),n16)
488 x11= x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
489 xmsr(idir) =max(xmsr(idir) ,x11-xsav(idir,j))
490 xmsr(idir+3)=min(xmsr(idir+3),x11-xsav(idir,j))
491 j = ixs16(iperm(8,iface),n16)
492 x12= x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
493 xmsr(idir) =max(xmsr(idir) ,x12-xsav(idir,j))
494 xmsr(idir+3)=min(xmsr(idir+3),x12-xsav(idir,j))
495C
496 xc = half*(x9+x10+x11+x12) - fourth*(x1+x2+x3+x4)
497C
498 d4 = fourth * abs(x1-x2)
499 an = min( x1 , x2 , x9-d4 )
500 ax = max( x1 , x2 , x9+d4 )
501C
502 d4 = fourth * abs(x3-x4)
503 bn = min( x3 , x4 , x11-d4 )
504 bx = max( x3 , x4 , x11+d4 )
505C
506 d4 = fourth * abs(x12-x10)
507 cn = min( x12 , x10 , xc-d4 )
508 cx = max( x12 , x10 , xc+d4 )
509C
510 d8 = one_over_8 * max( ax-bn , bx-an )
511 d4 = d8 + d8
512 dn = max(min( an , bn , cn-d4 ),min(an , bn , cn) - d8 )
513 dx = min(max( ax , bx , cx+d4 ),max( ax , bx , cx) + d8 )
514C
515 eminx(idir,i) = min( eminx(idir,i) , dn )
516 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
517C
518 SIZE = SIZE + dx - dn
519C
520 ENDDO
521 ENDDO
522 ENDDO
523C--------------------------------------------------------------
524C
525 RETURN
integer function iface(ip, n)
Definition iface.F:36

◆ i16crit()

subroutine i16crit ( x,
integer, dimension(*) nsv,
integer, dimension(*) nelem,
integer nsn,
eminx,
integer nme,
integer itask,
xsav,
integer, dimension(nixs,*) ixs,
integer, dimension(8,*) ixs16,
integer, dimension(12,*) ixs20,
integer, dimension(6,*) ixs10,
v,
a,
dimension(7) xmsrg,
dimension(7) xslvg )

Definition at line 36 of file i16crit.F.

41C-----------------------------------------------
42 USE intbufdef_mod
43 use element_mod , only : nixs
44C-----------------------------------------------
45C I m p l i c i t T y p e s
46C-----------------------------------------------
47#include "implicit_f.inc"
48#include "comlock.inc"
49C-----------------------------------------------
50C G l o b a l P a r a m e t e r s
51C-----------------------------------------------
52#include "mvsiz_p.inc"
53C-----------------------------------------------
54C C o m m o n B l o c k s
55C-----------------------------------------------
56#include "com04_c.inc"
57#include "com08_c.inc"
58#include "task_c.inc"
59 COMMON /i16tmp/SIZE
60 my_real
61 . SIZE
62C-----------------------------------------------
63C D u m m y A r g u m e n t s
64C-----------------------------------------------
65 my_real, DIMENSION(7) :: xmsrg
66 my_real, DIMENSION(7) :: xslvg
67 INTEGER NSN,ITASK,NME,
68 . NSV(*),NELEM(*),IXS(NIXS,*),IXS16(8,*),IXS20(12,*),
69 . IXS10(6,*)
71 . x(3,*),v(3,*),a(3,*),xsav(3,*),eminx(6,*)
72C-----------------------------------------------
73C L o c a l V a r i a b l e s
74C-----------------------------------------------
75 INTEGER NSNF,NMEF,NSNL,NMEL,I, J,I16,I20,LFT16,LLT16,
76 . LFT20,LLT20,LFT8,LLT8,LFT10,LLT10,I8,I10,
77 . INDEX16(MVSIZ),INDEX20(MVSIZ),INDEX8(MVSIZ),INDEX10(MVSIZ)
79 . xmsr(6),xslv(6),size_t ,xx,yy,zz
80C-----------------------------------------------
81C S o u r c e L i n e s
82C-----------------------------------------------
83 nsnf = 1 + itask*nsn / nthread
84 nsnl = (itask+1)*nsn / nthread
85 nmef = 1 + itask*nme / nthread
86 nmel = (itask+1)*nme / nthread
87C--------------------------------------------------------------
88C 0- criterion calculation to know if we need to sort or not
89C--------------------------------------------------------------
90 xslv(1) = -ep30
91 xslv(2) = -ep30
92 xslv(3) = -ep30
93 xslv(4) = ep30
94 xslv(5) = ep30
95 xslv(6) = ep30
96 xmsr(1) = -ep30
97 xmsr(2) = -ep30
98 xmsr(3) = -ep30
99 xmsr(4) = ep30
100 xmsr(5) = ep30
101 xmsr(6) = ep30
102C
103 size_t = zero
104C
105 DO i=nsnf,nsnl
106 j=nsv(i)
107 xx=x(1,j)+dt2*(v(1,j)+dt12*a(1,j))
108 yy=x(2,j)+dt2*(v(2,j)+dt12*a(2,j))
109 zz=x(3,j)+dt2*(v(3,j)+dt12*a(3,j))
110 xslv(1)=max(xslv(1),xx-xsav(1,j))
111 xslv(2)=max(xslv(2),yy-xsav(2,j))
112 xslv(3)=max(xslv(3),zz-xsav(3,j))
113 xslv(4)=min(xslv(4),xx-xsav(1,j))
114 xslv(5)=min(xslv(5),yy-xsav(2,j))
115 xslv(6)=min(xslv(6),zz-xsav(3,j))
116 END DO
117C------------------------------------
118C calculation of element bounds
119C------------------------------------
120 DO i=nmef,nmel
121 eminx(1,i) = ep30
122 eminx(2,i) = ep30
123 eminx(3,i) = ep30
124 eminx(4,i) = -ep30
125 eminx(5,i) = -ep30
126 eminx(6,i) = -ep30
127 ENDDO
128C
129 lft16=1
130 llt16=0
131 lft20=1
132 llt20=0
133 lft8 =1
134 llt8 =0
135 lft10=1
136 llt10=0
137 DO i=nmef,nmel
138 i8 = nelem(i)
139 i10 = i8-numels8
140 i20 = i10-numels10
141 i16 = i20-numels20
142 IF(i16>=1.AND.i16<=numels16)THEN
143 llt16=llt16+1
144 index16(llt16)=i
145 IF(llt16==mvsiz-1)THEN
146 CALL i16box(
147 1 lft16,llt16 ,nelem,eminx,nmef ,nmel ,
148 2 x ,v ,a ,ixs ,ixs16,size_t,
149 3 xmsr ,index16,xsav )
150 llt16=0
151 ENDIF
152 ELSEIF(i20>=1.AND.i20<=numels20)THEN
153 llt20=llt20+1
154 index20(llt20)=i
155 IF(llt20==mvsiz-1)THEN
156 CALL i20box(
157 1 lft20,llt20 ,nelem,eminx,nmef ,nmel ,
158 2 x ,v ,a ,ixs ,ixs20,size_t,
159 3 xmsr ,index20,xsav )
160 llt20=0
161 ENDIF
162 ELSEIF(i10>=1)THEN
163 llt10=llt10+1
164 index10(llt10)=i
165 IF(llt10==mvsiz-1)THEN
166 CALL i10box(
167 1 lft10,llt10 ,nelem,eminx,nmef ,nmel ,
168 2 x ,v ,a ,ixs ,ixs10,size_t,
169 3 xmsr ,index10,xsav )
170 llt10=0
171 ENDIF
172 ELSEIF(i8>=1)THEN
173 llt8=llt8+1
174 index8(llt8)=i
175 IF(llt8==mvsiz-1)THEN
176 CALL i8box(
177 1 lft8 ,llt8 ,nelem,eminx,nmef ,nmel ,
178 2 x ,v ,a ,ixs ,size_t,
179 3 xmsr ,index8 ,xsav )
180 llt8=0
181 ENDIF
182 ENDIF
183 END DO
184 IF(llt16>0)CALL i16box(
185 1 lft16,llt16 ,nelem,eminx,nmef ,nmel ,
186 2 x ,v ,a ,ixs ,ixs16,size_t,
187 3 xmsr ,index16,xsav )
188 IF(llt20>0)CALL i20box(
189 1 lft20,llt20 ,nelem,eminx,nmef ,nmel ,
190 2 x ,v ,a ,ixs ,ixs20,size_t,
191 3 xmsr ,index20,xsav )
192 IF(llt8>0)CALL i8box(
193 1 lft8 ,llt8 ,nelem,eminx,nmef ,nmel ,
194 2 x ,v ,a ,ixs ,size_t,
195 3 xmsr ,index8 ,xsav )
196 IF(llt10>0)CALL i10box(
197 1 lft10,llt10 ,nelem,eminx,nmef ,nmel ,
198 2 x ,v ,a ,ixs ,ixs10,size_t,
199 3 xmsr ,index10,xsav )
200C
201#include "lockon.inc"
202 xslvg(1)=max(xslvg(1),xslv(1))
203 xslvg(2)=max(xslvg(2),xslv(2))
204 xslvg(3)=max(xslvg(3),xslv(3))
205 xslvg(4)=min(xslvg(4),xslv(4))
206 xslvg(5)=min(xslvg(5),xslv(5))
207 xslvg(6)=min(xslvg(6),xslv(6))
208 xmsrg(1)=max(xmsrg(1),xmsr(1))
209 xmsrg(2)=max(xmsrg(2),xmsr(2))
210 xmsrg(3)=max(xmsrg(3),xmsr(3))
211 xmsrg(4)=min(xmsrg(4),xmsr(4))
212 xmsrg(5)=min(xmsrg(5),xmsr(5))
213 xmsrg(6)=min(xmsrg(6),xmsr(6))
214 SIZE = SIZE + size_t
215#include "lockoff.inc"
216C
217 RETURN
subroutine i16box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs16, size, xmsr, index, xsav)
Definition i16crit.F:230
subroutine i20box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs20, size, xmsr, index, xsav)
Definition i16crit.F:538
subroutine i10box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, ixs10, size, xmsr, index, xsav)
Definition i16crit.F:868
subroutine i8box(lft, llt, nelem, eminx, nmef, nmel, x, v, a, ixs, size, xmsr, index, xsav)
Definition i16crit.F:992

◆ i20box()

subroutine i20box ( integer lft,
integer llt,
integer, dimension(*) nelem,
eminx,
integer nmef,
integer nmel,
x,
v,
a,
integer, dimension(nixs,*) ixs,
integer, dimension(12,*) ixs20,
size,
xmsr,
integer, dimension(*) index,
xsav )

Definition at line 535 of file i16crit.F.

538 use element_mod , only : nixs
539C-----------------------------------------------
540C I m p l i c i t T y p e s
541C-----------------------------------------------
542#include "implicit_f.inc"
543C-----------------------------------------------
544C C o m m o n B l o c k s
545C-----------------------------------------------
546#include "com04_c.inc"
547#include "com08_c.inc"
548C-----------------------------------------------
549C D u m m y A r g u m e n t s
550C-----------------------------------------------
551 INTEGER LFT ,LLT,NMEF,NMEL,
552 . IXS(NIXS,*),IXS20(12,*),NELEM(*),INDEX(*)
553C REAL
554 my_real
555 . x(3,*),v(3,*),a(3,*),eminx(6,*),SIZE,xmsr(*),xsav(3,*)
556C-----------------------------------------------
557C L o c a l V a r i a b l e s
558C-----------------------------------------------
559 INTEGER I,J,L,NE,IDIR,N20
560 my_real
561 . an12,ax12,an34,ax34,an56,ax56,an78,ax78,cn,cx,dx,dn,d4,d8,
562 . x1,x2,x3,x4,x5,x6,x7,x8,
563 . x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,xc
564C------------------------------------
565C calculation of element bounds
566C------------------------------------
567C-----------------------------------------------------------------------
568C Face 1 2 3 4 or 5 6 7 8
569C-----------------------------------------------------------------------
570 DO idir=1,3
571C-----------------------------------------------------------------------
572C x y or z
573C-----------------------------------------------------------------------
574 DO l=lft,llt
575 i = index(l)
576 ne = nelem(i)
577 n20= ne - numels8 - numels10
578C-----------------------------------------------------------------------
579 j = ixs(2,ne)
580 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
581 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
582 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
583 j = ixs(3,ne)
584 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
585 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
586 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
587 j = ixs(4,ne)
588 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
589 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
590 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
591 j = ixs(5,ne)
592 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
593 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
594 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
595 j = ixs(6,ne)
596 x5 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
597 xmsr(idir) =max(xmsr(idir) ,x5-xsav(idir,j))
598 xmsr(idir+3)=min(xmsr(idir+3),x5-xsav(idir,j))
599 j = ixs(7,ne)
600 x6 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
601 xmsr(idir) =max(xmsr(idir) ,x6-xsav(idir,j))
602 xmsr(idir+3)=min(xmsr(idir+3),x6-xsav(idir,j))
603 j = ixs(8,ne)
604 x7 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
605 xmsr(idir) =max(xmsr(idir) ,x7-xsav(idir,j))
606 xmsr(idir+3)=min(xmsr(idir+3),x7-xsav(idir,j))
607 j = ixs(9,ne)
608 x8 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
609 xmsr(idir) =max(xmsr(idir) ,x8-xsav(idir,j))
610 xmsr(idir+3)=min(xmsr(idir+3),x8-xsav(idir,j))
611C
612 j = ixs20(1,n20)
613 IF(j/=0)THEN
614 x9 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
615 xmsr(idir) =max(xmsr(idir) ,x9-xsav(idir,j))
616 xmsr(idir+3)=min(xmsr(idir+3),x9-xsav(idir,j))
617 ELSE
618 x9=0.5*(x(idir,ixs(2,ne))+x(idir,ixs(3,ne)))
619 ENDIF
620 j = ixs20(2,n20)
621 IF(j/=0)THEN
622 x10 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
623 xmsr(idir) =max(xmsr(idir) ,x10-xsav(idir,j))
624 xmsr(idir+3)=min(xmsr(idir+3),x10-xsav(idir,j))
625 ELSE
626 x10=0.5*(x(idir,ixs(3,ne))+x(idir,ixs(4,ne)))
627 ENDIF
628 j = ixs20(3,n20)
629 IF(j/=0)THEN
630 x11 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
631 xmsr(idir) =max(xmsr(idir) ,x11-xsav(idir,j))
632 xmsr(idir+3)=min(xmsr(idir+3),x11-xsav(idir,j))
633 ELSE
634 x11=0.5*(x(idir,ixs(4,ne))+x(idir,ixs(5,ne)))
635 ENDIF
636 j = ixs20(4,n20)
637 IF(j/=0)THEN
638 x12 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
639 xmsr(idir) =max(xmsr(idir) ,x12-xsav(idir,j))
640 xmsr(idir+3)=min(xmsr(idir+3),x12-xsav(idir,j))
641 ELSE
642 x12=0.5*(x(idir,ixs(5,ne))+x(idir,ixs(2,ne)))
643 ENDIF
644 j = ixs20(5,n20)
645 IF(j/=0)THEN
646 x13 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
647 xmsr(idir) =max(xmsr(idir) ,x13-xsav(idir,j))
648 xmsr(idir+3)=min(xmsr(idir+3),x13-xsav(idir,j))
649 ELSE
650 x13=0.5*(x(idir,ixs(2,ne))+x(idir,ixs(6,ne)))
651 ENDIF
652 j = ixs20(6,n20)
653 IF(j/=0)THEN
654 x14 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
655 xmsr(idir) =max(xmsr(idir) ,x14-xsav(idir,j))
656 xmsr(idir+3)=min(xmsr(idir+3),x14-xsav(idir,j))
657 ELSE
658 x14=0.5*(x(idir,ixs(3,ne))+x(idir,ixs(6,ne)))
659 ENDIF
660 j = ixs20(7,n20)
661 IF(j/=0)THEN
662 x15 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
663 xmsr(idir) =max(xmsr(idir) ,x15-xsav(idir,j))
664 xmsr(idir+3)=min(xmsr(idir+3),x15-xsav(idir,j))
665 ELSE
666 x15=0.5*(x(idir,ixs(4,ne))+x(idir,ixs(8,ne)))
667 ENDIF
668 j = ixs20(8,n20)
669 IF(j/=0)THEN
670 x16 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
671 xmsr(idir) =max(xmsr(idir) ,x16-xsav(idir,j))
672 xmsr(idir+3)=min(xmsr(idir+3),x16-xsav(idir,j))
673 ELSE
674 x16=0.5*(x(idir,ixs(5,ne))+x(idir,ixs(9,ne)))
675 ENDIF
676 j = ixs20(9,n20)
677 IF(j/=0)THEN
678 x17 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
679 xmsr(idir) =max(xmsr(idir) ,x17-xsav(idir,j))
680 xmsr(idir+3)=min(xmsr(idir+3),x17-xsav(idir,j))
681 ELSE
682 x17=0.5*(x(idir,ixs(6,ne))+x(idir,ixs(7,ne)))
683 ENDIF
684 j = ixs20(10,n20)
685 IF(j/=0)THEN
686 x18 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
687 xmsr(idir) =max(xmsr(idir) ,x18-xsav(idir,j))
688 xmsr(idir+3)=min(xmsr(idir+3),x18-xsav(idir,j))
689 ELSE
690 x18=0.5*(x(idir,ixs(7,ne))+x(idir,ixs(8,ne)))
691 ENDIF
692 j = ixs20(11,n20)
693 IF(j/=0)THEN
694 x19 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
695 xmsr(idir) =max(xmsr(idir) ,x19-xsav(idir,j))
696 xmsr(idir+3)=min(xmsr(idir+3),x19-xsav(idir,j))
697 ELSE
698 x19=0.5*(x(idir,ixs(8,ne))+x(idir,ixs(9,ne)))
699 ENDIF
700 j = ixs20(12,n20)
701 IF(j/=0)THEN
702 x20 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
703 xmsr(idir) =max(xmsr(idir) ,x20-xsav(idir,j))
704 xmsr(idir+3)=min(xmsr(idir+3),x20-xsav(idir,j))
705 ELSE
706 x20=0.5*(x(idir,ixs(6,ne))+x(idir,ixs(9,ne)))
707 ENDIF
708C
709C-----------------------------------------------------------------------
710C Face 1 2 3 4
711C-----------------------------------------------------------------------
712 xc = half*(x9+x10+x11+x12) - fourth*(x1+x2+x3+x4)
713C
714 d4 = fourth * abs(x1-x2)
715 an12 = min( x1 , x2 , x9-d4 )
716 ax12 = max( x1 , x2 , x9+d4 )
717C
718 d4 = fourth * abs(x3-x4)
719 an34 = min( x3 , x4 , x11-d4 )
720 ax34 = max( x3 , x4 , x11+d4 )
721C
722 d4 = fourth * abs(x12-x10)
723 cn = min( x12 , x10 , xc-d4 )
724 cx = max( x12 , x10 , xc+d4 )
725C
726 d8 = one_over_8 * max( ax12-an34 , ax34-an12 )
727 d4 = d8 + d8
728 dn = max(min(an12 , an34 , cn-d4 ),
729 . min(an12 , an34 , cn) - d8 )
730 dx = min(max(ax12 , ax34 , cx+d4 ),
731 . max(ax12 , ax34 , cx) + d8 )
732C
733 eminx(idir,i) = min( eminx(idir,i) , dn )
734 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
735C-----------------------------------------------------------------------
736C Face 5 6 7 8
737C-----------------------------------------------------------------------
738 xc = half*(x17+x18+x19+x20) - fourth*(x5+x6+x7+x8)
739C
740 d4 = fourth * abs(x5-x6)
741 an56 = min( x5 , x6 , x17-d4 )
742 ax56 = max( x5 , x6 , x17+d4 )
743C
744 d4 = fourth * abs(x7-x8)
745 an78 = min( x7 , x8 , x19-d4 )
746 ax78 = max( x7 , x8 , x19+d4 )
747C
748 d4 = fourth * abs(x20-x18)
749 cn = min( x20 , x18 , xc-d4 )
750 cx = max( x20 , x18 , xc+d4 )
751C
752 d8 = one_over_8 * max( ax56-an78 , ax78-an56 )
753 d4 = d8 + d8
754 dn = max(min(an56 , an78 , cn-d4 ),
755 . min(an56 , an78 , cn) - d8 )
756 dx = min(max(ax56 , ax78 , cx+d4 ),
757 . max(ax56 , ax78 , cx) + d8 )
758C
759 eminx(idir,i) = min( eminx(idir,i) , dn )
760 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
761C-----------------------------------------------------------------------
762C Face 1 2 6 5
763C-----------------------------------------------------------------------
764 xc = half*(x9+x14+x17+x13) - fourth*(x1+x2+x6+x5)
765C
766 d4 = fourth * abs(x13-x14)
767 cn = min( x13 , x14 , xc-d4 )
768 cx = max( x13 , x14 , xc+d4 )
769C
770 d8 = one_over_8 * max( ax12-an56 , ax56-an12 )
771 d4 = d8 + d8
772 dn = max(min(an12 , an56 , cn-d4 ),
773 . min(an12 , an56 , cn) - d8 )
774 dx = min(max(ax12 , ax56 , cx+d4 ),
775 . max(ax12 , ax56 , cx) + d8 )
776C
777 eminx(idir,i) = min( eminx(idir,i) , dn )
778 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
779C-----------------------------------------------------------------------
780C Face 3 4 8 7
781C-----------------------------------------------------------------------
782 xc = half*(x11+x15+x19+x16) - fourth*(x3+x4+x8+x7)
783C
784 d4 = fourth * abs(x16-x15)
785 cn = min( x15 , x16 , xc-d4 )
786 cx = max( x15 , x16 , xc+d4 )
787C
788 d8 = one_over_8 * max( ax34-an78 , ax78-an34 )
789 d4 = d8 + d8
790 dn = max(min(an34 , an78 , cn-d4 ),
791 . min(an34 , an78 , cn) - d8 )
792 dx = min(max(ax34 , ax78 , cx+d4 ),
793 . max(ax34 , ax78 , cx) + d8 )
794C
795 eminx(idir,i) = min( eminx(idir,i) , dn )
796 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
797C-----------------------------------------------------------------------
798C Face 4 1 5 8
799C-----------------------------------------------------------------------
800 xc = half*(x12+x13+x20+x16) - fourth*(x4+x1+x5+x8)
801C
802 d4 = fourth * abs(x4-x1)
803 an12 = min( x4 , x1 , x12-d4 )
804 ax12 = max( x4 , x1 , x12+d4 )
805C
806 d4 = fourth * abs(x8-x5)
807 an34 = min( x8 , x5 , x20-d4 )
808 ax34 = max( x8 , x5 , x20+d4 )
809C
810 d4 = fourth * abs(x16-x13)
811 cn = min( x16 , x13 , xc-d4 )
812 cx = max( x16 , x13 , xc+d4 )
813C
814 d8 = one_over_8 * max( ax12-an34 , ax34-an12 )
815 d4 = d8 + d8
816 dn = max(min(an12 , an34 , cn-d4 ),
817 . min(an12 , an34 , cn) - d8 )
818 dx = min(max(ax12 , ax34 , cx+d4 ),
819 . max(ax12 , ax34 , cx) + d8 )
820C
821 eminx(idir,i) = min( eminx(idir,i) , dn )
822 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
823C-----------------------------------------------------------------------
824C Face 3 2 6 7
825C-----------------------------------------------------------------------
826 xc = half*(x10+x14+x18+x15) - fourth*(x3+x2+x6+x7)
827C
828 d4 = fourth * abs(x3-x2)
829 an12 = min( x3 , x2 , x10-d4 )
830 ax12 = max( x3 , x2 , x10+d4 )
831C
832 d4 = fourth * abs(x7-x6)
833 an34 = min( x7 , x6 , x18-d4 )
834 ax34 = max( x7 , x6 , x18+d4 )
835C
836 d4 = fourth * abs(x15-x14)
837 cn = min( x15 , x14 , xc-d4 )
838 cx = max( x15 , x14 , xc+d4 )
839C
840 d8 = one_over_8* max( ax12-an34 , ax34-an12 )
841 d4 = d8 + d8
842 dn = max(min(an12 , an34 , cn-d4 ),
843 . min(an12 , an34 , cn) - d8 )
844 dx = min(max(ax12 , ax34 , cx+d4 ),
845 . max(ax12 , ax34 , cx) + d8 )
846C
847 eminx(idir,i) = min( eminx(idir,i) , dn )
848 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
849C-----------------------------------------------------------------------
850 SIZE = SIZE + dx - dn
851C
852 ENDDO
853 ENDDO
854C--------------------------------------------------------------
855C
856 RETURN

◆ i8box()

subroutine i8box ( integer lft,
integer llt,
integer, dimension(*) nelem,
eminx,
integer nmef,
integer nmel,
x,
v,
a,
integer, dimension(nixs,*) ixs,
size,
xmsr,
integer, dimension(*) index,
xsav )

Definition at line 989 of file i16crit.F.

992 use element_mod , only : nixs
993C-----------------------------------------------
994C I m p l i c i t T y p e s
995C-----------------------------------------------
996#include "implicit_f.inc"
997C-----------------------------------------------
998C C o m m o n B l o c k s
999C-----------------------------------------------
1000#include "com08_c.inc"
1001C-----------------------------------------------
1002C D u m m y A r g u m e n t s
1003C-----------------------------------------------
1004 INTEGER LFT ,LLT,NMEF,NMEL,
1005 . IXS(NIXS,*),NELEM(*),INDEX(*)
1006C REAL
1007 my_real
1008 . x(3,*),v(3,*),a(3,*),eminx(6,*),SIZE,xmsr(*),xsav(3,*)
1009C-----------------------------------------------
1010C L o c a l V a r i a b l e s
1011C-----------------------------------------------
1012 INTEGER I,J,L,NE,IDIR
1013 my_real
1014 . dx,dn,
1015 . x1,x2,x3,x4,x5,x6,x7,x8
1016C------------------------------------
1017C calculation of element bounds
1018C------------------------------------
1019C-----------------------------------------------------------------------
1020C Face 1 2 3 4 or 5 6 7 8
1021C-----------------------------------------------------------------------
1022 DO idir=1,3
1023C-----------------------------------------------------------------------
1024C x y or z
1025C-----------------------------------------------------------------------
1026 DO l=lft,llt
1027 i = index(l)
1028 ne = nelem(i)
1029C-----------------------------------------------------------------------
1030 j = ixs(2,ne)
1031 x1 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1032 xmsr(idir) =max(xmsr(idir) ,x1-xsav(idir,j))
1033 xmsr(idir+3)=min(xmsr(idir+3),x1-xsav(idir,j))
1034 j = ixs(3,ne)
1035 x2 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1036 xmsr(idir) =max(xmsr(idir) ,x2-xsav(idir,j))
1037 xmsr(idir+3)=min(xmsr(idir+3),x2-xsav(idir,j))
1038 j = ixs(4,ne)
1039 x3 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1040 xmsr(idir) =max(xmsr(idir) ,x3-xsav(idir,j))
1041 xmsr(idir+3)=min(xmsr(idir+3),x3-xsav(idir,j))
1042 j = ixs(5,ne)
1043 x4 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1044 xmsr(idir) =max(xmsr(idir) ,x4-xsav(idir,j))
1045 xmsr(idir+3)=min(xmsr(idir+3),x4-xsav(idir,j))
1046 j = ixs(6,ne)
1047 x5 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1048 xmsr(idir) =max(xmsr(idir) ,x5-xsav(idir,j))
1049 xmsr(idir+3)=min(xmsr(idir+3),x5-xsav(idir,j))
1050 j = ixs(7,ne)
1051 x6 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1052 xmsr(idir) =max(xmsr(idir) ,x6-xsav(idir,j))
1053 xmsr(idir+3)=min(xmsr(idir+3),x6-xsav(idir,j))
1054 j = ixs(8,ne)
1055 x7 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1056 xmsr(idir) =max(xmsr(idir) ,x7-xsav(idir,j))
1057 xmsr(idir+3)=min(xmsr(idir+3),x7-xsav(idir,j))
1058 j = ixs(9,ne)
1059 x8 = x(idir,j)+dt2*(v(idir,j)+dt12*a(idir,j))
1060 xmsr(idir) =max(xmsr(idir) ,x8-xsav(idir,j))
1061 xmsr(idir+3)=min(xmsr(idir+3),x8-xsav(idir,j))
1062C
1063
1064 dx=max(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 )
1065 dn=min(x1,x2 ,x3 ,x4 ,x5 ,x6 ,x7 ,x8 )
1066C
1067 eminx(idir,i) = min( eminx(idir,i) , dn )
1068 eminx(idir+3,i) = max( eminx(idir+3,i), dx )
1069C
1070 SIZE = SIZE + dx - dn
1071C
1072 ENDDO
1073 ENDDO
1074C--------------------------------------------------------------
1075C
1076 RETURN