Диссертация (531291), страница 36
Текст из файла (страница 36)
IQT372')103 FORMAT(1X,5X,'Nelem = ',I3,5X,' Velem =',E12.5)STOPENDSUBROUTINE ISOQ72(X,Y,Z,VOL,B,C,D,ERG,A)COMMON/CONTR/NP,NE,NB,NDF,NCN,NSZF,CMY,XMAS,YMAS,RO,ZMAS,EMATDIMENSION AA(3,24),L(3),M(3),ERG(72,72),A(72,72)DIMENSION S(3,3),B(1),C(1),D(1),ANI(24),AB(3,72)272NSTIF=72C Определение матрицы AB=N (Ф.формы в точке интегрирования)AB = 0.CALL NFOR72(X,Y,Z,ANI)DO 30 I=1,NCNII=(I-1)*3+1AB(1,II)=ANI(I)AB(2,II+1)=ANI(I)30AB(3,II+2)=ANI(I)c Определение A=Nт*N в (.) интегрированияDO 40 I=1,NSTIFDO 40 J=1,NSTIFA(I,J)=0.DO 40 K=1,340A(I,J)=A(I,J)+AB(K,I)*AB(K,J)CОпрелеление производных функций формы по естественным координатамCALL GRAD72(X,Y,Z,AA)Определение матрицы [J]DO 45 I=1,3DO 45 J=1,345S(I,J)=0.DO 46 J=1,NCNS(1,1)=S(1,1)+AA(1,J)*B(J)S(1,2)=S(1,2)+AA(1,J)*C(J)S(1,3)=S(1,3)+AA(1,J)*D(J)S(2,1)=S(2,1)+AA(2,J)*B(J)S(2,2)=S(2,2)+AA(2,J)*C(J)S(2,3)=S(2,3)+AA(2,J)*D(J)S(3,1)=S(3,1)+AA(3,J)*B(J)S(3,2)=S(3,2)+AA(3,J)*C(J)46S(3,3)=S(3,3)+AA(3,J)*D(J)CALL MINV(S,3,VOL,L,M,.TRUE.)c обратный переход к глобальным координатамAB = 0.DO 60 I=1,NCNII=(I-1)*3+1AB(1,II )=S(1,1)*AA(1,I)+S(1,2)*AA(2,IAB(2,II+1)=S(2,1)*AA(1,I)+S(2,2)*AA(2,I60AB(3,II+2)=S(3,1)*AA(1,I)+S(3,2)*AA(2,IC Определение матрицы градиентов ERG=BERG = 0.DO 70 I=1,3DO 70 J=1,NSTIF70ERG(I,J)=AB(I,J)J=NCN*3-2DO 80 I=1,J,3ERG(4,I )=ERG(2,I+1)ERG(4,I+1)=ERG(1,I )ERG(5,I )=ERG(3,I+2)ERG(5,I+2)=ERG(1,I)ERG(6,I+1)=ERG(3,I+2)80ERG(6,I+2)=ERG(2,I+1)RETURNENDcC)+S(1,3)*AA(3,I)+S(2,3)*AA(3,I)+S(3,3)*AA(3,I)))SUBROUTINE GRAD72(X,Y,Z,AA)DIMENSION AA(3,24)Опрелеление производных функций формы по естественным коорди-натамE =1.X1=E-XX2=E+XY1=E-Y273CCY2=E+YZ1=E-ZZ2=E+ZX22=E-X*XZ22=E-Z*ZОпределение dN/dxF1=(10.+18.*X-27.*X*X-9.*Z*Z)/64.F2=(-10.+18.*X+27.*X*X+9.*Z*Z)/64.AA(1, 1)= Y1 *Z1 *F1AA(1, 4)= Y1 *Z1 *F2AA(1, 5)= Y2 *Z1 *F2AA(1,08)= Y2 *Z1 *F1AA(1,17)= Y1 *Z2 *F1AA(1,20)= Y1 *Z2 *F2AA(1,21)= Y2 *Z2 *F2AA(1,24)= Y2 *Z2 *F1F1=9.*(1.-3.*Z)*Z22/64.F2=9.*(1.+3.*Z)*Z22/64.AA(1,10)= Y1*F1AA(1,14)= Y1*F2AA(1,11)= Y2*F1AA(1,15)= Y2*F2AA(1, 9)= -AA(1,10)AA(1,13)= -AA(1,14)AA(1,12)= -AA(1,11)AA(1,16)= -AA(1,15)F1=9.*(-2.*X-3.+9.*X*X)/64.F2=9.*(-2.*X+3.-9.*X*X)/64.AA(1, 2)= Y1 * Z1 *F1AA(1, 3)= Y1 * Z1 *F2AA(1, 7)= Y2 * Z1 *F1AA(1, 6)= Y2 * Z1 *F2AA(1,18)= Y1 * Z2 *F1AA(1,19)= Y1 * Z2 *F2AA(1,23)= Y2 * Z2 *F1AA(1,22)= Y2 * Z2 *F2Определение dN/dyF1=(-10.+9.*(X*X+Z*Z))/64.F2=F1*Z2F1=F1*Z1AA(2, 5)= F1 * X2AA(2, 8)= F1 * X1AA(2,21)= F2 * X2AA(2,24)= F2 * X1AA(2, 4)= -AA(2, 5)AA(2, 1)= -AA(2, 8)AA(2,20)= -AA(2,21)AA(2,17)= -AA(2,24)F1=Z22*(1.-3.*Z)*9./64.F2=Z22*(1.+3.*Z)*9./64.AA(2,12)= F1 * X1AA(2,16)= F2 * X1AA(2,11)= F1 * X2AA(2,15)= F2 * X2AA(2, 9)= -AA(2,12)AA(2,13)= -AA(2,16)AA(2,10)= -AA(2,11)AA(2,14)= -AA(2,15)F1=9.*X22*(1.-3.*X)/64.F2=9.*X22*(1.+3.*X)/64.AA(2, 6)= F2 * Z1AA(2, 7)= F1 * Z1AA(2,22)= F2 * Z2AA(2,23)= F1 * Z2274CAA(2, 3)= -AA(2, 6)AA(2, 2)= -AA(2, 7)AA(2,19)= -AA(2,22)AA(2,18)= -AA(2,23)Определение dN/dzF1=(10.+18.*Z-9.*X*X-27.*Z*Z)/64.F2=(-10.+18.*Z+9.*X*X+27.*Z*Z)/64.AA(3, 1)= X1 *Y1 * F1AA(3, 4)= X2 *Y1 * F1AA(3, 5)= X2 *Y2 * F1AA(3, 8)= X1 *Y2 * F1AA(3,17)= X1 *Y1 * F2AA(3,20)= X2 *Y1 * F2AA(3,21)= X2 *Y2 * F2AA(3,24)= X1 *Y2 * F2F1=9.*(-3.-2.*Z+9.*Z*Z)/64.F2=9.*(3.-2.*Z-9.*Z*Z)/64.AA(3, 9)= X1 * Y1 * F1AA(3,13)= X1 * Y1 * F2AA(3,12)= X1 * Y2 * F1AA(3,16)= X1 * Y2 * F2AA(3,10)= X2 * Y1 * F1AA(3,14)= X2 * Y1 * F2AA(3,11)= X2 * Y2 * F1AA(3,15)= X2 * Y2 * F2F1=9.*X22*(1.-3.*X)/64.F2=9.*X22*(1.+3.*X)/64.AA(3,18)= Y1 * F1AA(3,19)= Y1 * F2AA(3,22)= Y2 * F2AA(3,23)= Y2 * F1AA(3, 2)= -AA(3,18)AA(3, 3)= -AA(3,19)AA(3, 6)= -AA(3,22)AA(3, 7)= -AA(3,23)RETURNENDSUBROUTINE NFOR72(X,Y,Z,AA)DIMENSION AA(1)E=1.X1=E-XX2=E+XY1=E-YY2=E+YZ1=E-ZZ2=E+ZA =(9.*(X*X+Z*Z)-10.)/64.F1=A*Z1F2=A*Z2AA( 1) = X1 * Y1 * F1AA( 4) = X2 * Y1 * F1AA( 5) = X2 * Y2 * F1AA( 8) = X1 * Y2 * F1AA( 17) = X1 * Y1 * F2AA( 20) = X2 * Y1 * F2AA( 21) = X2 * Y2 * F2AA( 24) = X1 * Y2 * F2F1=9.*(E-Z*Z)/64.F2=F1*(E+3.*Z)F1=F1*(E-3.*Z)AA( 9) = F1 * X1 * Y1AA( 13) = F2 * X1 * Y1AA( 12) = F1 * X1 * Y2275AA(AA(AA(AA(AA(16) = F2 * X1 * Y210) = F1 * X2 * Y114) = F2 * X2 * Y111) = F1 * X2 * Y215) = F2 * X2 * Y2F1=9.*(E-X*X)/64.F2=F1*(E+3.*X)F1=F1*(E-3.*X)AA( 2) = F1 * Z1 * Y1AA( 3) = F2 * Z1 * Y1AA( 6) = F2 * Z1 * Y2AA( 7) = F1 * Z1 * Y2AA( 18) = F1 * Z2 * Y1AA( 19) = F2 * Z2 * Y1AA( 22) = F2 * Z2 * Y2AA( 23) = F1 * Z2 * Y2RETURNENDSUBROUTINE DDTED(C,D,E,NCD,NE,NCDC,NCDD,NCDE,S)C Определение [C]=[D]T*[E]*[D]C Входные параметры:C D -матрица, размерностью (NE x NCD)C E -матрица, размерностью (NE x NE]C S –рабочий массив, размерностью не менее (NCD x NE)C!!! NCDC - ЧИСЛО СТРОК МАТРИЦЫ С В ВЫЗЫВАЮЩЕЙ ПРОГРАММЕC!!! NCDD - ЧИСЛО СТРОК МАТРИЦЫ D В ВЫЗЫВАЮЩЕЙ ПРОГРАММЕC!!! NCDE - ЧИСЛО СТРОК МАТРИЦЫ E В ВЫЗЫВАЮЩЕЙ ПРОГРАММЕC Выходные параметры:C C -матрица, размерностью (NCD x NCD)double precision rsDIMENSION C(NCDC,1),D(NCDD,1),E(NCDE,1),S(NCD,1)DO 2 I=1,NCDDO 2 J=1,NERS=0.DO 1 L=1,NE1 RS=D(L,I)*E(L,J)+RS2 S(I,J)=RSDO 4 I=1,NCDDO 4 J=1,NCDRS=0.DO 3 L=1,NE3 RS=S(I,L)*D(L,J)+RS4 C(I,J)=RSRETURNENDCcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccCminv subsystem (обращение матриц)CSUBROUTINE MINV(A,N,D,L,M,FLAG)DIMENSION A(1),L(1),M(1)LOGICAL FLAGCALL MINV1(A,N,D,L,M,FLAG)IF(D.NE.0.) CALL MINV2(A,N,L,M)RETURNENDSUBROUTINE MINV1(A,N,D,L,M,FLAG)C Часть программы minv03, определяющая det|A| (вых-ой параметр D)C если FLAG = .TRUE.
и не определяющая если FLAG = .FALSE.LOGICAL FLAG27620303540454855657580DIMENSION A(1),L(1),M(1)D=1.NK=-NDO 80 K=1,NNK=NK+NL(K)=KM(K)=KKK=NK+KBIGA=A(KK)DO 20 J=K,NIZ=N*(J-1)DO 20 I=K,NIJ=IZ+IIF(ABS(BIGA)-ABS(A(IJ)).GE.0.) GO TO 20BIGA=A(IJ)L(K)=IM(K)=JCONTINUEJ=L(K)IF(J-K.LE.0) GO TO 35KI=K-NDO 30 I=1,NKI=KI+NHOLD=-A(KI)JI=KI-K+JA(KI)=A(JI)A(JI)=HOLDI=M(K)IF(I-K.LE.0) GO TO 45JP=N*(I-1)DO 40 J=1,NJK=NK+JJI=JP+JHOLD=-A(JK)A(JK)=A(JI)A(JI)=HOLDIF(ABS(BIGA).GT.0.) GO TO 48D=0.RETURNDO 55 I=1,NIF(I-K.EQ.0) GO TO 55IK=NK+IA(IK)=A(IK)/(-BIGA)CONTINUEDO 65 I=1,NIK=NK+IHOLD=A(IK)IJ=I-NDO 65 J=1,NIJ=IJ+NIF(I-K.EQ.0) GO TO 65IF(J-K.EQ.0) GO TO 65KJ=IJ-I+KA(IJ)=HOLD*A(KJ)+A(IJ)CONTINUEKJ=K-NDO 75 J=1,NKJ=KJ+NIF(J-K.EQ.0) GO TO 75A(KJ)=A(KJ)/BIGACONTINUEIF(FLAG) D=D*BIGAA(KK)=1./BIGACONTINUE277RETURNEND100110120130150SUBROUTINE MINV2(A,N,L,M)DIMENSION A(1),L(1),M(1)K=NK=K-1IF(K.LE.0) GO TO 150I=L(K)IF(I-K.LE.0) GO TO 120JQ=N*(K-1)JR=N*(I-1)DO 110 J=1,NJK=JQ+JHOLD=A(JK)JI=JR+JA(JK)=-A(JI)A(JI)=HOLDJ=M(K)IF(J-K.LE.0) GO TO 100KI=K-NDO 130 I=1,NKI=KI+NHOLD=A(KI)JI=KI-K+JA(KI)=-A(JI)A(JI)=HOLDGO TO 100CONTINUERETURNEND278.