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