C********************************************************************* C * C MAIN PROGRAM OF POTENTIAL FLOW TO SOLVE LAPLACE EQUATIONS * C BY DUAL INTEGRAL EQUATIONS WITH U(s,x),T(s,x),L(s,x),M(s,x) * C KERNEL FUNCTIONS THROUGTH BEM BY USING CONSTANT ELEMENTS. * C * C This program can be executed in VAX and CRAY system * C * C modified by J. T. Chen Aug.23,1991 * C using invariant integral and complete theory of * C potential theory for normal and tangent derivative* C of single and double potential * C input for001.dat(neumann) output for016.dat(boundary data)* C input for002.dat(dirichlet) output for077.dat(interior pote)* C input for003.dat(BC) output for078.dat(interior flux)* C input for015.dat(mesh geom) * C input for080.dat(interior pt. geom) * C********************************************************************* C IMPLICIT REAL*8 (A-H,O-Z) COMMON/AA/ NELM,NNODE COMMON/NORMAL/ XNORMAL,YNORMAL,THETA ! COMMON/BAR/ XBAR,YBAR,XROT,YROT DIMENSION INC(2,1000),ELENGTH(1000),THETA(1000),COLY(1000) DIMENSION X(1000),Y(1000),NINDEX(1000),COLX(1000) DIMENSION U(1000,1000),T(1000,1000),VL(1000,1000),VM(1000,1000) DIMENSION XNORMAL(1000),YNORMAL(1000),VKNOWN(1000) DIMENSION BCU(1000),BCT(1000),B(1000),L(1000) DIMENSION TMOD(1000,1000),UMOD(1000,1000),B2(1000) * ,VLMOD(1000,1000),VMMOD(1000,1000),BCULM(1000),BCTLM(1000) DIMENSION P1(1000),P2X(1000),P2Y(1000),M1(1000) OPEN(UNIT=11,file='f01.dat',status='OLD') OPEN(UNIT=12,file='f02.dat',status='OLD') OPEN(UNIT=13,file='f03.dat',status='OLD') OPEN(UNIT=15,file='f15.dat',status='OLD') C OPEN(UNIT=5,file='f05.dat') C OPEN(UNIT=6,file='f06.dat') OPEN(UNIT=80,file='f80.dat',status='OLD') OPEN(UNIT=16,file='f16.dat') OPEN(UNIT=77,file='f77.dat') OPEN(UNIT=78,file='F78.dat') PRINT *,'INPUT NELM' READ(5,*)NELM PRINT *,'INPUT NINTER' READ(5,*)NINTER DO 901 I=1,NELM BCU(I)=9999. BCT(I)=9999. 901 CONTINUE CALL BC(BCU,BCT,VKNOWN) CALL GEOM(X,Y,INC,NINDEX,NELM,NNODE) C DO I=1,NNODE C WRITE(91,101) (I,X(I),Y(I),Z(I)) C 101 FORMAT(1X,'NODE NO.',I3,5X,'X COORD',1X,f13.5,5X,'Y COORD',1X, C * f13.5, 5X,'Z COORD',1X,f13.5) C ENDDO DO 902 I=1,NELM J=NINDEX(I) COLX(J)=0.5*(X(INC(1,J))+X(INC(2,J))) !±`¼Æ¤¸¯À§¤®y¼ÐÀ COLY(J)=0.5*(Y(INC(1,J))+Y(INC(2,J))) DELTX=X(INC(2,J))-X(INC(1,J)) DELTY=Y(INC(2,J))-Y(INC(1,J)) ELENGTH(J)=SQRT(DELTX**2+DELTY**2) IF(ABS(DELTY).LE.0.00000001) GO TO 121 ARCCOT=90.0-ATAN(DELTX/DELTY)*180./3.141592654 THETA(J)=ARCCOT IF(DELTY.LT.0.0) THETA(J)=ARCCOT+180. 121 IF(ABS(DELTY).LE.0.00000001.AND.DELTX.GE.0.0) THETA(J)=0. IF(ABS(DELTY).LE.0.00000001.AND.DELTX.LE.0.0) THETA(J)=180. THETT=THETA(J)*3.141592654/180. SS=THETT-0.5*3.141592654 XNORMAL(J)=COS(SS) YNORMAL(J)=SIN(SS) C WRITE(94,106) J,XNORMAL(J),YNORMAL(J) 106 FORMAT(1X,'NEL=',I3,1X,'XNORM COMP',F10.5,1X,'YNORM COMP',F10.5) C 200 WRITE(91,102) (J,INC(1,J),INC(2,J)) C WRITE(91,103) (J,COLX(J),COLY(J),ELENGTH(J),THETA(J)) 102 FORMAT(1X,'ELEM NO.= ',I3,5X,'INC(1,J)=',1X,I3,5X,'INC(2,J)=' * ,1X,I3) 103 FORMAT(1X,'COLO NO.= ',I3,5X,'COL X COORD=',1X,F13.5,5X, * 'COL Y COORD=' * ,1X,F13.5,5X,/,'ELEM LENGTH =',1X,F13.5,5X,'THETA =',1X,F7.2) 902 CONTINUE C********************************************************* DO 111 ICOL=1,NELM DO 222 NEL=1,NELM XBAR=COLX(ICOL)-COLX(NEL) YBAR=COLY(ICOL)-COLY(NEL) THET=THETA(NEL)*3.141592654/180. XROT=XBAR*COS(THET)+YBAR*SIN(THET) YROT=-XBAR*SIN(THET)+YBAR*COS(THET) RL=-0.5*ELENGTH(NEL) RH=0.5*ELENGTH(NEL) C WRITE(91,104) ICOL,NEL,XROT,YROT,RL,RH CALL UKERNEL(NEL,XROT,YROT,RL,RH,ELENGTH,VALUEU) CALL TKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUET) CALL LKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEL) CALL MKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEM) U(ICOL,NEL)=VALUEU T(ICOL,NEL)=VALUET VL(ICOL,NEL)=VALUEL VM(ICOL,NEL)=VALUEM 104 FORMAT(1X,'ICOL=',1X,I3,1X,'NEL=',1X,I3,1X,'XROT=',1X,F10.5,1X, * 'YROT=',1X,F10.5,1X,'RL=',F5.2,1X,'RH=',F5.2) 222 CONTINUE 111 CONTINUE WRITE(16,1000) DO IP=1,NELM WRITE(16,333) (U(IP,IQ),IQ=1,NELM) ENDDO WRITE(16,2000) DO IP=1,NELM WRITE(16,333) (T(IP,IQ),IQ=1,NELM) ENDDO WRITE(16,3000) DO IP=1,NELM WRITE(16,333) (VL(IP,IQ),IQ=1,NELM) ENDDO WRITE(16,4000) DO IP=1,NELM WRITE(16,333) (VM(IP,IQ),IQ=1,NELM) ENDDO 333 FORMAT(10X,16F10.5) 1000 FORMAT(10X,//,' U(I,J) OF U(s,x) INTEGRATION ') 2000 FORMAT(10X,//,' T(I,J) OF T(s,x) INTEGRATION ') 3000 FORMAT(10X,//,' L(I,J) OF L(s,x) INTEGRATION ') 4000 FORMAT(10X,//,' M(I,J) OF M(s,x) INTEGRATION ') C************************************************************** C* FOR DEGENERATE BOUNDARY TREATMENT * C* * C************************************************************* DO 77 I=1,NELM DO 77 J=I+1,NELM IF(COLX(I).NE.COLX(J).OR.COLY(I).NE.COLY(J)) GO TO 77 DO 903 JJ=1,NELM TEMP1=VM(I,JJ) TEMP2=T(I,JJ) TEMP3=VL(I,JJ) TEMP4=U(I,JJ) T(I,JJ)=TEMP1 VM(I,JJ)=TEMP2 U(I,JJ)=TEMP3 VL(I,JJ)=TEMP4 903 continue 77 CONTINUE C WRITE(6,1000) C DO IP=1,NELM C WRITE(6,333) (U(IP,IQ),IQ=1,NELM) C ENDDO C WRITE(6,2000) C DO IP=1,NELM C WRITE(6,333) (T(IP,IQ),IQ=1,NELM) C ENDDO C WRITE(6,3000) C DO IP=1,NELM C WRITE(6,333) (VL(IP,IQ),IQ=1,NELM) C ENDDO C WRITE(6,4000) C DO IP=1,NELM C WRITE(6,333) (VM(IP,IQ),IQ=1,NELM) C ENDDO C************************************************************** C MATRIX OPERATION * C************************************************************** ibcu=0 nofu=0 nofk=0 DO 904 I=1,NELM IF(BCU(I).EQ.9999.) IBCU=IBCU+1 904 continue M=IBCU !UNKNOWN OF BCU N=NELM-IBCU !KNOWN OF BCU DO 905 J=1,NELM IF(BCU(J).EQ.9999.) THEN NOFU=NOFU+1 DO 906 I=1,NELM TMOD(I,NOFU)=T(I,J) VMMOD(I,NOFU)=VM(I,J) UMOD(I,NOFU)=U(I,J) VLMOD(I,NOFU)=VL(I,J) 906 continue ELSE NOFK=NOFK+1 lH1=NOFK+M DO 907 I=1,NELM TMOD(I,lH1)=T(I,J) VMMOD(I,lH1)=VM(I,J) UMOD(I,lH1)=U(I,J) VLMOD(I,lH1)=VL(I,J) 907 continue ENDIF 905 continue c WRITE(6,2000) c DO IP=1,NELM c WRITE(6,333) (TMOD(IP,IQ),IQ=1,NELM) c ENDDO C*********PARTITIONING************ C WRITE(6,1021) M,N,NELM C1021 FORMAT(1X,' M AND N AND NELM = ',3I5) IF(M.EQ.0) GO TO 7 DO 910 I=1,M DO 910 J=1,M T(I,J)=TMOD(I,J) U(I,J)=UMOD(I,J) VM(I,J)=VMMOD(I,J) VL(I,J)=VLMOD(I,J) 910 continue C---------------------------------- IF(M.EQ.NELM) GO TO 8 DO 911 I=M+1,NELM DO 911 J=1,M T(I,J)=TMOD(I,J) U(I,J)=UMOD(I,J) VM(I,J)=VMMOD(I,J) VL(I,J)=VLMOD(I,J) 911 continue C--------------------------------- DO 912 I=1,M DO 912 J=M+1,NELM T(I,J)=-UMOD(I,J) U(I,J)=-TMOD(I,J) VM(I,J)=-VLMOD(I,J) VL(I,J)=-VMMOD(I,J) 912 continue C----------------------------------- 7 DO 913 I=M+1,NELM DO 913 J=M+1,NELM T(I,J)=-UMOD(I,J) U(I,J)=-TMOD(I,J) VM(I,J)=-VLMOD(I,J) VL(I,J)=-VMMOD(I,J) 913 continue C------------------------------- 8 CONTINUE DO 314 IB=1,NELM B(IB)=0. B2(IB)=0. 314 continue DO 914 IB=1,NELM DO 914 IL=1,NELM B(IB)=B(IB)+U(IB,IL)*VKNOWN(IL) B2(IB)=B2(IB)+VL(IB,IL)*VKNOWN(IL) 914 continue C WRITE(6,5000) C5000 FORMAT(1X,'U(I,J) MATRIX') c DO IP=1,NELM c WRITE(6,333) (T(IP,IQ),IQ=1,NELM) c ENDDO CALL LUPPDC(VM,L,NELM) CALL LUPPSB(VM,B2,L,NELM) CALL LUPPDC(T,L,NELM) CALL LUPPSB(T,B,L,NELM) C*********** ANS QUEUE FOR U T KERNEL ******************8 ik=0 DO 915 I=1,NELM IF(BCU(I).EQ.9999.) THEN IK=IK+1 BCULM(I)=B2(IK) BCU(I)=B(IK) ELSE BCULM(I)=BCU(I) ENDIF 915 continue it=0 DO 916 J=1,NELM IF(BCT(J).EQ.9999.) THEN IT=IT+1 BCT(J)=B(IT+M) BCTLM(J)=B2(IT+M) ELSE BCTLM(J)=BCT(J) ENDIF 916 continue WRITE(16,1100) (I,BCU(I),I=1,NELM) WRITE(16,1200) (I,BCT(I),I=1,NELM) 1100 FORMAT(1X,'U T KERNEL ANS ','BCU(',I3,')=',F10.3) 1200 FORMAT(1X,'U T KERNEL ANS ','BCT(',I3,')=',F10.3) WRITE(16,2100) (I,BCULM(I),I=1,NELM) WRITE(16,2200) (I,BCTLM(I),I=1,NELM) 2100 FORMAT(1X,'L M KERNEL ANS ','BCU(',I3,')=',F10.3) 2200 FORMAT(1X,'L M KERNEL ANS ','BCT(',I3,')=',F10.3) C*********** INTERIOR POINT ********************************* DO 7001 K=1,NINTER KB=NELM+K READ(80,100) M1(K),M2,M3,M4,COLX(KB),COLY(KB),R 100 FORMAT(4I10,3E13.5) 501 FORMAT(80A) 7001 CONTINUE C************************************** C********** FOR X GRADIENT ************ DO 920 ICOL=NELM+1,NELM+NINTER XNORMAL(ICOL)=1.0 YNORMAL(ICOL)=0.0 DO 920 NEL=1,NELM XBAR=COLX(ICOL)-COLX(NEL) YBAR=COLY(ICOL)-COLY(NEL) THET=THETA(NEL)*3.141592654/180. XROT=XBAR*COS(THET)+YBAR*SIN(THET) YROT=-XBAR*SIN(THET)+YBAR*COS(THET) RL=-0.5*ELENGTH(NEL) RH=0.5*ELENGTH(NEL) CALL UKERNEL(NEL,XROT,YROT,RL,RH,ELENGTH,VALUEU) CALL TKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUET) CALL LKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEL) CALL MKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEM) U(ICOL,NEL)=VALUEU T(ICOL,NEL)=VALUET VL(ICOL,NEL)=VALUEL VM(ICOL,NEL)=VALUEM C write(1,*)u(icol,nel),t(icol,nel),vl(icol,nel),vm(icol,nel) 920 CONTINUE DO 930 ICOL=1,NINTER P1(ICOL)=0.0 P2X(ICOL)=0.0 DO 930 NEL=1,NELM P1(ICOL)=P1(ICOL)+T(ICOL+NELM,NEL)*BCU(NEL) * -U(ICOL+NELM,NEL)*BCT(NEL) P2X(ICOL)=P2X(ICOL)+VM(ICOL+NELM,NEL)*BCULM(NEL) * -VL(ICOL+NELM,NEL)*BCTLM(NEL) c write(*,*)p1(icol),p2x(icol) 930 continue C*********************************************** C************** FOR Y GRADIENT ******* DO 940 ICOL=NELM+1,NELM+NINTER XNORMAL(ICOL)=0.0 YNORMAL(ICOL)=1.0 DO 940 NEL=1,NELM XBAR=COLX(ICOL)-COLX(NEL) YBAR=COLY(ICOL)-COLY(NEL) THET=THETA(NEL)*3.141592654/180. XROT=XBAR*COS(THET)+YBAR*SIN(THET) YROT=-XBAR*SIN(THET)+YBAR*COS(THET) RL=-0.5*ELENGTH(NEL) RH=0.5*ELENGTH(NEL) CALL LKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEL) CALL MKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEM) VL(ICOL,NEL)=VALUEL VM(ICOL,NEL)=VALUEM 940 continue DO 950 ICOL=1,NINTER P2Y(ICOL)=0.0 DO 950 NEL=1,NELM P2Y(ICOL)=P2Y(ICOL)+VM(ICOL+NELM,NEL)*BCULM(NEL) * -VL(ICOL+NELM,NEL)*BCTLM(NEL) 950 continue AIN=2. DO 951 MOD=1,NINTER P1(MOD)=P1(MOD)/(2.0*3.141592654) P2X(MOD)=P2X(MOD)/(2.0*3.141592654) P2Y(MOD)=P2Y(MOD)/(2.0*3.141592654) WRITE(77,130)M1(MOD),P1(MOD) WRITE(78,140)M1(MOD),P2X(MOD),P2Y(MOD) AIN=AIN+1. 951 continue 130 FORMAT(I10,E13.5) 140 FORMAT(I10,2E13.5) STOP END C************************************************************** C UKERNEL SUBROUTINE CAN CALCULATE THE U(s,x) INTEGRATION * C INCLUDING SINGULAR AND NONSINGULAR VALUE * C * C************************************************************** SUBROUTINE UKERNEL(NEL,XROT,YROT,RL,RH,ELENGTH,VALUEU) c IMPLICIT REAL*8 (A-H,O-Z) EXTERNAL FUNOSI1,FUNOSI2 DIMENSION ELENGTH(1000) save c write(6,*) Icol, nel, rh, rl IF(ABS(XROT).LE.0.00000001.AND.ABS(YROT).LE.0.00000001) THEN VALUEU=RH*LOG(ABS(RH))-RH-(RL*LOG(ABS(RL))-RL) c write(6,*) icol, nel, rh, rl, yrot CHECK=FUNOSI1(RH,YROT)-FUNOSI1(RL,YROT) ELSE c CALL GRULE(50,XW,W) c CALL ARANGE(50,XW,W) c UVALUE=0.0 C-----------IF YOU CHECK THE INTEGRAL USING GAUSSIAN QUATURE -------- C DO 555 I=1,100 C UVALUE=UVALUE+FUNOSI2(ELENGTH(NEL),XW(I),XROT,YROT) C * *W(I)*0.5*ELENGTH(NEL) C555 CONTINUE C-------------------------------------------------------------------- C CHECK=UVALUE RTOP=0.5*ELENGTH(NEL)-XROT RLOW=-0.5*ELENGTH(NEL)-XROT VALUEU=FUNOSI1(RTOP,YROT)-FUNOSI1(RLOW,YROT) ENDIF c IF(CHECK.EQ.0.0) RETURN TESTVAL=(CHECK-VALUEU)/CHECK IF(ABS(TESTVAL).GE.0.001) THEN c WRITE(6,190) 190 FORMAT(1X,'ERROR OF U(S,X) INTEGRATION') ELSE ENDIF RETURN END C**************************************************** FUNCTION FUNOSI1(V,Y0) c IMPLICIT REAL*8 (A-H,O-Z) save c write(6,*) v, y0 IF(Y0.NE.0.0) THEN FUNOSI1=V*LOG(SQRT(V*V+Y0*Y0))-V+Y0*ATAN(V/Y0) ELSE c write(6,*) v FUNOSI1=V*LOG(ABS(V))-V ENDIF RETURN END C**************************************************** C******** FOR GAUSS INTEGRATION USE ************* FUNCTION FUNOSI2(EL,T,XROT,YROT) c IMPLICIT REAL*8 (A-H,O-Z) FUNOSI2=LOG(SQRT(YROT**2+((0.5*EL*T-XROT)**2))) RETURN END CC************************************************************** C TKERNEL SUBROUTINE CAN CALCULATE THE T(s,x) INTEGRATION * C INCLUDING SINGULAR AND NONSINGULAR VALUE * C * C************************************************************** SUBROUTINE TKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUET) c IMPLICIT REAL*8 (A-H,O-Z) COMMON /NORMAL/ XNORMAL,YNORMAL,THETA ! COMMON/BAR/ XBAR,YBAR COMMON/AA/ NELM,NNODE EXTERNAL T1 DIMENSION ELENGTH(1000) * ,XNORMAL(1000),YNORMAL(1000) * ,THETA(1000) IF(ABS(YROT).LE.0.00000001) THEN VALUET=-3.141592654 IF(ICOL.GT.NELM) VALUET=3.141592654 IF(ABS(XROT).GE.0.00000001) VALUET=0.0 ELSE RTOP=0.5*ELENGTH(NEL)-XROT RLOW=-0.5*ELENGTH(NEL)-XROT XN=XNORMAL(NEL) YN=YNORMAL(NEL) CDAR=THETA(NEL)*3.141592654/180. VALUET=T1(RTOP,YROT)-T1( * RLOW,YROT) ENDIF RETURN END C**************************************************** FUNCTION T1(V,Y0) c IMPLICIT REAL*8 (A-H,O-Z) T1=ATAN(V/Y0) RETURN END CC************************************************************* C LKERNEL SUBROUTINE CAN CALCULATE THE L(s,x) INTEGRATION * C INCLUDING SINGULAR AND NONSINGULAR VALUE * C BY CLOSED FORM SOLUTION * C************************************************************** SUBROUTINE LKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEL) c IMPLICIT REAL*8 (A-H,O-Z) COMMON /NORMAL/ XNORMAL,YNORMAL,THETA ! COMMON/BAR/ XBAR,YBAR COMMON/AA/ NELM,NNODE EXTERNAL VL1 DIMENSION ELENGTH(1000) * ,XNORMAL(1000),YNORMAL(1000) * ,THETA(1000) RTOP=0.5*ELENGTH(NEL)-XROT RLOW=-0.5*ELENGTH(NEL)-XROT XN=XNORMAL(ICOL) XS=XNORMAL(NEL) YN=YNORMAL(ICOL) YS=YNORMAL(NEL) CDAR=THETA(NEL)*3.141592654/180. DOT=XN*XS+YN*YS CROSS=-1.*XN*YS+YN*XS IF(ABS(YROT).LE.0.000001) THEN VALUEL=3.141592654 IF(ABS(XROT).LE.0.00000001.AND.ICOL.NE.NEL) VALUEL=-3.141592654 TAGENT1=-0.5*CROSS*LOG(RTOP*RTOP+YROT*YROT) TAGENT2=-0.5*CROSS*LOG(RLOW*RLOW+YROT*YROT) TAGENT=TAGENT1-TAGENT2 TAGENTA=-3.141592654*DOT+TAGENT IF(ABS(XROT).GE.0.00000001) VALUEL=TAGENT DISTANCE=0.5*ELENGTH(NEL)-ABS(XROT) IF(DISTANCE.GE.0.AND.ICOL.GT.NELM)VALUEL=1000000 ELSE VALUEL=VL1(DOT,CROSS,RTOP,YROT)-VL1(DOT,CROSS,RLOW,YROT) ENDIF RETURN END C**************************************************** FUNCTION VL1(DOT,CROSS,V,Y0) c IMPLICIT REAL*8 (A-H,O-Z) VL1=-1.*DOT*ATAN(V/Y0)-0.5*CROSS*LOG(V*V+Y0*Y0) RETURN END CC************************************************************* C MKERNEL SUBROUTINE CAN CALCULATE THE M(s,x) INTEGRATION * C INCLUDING SINGULAR AND NONSINGULAR VALUE * C BY CLOSED FORM SOLUTION * C************************************************************** SUBROUTINE MKERNEL(ICOL,NEL,XROT,YROT,ELENGTH,VALUEM) c IMPLICIT REAL*8 (A-H,O-Z) COMMON /NORMAL/ XNORMAL,YNORMAL,THETA ! COMMON/BAR/ XBAR,YBAR COMMON/AA/ NELM,NNODE DIMENSION ELENGTH(1000) * ,XNORMAL(1000),YNORMAL(1000) * ,THETA(1000) save RTOP=0.5*ELENGTH(NEL)-XROT RLOW=-0.5*ELENGTH(NEL)-XROT XN=XNORMAL(ICOL) YN=YNORMAL(ICOL) XS=XNORMAL(NEL) YS=YNORMAL(NEL) DOT=XN*XS+YN*YS CROSS=-1.*XN*YS+YN*XS C***************************************************** C FOR SINGULAR ELEMENT * C***************************************************** IF(ABS(YROT).LE.0.000001) THEN VALUEM=DOT*(1.0/RTOP-1.0/RLOW) ELSE C************MACSYMA INTEGRATION ************************ C FOR REGULAR ELEMENT * C************UP AND LOWER !****************************** UP=DOT*RTOP/(RTOP*RTOP+YROT*YROT) * -YROT*CROSS/(RTOP*RTOP+YROT*YROT) DOWN=DOT*RLOW/(RLOW*RLOW+YROT*YROT) * -YROT*CROSS/(RLOW*RLOW+YROT*YROT) VALUEM=UP-DOWN ENDIF RETURN END C******************************************************** C**** SUBROUTINE GEOM PREPROCESSOR PROGRAM BY UNV FILE* C**** GEOMTRY DATA PROCESSING OF X,Y,Z COORDINATE * C**** AND CONNECTIVITY * C******************************************************** SUBROUTINE GEOM(X,Y,INC,NINDEX,NELM,NNODE) c IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(1000),Y(1000),Z(1000),NINDEX(1000) DIMENSION INC(2,1000) CHARACTER*80 TITLE0,TITLE1 PRINT *,'INPUT NELM,NNODE' READ(5,*) NELM,NNODE 121 FORMAT(36A) 101 READ(15,501)TITLE0 IF(TITLE0(4:6).NE.' 15')GO TO 101 !¨ú³¡¥÷¦r¦ê C*****READ THE CORNER POINT COORDINATES X,Y,Z ********************** DO 7001 K=1,NNODE READ(15,100,END=7001) M1,M2,M3,M4,PP,Q,R X(M1)=PP Y(M1)=Q Z(M1)=R 100 FORMAT(4I10,3E13.5) 501 FORMAT(80A) !A:INPUT ¤å¦r«¬¸ê®Æ 7001 CONTINUE 202 READ(15,501)TITLE1 IF(TITLE1(4:6).NE.' 71')GO TO 202 C*****READ THE INCIDENCE OF EACH ELEMENT *************************** DO 8001 J=1,NELM READ(15,800) II,JJ,KK,LL,MM,NN,IL NINDEX(J)=II 800 FORMAT(7I10) READ(15,200,END=8001) NINC1,NINC2 INC(1,II)=NINC1 INC(2,II)=NINC2 200 FORMAT(2I10) 8001 CONTINUE RETURN END C C C C******************************************************* C SUBROUTINE GRULE * C THIS IS A SUBROUTINE TO COMPUTE THE (N+1)/2 * C NONNEGATIVE AB-SCISSAS X(I),W(I) OF THE N-PT * C GAUSSIAN-LENGENDRE INTEGRATION RULE NORMALIZED * C TO THE INTERVAL [-1,1] * C******************************************************* SUBROUTINE GRULE(N,X,W) c IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(N),W(N) M=(N+1)/2 E1=N*(N+1) DO 1 I=1,M T=(4.0*I-1.0)*3.1415926536/(4*N+2) X0=(1.-(1.-1./N)/(8.*N*N))*COS(T) PKM1=1.0 PK=X0 DO 3 K=2,N T1=PK*X0 PKP1=T1-PKM1-(T1-PKM1)/K+T1 PKM1=PK PK=PKP1 3 CONTINUE DEN=1.-X0*X0 D1=N*(PKM1-X0*PK) DPN=D1/DEN D2PN=(2.*X0*DPN-E1*PK)/DEN D3PN=(4.*X0*D2PN+(2.-E1)*DPN)/DEN D4PN=(6.*X0*D3PN+(6.-E1)*D2PN)/DEN U=PK/DPN V=D2PN/DPN H=-U*(1.+.5*U*(V+U*(V*V-U*D3PN/(3.*DPN)))) P=PK+H*(DPN+0.5*H*(D2PN+H/3.*(D3PN+.25*H*D4PN))) DP=DPN+H*(D2PN+.5*H*(D3PN+H*D4PN/3.)) H=H-P/DP X(I)=X0+H FX=D1-H*E1*(PK+.5*H*(DPN+H/3.*(D2PN+.25*H*(D3PN+.2*H*D4PN)))) W(I)=2.*(1.-X(I)*X(I))/(FX*FX) 1 CONTINUE IF(M+M.GT.N)X(M)=0. RETURN END C C C C*********************************************** C SUBROUTINE ARANGE * C THE SUBROUTINE IS TO ARANGE THE GAUSSIAN * C QUA COORD. AND WEIGHTING FUNCTION * C*********************************************** SUBROUTINE ARANGE(N,X,W) c IMPLICIT REAL*8 (A-H,O-Z) REAL X(N),W(N) M=N/2 DO 1 I=1,M X(N-I+1)=-X(I) W(N-I+1)=W(I) 1 CONTINUE RETURN END C************************************************ C SUBROUTINE LUPPDC MATRIX DECOMPOSITION * C************************************************ SUBROUTINE LUPPDC(A,L,N) c IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(1000,1000),L(N) C DIMENSION A(NN,N),L(N) DO 50 J=1,N J1=J-1 AJJ=0.0 DO 30 I=1,N I1=MIN0(I-1,J1) AIJ=A(I,J) IF(I1.LE.0) GO TO 15 DO 10 K=1,I1 10 AIJ=AIJ-A(I,K)*A(K,J) 15 IF(I.GE.J) GO TO 20 A(I,J)=AIJ/A(I,I) GO TO 30 20 A(I,J)=AIJ IF(AJJ.GE.ABS(AIJ)) GO TO 30 M=I AJJ=ABS(AIJ) 30 CONTINUE IF(AJJ.EQ.0.0) GO TO 55 L(J)=M IF(M.EQ.J) GO TO 50 DO 40 I=1,N AIJ=A(M,I) A(M,I)=A(J,I) 40 A(J,I)=AIJ 50 CONTINUE RETURN 55 PRINT 95,J,J 95 FORMAT(3H A(I3,1H,I3,6H)=0.0/) RETURN END C********************************************************* C SUBROUTINE LUPPSB SOLVE Ax=B * C********************************************************* SUBROUTINE LUPPSB(A,B,L,N) c IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A(1000,1000),B(N),L(N) C DIMENSION A(NN,N),B(N),L(N) DO 70 I=1,N I1=I-1 M=L(I) BI=B(M) B(M)=B(I) IF(I1.LE.0) GO TO 70 DO 60 K=1,I1 60 BI=BI-A(I,K)*B(K) 70 B(I)=BI/A(I,I) N1=N-1 DO 90 IN=1,N1 I=N-IN I1=I+1 BI=B(I) DO 80 K=I1,N 80 BI=BI-A(I,K)*B(K) 90 B(I)=BI RETURN END C********************************************************** C SUBROUTINE BC INPUT BOUNDARY CONDITION * C********************************************************** SUBROUTINE BC(BCU,BCT,VKNOWN) c IMPLICIT REAL*8 (A-H,O-Z) COMMON/AA/ NELM,NNODE DIMENSION BCU(1000),BCT(1000),VKNOWN(1000) PRINT *,'INPUT NU' READ(5,*) NU PRINT *,'INPUT NT' READ(5,*) NT DO 960 I=1,NU READ(11,*) K,BDU BCU(K)=BDU !K:¤º³¡«ü¼Ð WRITE(*,10) K,BCU(K) 960 CONTINUE DO 970 J=1,NT READ(12,*) K,BDT BCT(K)=BDT !K:¤º³¡«ü¼Ð WRITE(*,20) K,BCT(K) 970 CONTINUE DO 980 K=1,NELM READ(13,*) VKNOWN(K) 980 CONTINUE C WRITE(6,*) (VKNOWN(I),I=1,8) 10 FORMAT(1X,'BCU(',I3,')=',F10.5) 20 FORMAT(1X,'BCT(',I3,')=',F10.5) RETURN END