546
547
548
549#include "implicit_f.inc"
550
551
552
553#include "param_c.inc"
554
555
556
557 INTEGER INRBE3(*),ILRBE3(*),NG, NS,IROT,IMODIF,IERR,IPEN
558
560 . xyz(3,*), frbe3(6,*), skew(lskew,*),wmin
561
562
563
564 INTEGER I, J, K,N, ,NML, IAD,JJ,KG,,IELSUB,KDIAG
565
567 * tw(3,ng), rw(3,ng),
568 * fufxlc(3,ng), fufylc(3,ng), fufzlc(3,ng),
569 * fumxlc(3,ng), fumylc(3,ng), fumzlc(3,ng),
570 * mxlc(3,ng), mylc(3,ng), mzlc(3,ng),
571 * fufx(3,ng), fufy(3,ng), fufz(3,ng),
572 * mufx(3,ng), mufy(3,ng), mufz(3,ng),
573 * fumx(3,ng), fumy(3,ng), fumz(3,ng),
574 * mx(3,ng), my(3,ng), mz(3,ng),
575 * mumx(3,ng), mumy(3,ng), mumz(3,ng),
576 * flocal(3,ng,6), mlocal(3,ng,6),
577 * fbasic(3,ng,6), mbasic(3,ng,6),
578 * fdstnl(3,ng,6), mdstnl(3,ng,6),
579 * fdstnb(3,ng,6), mdstnb(3,ng,6),el(3,3,ng)
581 * denfx, denfy, denfz, denmx, denmy, denmz,
582 * refpt(3), cgmx(3), cgmy(3), cgmz(3), averef,
583 * tfufx(3), tfufy(3), tfufz(3),
584 * tmufx(3), tmufy(3), tmufz(3),
585 * tfumx(3), tfumy(3), tfumz(3),
586 * tmumx(3), tmumy(3), tmumz(3),
587 * a(6,6), c(6,6), t(3,3),smin,smax,mmax,tmax,
588 * xbar(3),rn(3),gamma(9),wi(ng),gamma_max,rndotrn,det,arm
589
590
591
592 CALL zero1(flocal,3*ng*6)
593 CALL zero1(mlocal,3*ng*6)
594 CALL zero1(fbasic,3*ng*6)
595 CALL zero1(mbasic,3*ng*6)
596 CALL zero1(fdstnl,3*ng*6)
597 CALL zero1(mdstnl,3*ng*6)
598 CALL zero1(fdstnb,3*ng*6)
599 CALL zero1(mdstnb,3*ng*6)
605 ierr = 0
606 wmin =zero
607
608 kdiag = 0
609 refpt(1) = xyz(1,ns)
610 refpt(2) = xyz(2,ns)
611 refpt(3) = xyz(3,ns)
612 DO k = 1, ng
613 DO i = 1, 3
614 tw(i,k) = frbe3(i,k)
615 rw(i,k) = frbe3(i+3,k)
616 ENDDO
617 ENDDO
618
619
620
621
622
623 IF (ng == 2.AND.irot==0) THEN
624 ierr = 322
625 GOTO 999
626 ENDIF
627
628
629
630 DO k = 1, ng
631 ielsub = ilrbe3(k)
632 IF (ielsub > 0) THEN
633 DO i = 1, 3
634 el(i,1,k) = skew(i,ielsub)
635 el(i,2,k) = skew(i+3,ielsub)
636 el(i,3,k) = skew(i+6,ielsub)
637 ENDDO
638 ENDIF
639 ENDDO
640
641
642
643 denfx = zero
644 denfy = zero
645 denfz = zero
646 averef = zero
647
648 DO 70 k = 1, ng
649 kg = inrbe3(k)
650 ielsub = ilrbe3(k)
651 IF (ielsub > 0) THEN
652
653
654
655 DO 60 i = 1, 3
656 denfx = denfx + tw(i,k)*el(i,1,k)**2
657 denfy = denfy + tw(i,k)*el(i,2,k)**2
658 denfz = denfz + tw(i,k)*el(i,3,k)**2
659 60 CONTINUE
660 ELSE
661 denfx = denfx + tw(1,k)
662 denfy = denfy + tw(2,k)
663 denfz = denfz + tw(3,k)
664 END IF
665
666 averef = averef + sqrt( (xyz(1,kg) - refpt(1))**2 +
667 * (xyz(2,kg) - refpt(2))**2 +
668 * (xyz(3,kg) - refpt(3))**2 )
669 70 CONTINUE
670
671 IF (abs(denfx) <= em20) THEN
672 ierr = 326
673 ENDIF
674
675 IF (abs(denfy) <= em20) THEN
676 ierr = 327
677 ENDIF
678
679 IF (abs(denfz) <= em20) THEN
680 ierr = 328
681 ENDIF
682 IF (ierr > 0) GOTO 999
683 averef = averef/ng
684 IF (averef == zero) averef = 1.0d0
685
686 IF (imodif==4.OR.ipen>0) THEN
687 DO k = 1, ng
688 frbe3(1,k) = frbe3(1,k)/denfx
689 frbe3(2,k) = frbe3(2,k)/denfy
690 frbe3(3,k) = frbe3(3,k)/denfz
691 frbe3(4,k) = frbe3(4,k)/denfx
692 frbe3(5,k) = frbe3(5,k)/denfy
693 frbe3(6,k) = frbe3(6,k)/denfz
694 ENDDO
695 END IF
696 IF (ipen > 0) THEN
697 xbar(1:3) = zero
698 DO k = 1, ng
699 kg = inrbe3(k)
700 wi(k) = frbe3(1,k)
701 xbar(1:3) = xbar(1:3) + wi(k)*xyz(1:3,kg)
702 ENDDO
703
704 rn(1:3) = refpt(1:3)-xbar(1:3)
705 arm = rn(1)*rn(1)+rn(2)*rn(2)+rn(3)*rn(3)
706 rndotrn = zero
707 DO k = 1, ng
708 kg = inrbe3(k)
709 rn(1:3) = xyz(1:3,kg)-xbar(1:3)
710 rndotrn =
max(rndotrn,rn(1)*rn(1)+rn(2)*rn(2)+rn(3)*rn(3))
711 END DO
712 IF (arm/rndotrn < em06) ipen =-2
713
714 IF (irot==0) THEN
715 gamma(1:9) = zero
716 DO k = 1, ng
717 kg = inrbe3(k)
718 rn(1:3) = xyz(1:3,kg)-xbar(1:3)
719 rndotrn = rn(1)*rn(1)+rn(2)*rn(2)+rn(3)*rn(3)
720
721 gamma(1) = gamma(1)+wi(k)*(rndotrn-rn(1)*rn(1))
722 gamma(2) = gamma(2)+wi(k)*( -rn(2)*rn(1))
723 gamma(3) = gamma(3)+wi(k)*( -rn(3)*rn(1))
724 gamma(4) = gamma(4)+wi(k)*( -rn(1)*rn(2))
725 gamma(5) = gamma(5)+wi(k)*(rndotrn-rn(2)*rn(2))
726 gamma(6) = gamma(6)+wi(k)*( -rn(3)*rn(2))
727 gamma(7) = gamma(7)+wi(k)*( -rn(1)*rn(3))
728 gamma(8) = gamma(8)+wi(k)*( -rn(2)*rn(3))
729 gamma(9) = gamma(9)+wi(k)*(rndotrn-rn(3)*rn(3))
730 ENDDO
731 det = (gamma(1)*(gamma(5)*gamma(9)-gamma(6)*gamma(8))-
732 * gamma(2)*(gamma(4)*gamma(9)-gamma(6)*gamma(7))+
733 * gamma(3)*(gamma(4)*gamma(8)-gamma(5)*gamma(7)))
734
735 gamma_max =
max(em20,gamma(1),gamma(5),gamma(9))
736 IF(abs(det/(gamma_max*gamma_max*gamma_max)) < em6) ierr = 400
737 END IF
738 END IF
739 IF (ierr > 0) GOTO 999
740
741
742
743
744
745 DO 40 k = 1, ng
746 kg = inrbe3(k)
747 ielsub = ilrbe3(k)
748 IF (ielsub > 0) THEN
749
750
751
752 DO 10 i = 1, 3
753 cgmx(2) = cgmx(2) + tw(i,k)*el(i,3,k)**2*xyz(2,kg)
754 cgmx(3) = cgmx(3) + tw(i,k)*el(i,2,k)**2*xyz(3,kg)
755 10 CONTINUE
756
757 DO 20 i = 1, 3
758 cgmy(3) = cgmy(3) + tw(i,k)*el(i,1,k)**2*xyz(3,kg)
759 cgmy(1) = cgmy(1) + tw(i,k)*el(i,3,k)**2*xyz(1,kg)
760 20 CONTINUE
761
762 DO 30 i = 1, 3
763 cgmz(1) = cgmz(1) + tw(i,k)*el(i,2,k)**2*xyz(1,kg)
764 cgmz(2) = cgmz(2) + tw(i,k)*el(i,1,k)**2*xyz(2,kg)
765 30 CONTINUE
766
767 ELSE
768 cgmx(2) = cgmx(2) + tw(3,k)*xyz(2,kg)
769 cgmx(3) = cgmx(3) + tw(2,k)*xyz(3,kg)
770
771 cgmy(3) = cgmy(3) + tw(1,k)*xyz(3,kg)
772 cgmy(1) = cgmy(1) + tw(3,k)*xyz(1,kg)
773
774 cgmz(1) = cgmz(1) + tw(2,k)*xyz(1,kg)
775 cgmz(2) = cgmz(2) + tw(1,k)*xyz(2,kg)
776 END IF
777 40 CONTINUE
778 cgmx(2) = cgmx(2)/denfz
779 cgmx(3) = cgmx(3)/denfy
780
781 cgmy(3) = cgmy(3)/denfx
782 cgmy(1) = cgmy(1)/denfz
783
784 cgmz(1) = cgmz(1)/denfy
785 cgmz(2) = cgmz(2)/denfx
786
787 denmx = zero
788 denmy = zero
789 denmz = zero
790
791 DO 90 k = 1, ng
792 kg = inrbe3(k)
793 ielsub = ilrbe3(k)
794
795
796
797
798
799
800 IF (ielsub > 0) THEN
801
802
803
804 DO 80 i = 1, 3
805 denmx = denmx + rw(i,k)*el(i,1,k)**2*averef**2 +
806 * tw(i,k)*( el(i,3,k)*(xyz(2,kg) - cgmx(2)) -
807 * el(i,2,k)*(xyz(3,kg) - cgmx(3))
808 * ) **2
809 denmy = denmy + rw(i,k)*el(i,2,k)**2*averef**2 +
810 * tw(i,k)*( el(i,1,k)*(xyz(3,kg) - cgmy(3)) -
811 * el(i,3,k)*(xyz(1,kg) - cgmy(1))
812 * ) **2
813 denmz = denmz + rw(i,k)*el(i,3,k)**2*averef**2 +
814 * tw(i,k)*( el(i,2,k)*(xyz(1,kg)
815 * el(i,1,k)*(xyz(2,kg) - cgmz(2))
816 * ) **2
817 80 CONTINUE
818 ELSE
819 denmx = denmx + rw(1,k)*averef**2 +
820 * tw(2,k)*(xyz(3,kg) - cgmx(3))**2 +
821 * tw(3,k)*(xyz(2,kg) - cgmx(2))**2
822 denmy = denmy + rw(2,k)*averef**2 +
823 * tw(1,k)*(xyz(3,kg) - cgmy(3))**2 +
824 * tw(3,k)*(xyz(1,kg) - cgmy(1))**2
825 denmz = denmz + rw(3,k)*averef**2 +
826 * tw(2,k)*(xyz(1,kg) - cgmz(1))**2 +
827 * tw(1,k)*(xyz(2,kg) - cgmz(2))**2
828 END IF
829 90 CONTINUE
830
831
832
833
834
835 IF (abs(denmx) <= em20) THEN
836 ierr = 329
837 ENDIF
838
839 IF (abs(denmy) <= em20) THEN
840 ierr = 330
841 ENDIF
842
843 IF (abs(denmz) <= em20) THEN
844 ierr = 331
845 ENDIF
846
847 smin =
min(abs(denmx),abs(denmy),abs(denmz))
848 smax =
max(abs(denmx),abs(denmy),abs(denmz))
849
850 IF (ierr > 0) GOTO 999
851
852 IF (irot==0 .AND.(smax/smin)>thirty) ierr = -100
853
854
855
856 CALL rbe3uf(inrbe3,ilrbe3,el,tw,xyz,refpt,
857 * fufxlc,fufylc,fufzlc,fufx,fufy,fufz,mufx,mufy,mufz,
858 * tfufx,tfufy,tfufz,tmufx,tmufy,tmufz,
859 * denfx,denfy,denfz,ng)
860
861
862
863
864
865 CALL rbe3um(inrbe3,ilrbe3,el,tw,rw,xyz,refpt,cgmx,cgmy,cgmz,
866 * fumxlc,fumylc,fumzlc,mxlc,mylc,mzlc,
867 * fumx,fumy,fumz,mx,my,mz,mumx,mumy,mumz,
868 * tfumx,tfumy,tfumz,tmumx,tmumy,tmumz,
869 * averef,denmx,denmy,denmz,ng,irot )
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891 DO 120 i = 1, 3
892 k = i + 3
893 a(i,1) = tfufx(i)
894 a(k,1) = tmufx(i)
895 a(i,2) = tfufy(i)
896 a(k,2) = tmufy(i)
897 a(i,3) = tfufz(i)
898 a(k,3) = tmufz(i)
899 a(i,4) = tfumx(i)
900 a(k,4) = tmumx(i)
901 a(i,5) = tfumy(i)
902 a(k,5) = tmumy(i)
903 a(i,6) = tfumz(i)
904 a(k,6) = tmumz(i)
905 120 CONTINUE
906
907
908
909 nsnglr = 0
911 IF (nsnglr /= 0) THEN
912 ierr = 332
913 GOTO 999
914 ENDIF
915 IF (kdiag >= 1) THEN
916 CALL wrrinf(
'C(i,1)=',c(1,1),3)
917 CALL wrrinf(
'C(i,2)=',c(1,2),3)
918 CALL wrrinf(
'C(i,3)=',c(1,3),3)
919 ENDIF
920 IF (kdiag==0.AND.ierr==0) RETURN
921
922
923
924 DO 170 j = 1, 6
925 DO 160 k = 1, ng
926 DO 150 i = 1, 3
927 flocal(i,k,j) = c(1,j)*fufxlc(i,k) + c(2,j)*fufylc(i,k) +
928 * c(3,j)*fufzlc(i,k) + c(4,j)*fumxlc(i,k) +
929 * c(5,j)*fumylc(i,k) + c(6,j)*fumzlc(i,k)
930 mlocal(i,k,j) = c(4,j)*mxlc(i,k) + c(5,j)*mylc(i,k) +
931 * c(6,j)*mzlc(i,k)
932 fbasic(i,k,j) = c(1,j)*fufx(i,k) + c(2,j)*fufy(i,k) +
933 * c(3,j)*fufz(i,k) + c(4,j)*fumx(i,k) +
934 * c(5,j)*fumy(i,k) + c(6,j)*fumz(i,k)
935 mbasic(i,k,j) = c(4,j)*mx(i,k) + c(5,j)*my(i,k) +
936 * c(6,j)*mz(i,k)
937 150 CONTINUE
938 16CONTINUE
939CONTINUE
940
941
942
943
944 DO 270 j = 1, 6
945 DO 260 k = 1, ng
946 DO 250 i = 1, 3
947 fdstnl(i,k,j) = flocal(i,k,j)
948 mdstnl(i,k,j) = mlocal(i,k,j)
949 fdstnb(i,k,j) = fbasic(i,k,j)
950 mdstnb(i,k,j) = mbasic(i,k,j)
951 250 CONTINUE
952 260 CONTINUE
953 270 CONTINUE
954
955 IF (ierr==-100) THEN
956 mmax=zero
957 DO j = 4, 6
958 DO k = 1, ng
959 DO i = 1, 3
960 IF (mmax<abs(fdstnb(i,k,j))) mmax = abs(fdstnb(i,k,j))
961 END DO
962 END DO
963 END DO
964 IF (mmax<=one) THEN
965 ierr=0
966 ELSE
967 tmax=zero
968 IF (imodif/=2) THEN
969 DO k = 1, ng
970 DO i = 1, 3
971 IF (tmax<tw(i,k)) tmax=tw(i,k)
972 ENDDO
973 ENDDO
974 ENDIF
975 wmin=tmax*em04
976 DO k = 1, ng
977 frbe3(1,k) =
max(wmin,frbe3(1,k))
978 frbe3(2,k) =
max(wmin,frbe3(2,k))
979 frbe3(3,k) =
max(wmin,frbe3(3,k))
980 ENDDO
981 ENDIF
982 END IF
983
984 999 CONTINUE
985
986
987
988 IF (kdiag >= 2) THEN
989
990
991
992 CALL wrrinf(
'TRAN_WGHTS',tw,3*ng)
993 CALL wrrinf(
'ROT_WGHTS',rw,3*ng)
994 CALL wrrinf(
'CGMX',cgmx,3)
995 CALL wrrinf(
'CGMY',cgmy,3)
996 CALL wrrinf(
'CGMZ',cgmz,3)
997 CALL wrrinf(
'DENFX',denfx,1)
998 CALL wrrinf(
'DENFY',denfy,1)
999 CALL wrrinf(
'DENFZ',denfz,1)
1000 CALL wrrinf(
'DENMX',denmx,1)
1001 CALL wrrinf(
'DENMY',denmy,1)
1002 CALL wrrinf(
'DENMZ',denmz,1)
1003 CALL wrrinf(
'AVEREF',averef,1)
1004
1005 IF (kdiag == 9.or.ierr/=0) THEN
1006 CALL wrrinf(
'FDSTNB_ULFX@REF',fdstnb(1,1,1),3*ng)
1007 CALL wrrinf(
'FDSTNB_ULFY@REF',fdstnb(1,1,2),3*ng)
1008 CALL wrrinf(
'FDSTNB_ULFZ@REF',fdstnb(1,1,3),3*ng)
1009 CALL wrrinf(
'FDSTNB_ULMX@REF',fdstnb(1,1,4),3*ng)
1010 CALL wrrinf(
'FDSTNB_ULMY@REF',fdstnb(1,1,5),3*ng)
1011 CALL wrrinf(
'FDSTNB_ULMZ@REF',fdstnb(1,1,6),3*ng)
1012 CALL wrrinf(
'MDSTNB_ULFX@REF',mdstnb(1,1,1),3*ng)
1013 CALL wrrinf(
'MDSTNB_ULFY@REF',mdstnb(1,1,2),3*ng)
1014 CALL wrrinf(
'MDSTNB_ULFZ@REF',mdstnb(1,1,3),3*ng)
1015 CALL wrrinf(
'MDSTNB_ULMX@REF',mdstnb(1,1,4),3*ng)
1016 CALL wrrinf(
'MDSTNB_ULMY@REF',mdstnb(1,1,5),3*ng)
1017 CALL wrrinf(
'MDSTNB_ULMZ@REF',mdstnb(1,1,6),3*ng)
1018 ENDIF
1019 IF (kdiag >= 30) THEN
1020 CALL wrrinf(
'FDSTNL_ULFX@REF',fdstnl(1,1,1),3*ng)
1021 CALL wrrinf(
'FDSTNL_ULFY@REF',fdstnl(1,1,2),3*ng)
1022 CALL wrrinf(
'FDSTNL_ULFZ@REF',fdstnl(1,1,3),3*ng)
1023 CALL wrrinf(
'FDSTNL_ULMX@REF',fdstnl(1,1,4),3*ng)
1024 CALL wrrinf(
'FDSTNL_ULMY@REF',fdstnl(1,1,5),3*ng)
1025 CALL wrrinf('fdstnl_ulmz@ref
',FDSTNL(1,1,6),3*NG)
1026 CALL WRRINF('mdstnl_ulfx@ref',mdstnl(1,1,1),3*ng)
1027 CALL wrrinf(
'MDSTNL_ULFY@REF',mdstnl(1,1,2),3*ng)
1028 CALL wrrinf(
'MDSTNL_ULFZ@REF',mdstnl(1,1,3),3*ng)
1029 CALL wrrinf(
'MDSTNL_ULMX@REF',mdstnl(1,1,4),3*ng)
1030 CALL wrrinf(
'MDSTNL_ULMY@REF',mdstnl(1,1,5),3*ng)
1031 CALL wrrinf(
'MDSTNL_ULMZ@REF',mdstnl(1,1,6),3*ng)
1032
1033
1034 ENDIF
1035 ENDIF
1036
1037 RETURN
subroutine invert(matrix, inverse, n, errorflag)
subroutine wrrinf(title, r, n)
subroutine rbe3um(inrbe3, ilrbe3, el, tw, rw, xyz, refpt, cgmx, cgmy, cgmz, fumxlc, fumylc, fumzlc, mxlc, mylc, mzlc, fumx, fumy, fumz, mx, my, mz, mumx, mumy, mumz, tfumx, tfumy, tfumz, tmumx, tmumy, tmumz, averef, denmx, denmy, denmz, ng, irot)
subroutine rbe3uf(inrbe3, ilrbe3, el, tw, xyz, refpt, fufxlc, fufylc, fufzlc, fufx, fufy, fufz, mufx, mufy, mufz, tfufx, tfufy, tfufz, tmufx, tmufy, tmufz, denfx, denfy, denfz, ng)