538 use element_mod , only : nixs
539
540
541
542#include "implicit_f.inc"
543
544
545
546#include "com04_c.inc"
547#include "com08_c.inc"
548
549
550
551 INTEGER ,LLT,NMEF,NMEL,
552 . IXS(NIXS,*),IXS20(12,*),NELEM(*),INDEX(*)
553
555 . x(3,*),v(3,*),a(3,*),eminx(6,*),SIZE,xmsr(*),xsav(3,*)
556
557
558
559 INTEGER I,J,L,NE,IDIR,N20
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
564
565
566
567
568
569
570 DO idir
571
572
573
574 DO l=lft,llt
575 i = index(l)
576 ne = nelem(i)
577 n20= ne - numels8 - numels10
578
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))
611
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
708
709
710
711
712 xc = half*(x9+x10+x11+x12) - fourth*(x1+x2+x3+x4)
713
714 d4 = fourth * abs(x1-x2)
715 an12 =
min( x1 , x2 , x9-d4 )
716 ax12 =
max( x1 , x2 , x9+d4 )
717
718 d4 = fourth * abs(x3-x4)
719 an34 =
min( x3 , x4 , x11-d4 )
720 ax34 =
max( x3 , x4 , x11+d4 )
721
722 d4 = fourth * abs(x12-x10)
723 cn =
min( x12 , x10 , xc-d4 )
724 cx =
max( x12 , x10 , xc+d4 )
725
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 )
732
733 eminx(idir,i) =
min( eminx(idir,i) , dn )
734 eminx(idir+3,i) =
max( eminx(idir+3,i), dx )
735
736
737
738 xc = half*(x17+x18+x19+x20) - fourth*(x5+x6+x7+x8)
739
740 d4 = fourth * abs(x5-x6)
741 an56 =
min( x5 , x6 , x17-d4 )
742 ax56 =
max( x5 , x6 , x17+d4 )
743
744 d4 = fourth * abs(x7-x8)
745 an78 =
min( x7 , x8 , x19-d4 )
746 ax78 =
max( x7 , x8 , x19+d4 )
747
748 d4 = fourth * abs(x20-x18)
749 cn =
min( x20 , x18 , xc-d4 )
750 cx =
max( x20 , x18 , xc+d4 )
751
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 )
758
759 eminx(idir,i) =
min( eminx(idir,i) , dn )
760 eminx(idir+3,i) =
max( eminx(idir+3,i), dx )
761
762
763
764 xc = half*(x9+x14+x17+x13) - fourth*(x1+x2+x6+x5)
765
766 d4 = fourth * abs(x13-x14)
767 cn =
min( x13 , x14 , xc-d4 )
768 cx =
max( x13 , x14 , xc+d4 )
769
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 )
776
777 eminx(idir,i) =
min( eminx(idir,i) , dn )
778 eminx(idir+3,i) =
max( eminx(idir+3,i), dx )
779
780
781
782 xc = half*(x11+x15+x19+x16) - fourth*(x3+x4+x8+x7)
783
784 d4 = fourth * abs(x16-x15)
785 cn =
min( x15 , x16 , xc-d4 )
786 cx =
max( x15 , x16 , xc+d4 )
787
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 )
794
795 eminx(idir,i) =
min( eminx(idir,i) , dn )
796 eminx(idir+3,i) =
max( eminx(idir+3,i), dx )
797
798
799
800 xc = half*(x12+x13+x20+x16) - fourth*(x4+x1+x5+x8)
801
802 d4 = fourth * abs(x4-x1)
803 an12 =
min( x4 , x1 , x12-d4 )
804 ax12 =
max( x4 , x1 , x12+d4 )
805
806 d4 = fourth * abs(x8-x5)
807 an34 =
min( x8 , x5 , x20-d4 )
808 ax34 =
max( x8 , x5 , x20+d4 )
809
810 d4 = fourth * abs(x16-x13)
811 cn =
min( x16 , x13 , xc-d4 )
812 cx =
max( x16 , x13 , xc+d4 )
813
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 )
820
821 eminx(idir,i) =
min( eminx(idir,i) , dn )
822 eminx(idir+3,i) =
max( eminx(idir+3,i), dx )
823
824
825
826 xc = half*(x10+x14+x18+x15) - fourth*(x3+x2+x6+x7)
827
828 d4 = fourth * abs(x3-x2)
829 an12 =
min( x3 , x2 , x10-d4 )
830 ax12 =
max( x3 , x2 , x10+d4 )
831
832 d4 = fourth * abs(x7-x6)
833 an34 =
min( x7 , x6 , x18-d4 )
834 ax34 =
max( x7 , x6 , x18+d4 )
835
836 d4 = fourth * abs(x15-x14)
837 cn =
min( x15 , x14 , xc-d4 )
838 cx =
max( x15 , x14 , xc+d4 )
839
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 )
846
847 eminx(idir,i) =
min( eminx(idir,i) , dn )
848 eminx(idir+3,i) =
max( eminx(idir+3,i), dx )
849
850 SIZE = SIZE + dx - dn
851
852 ENDDO
853 ENDDO
854
855
856 RETURN