PROGRAM HW1MAIN IMPLICIT NONE INTEGER II,JJ,KK,LL,NNODES,NUMEL,NUMGAUSS,JL(2) INTEGER LDLHS REAL*8 RADIUS, DSTEP,DIST,WK,PX,PY,PI,CIRCSTEP INTEGER CIRCNODES,LEVELS,INTNODES PARAMETER(INTNODES=450,CIRCNODES=50,LEVELS=9) C PARAMETER(CIRCSTEP=0.001,INTRADIUS,0.95) PARAMETER(RADIUS=1.0,WK=1.0) PARAMETER (NNODES=24,NUMEL=24,NUMGAUSS=2) PARAMETER(LDLHS=48) PARAMETER(PX=0.0,PY= - 0.5,PI=3.1415926535) INTEGER IN(NUMEL,2),IPIV(2*NNODES) REAL*8 XYCOORD(NNODES,2),INFO REAL*8 XL(2),YL(2),ALPHA(NNODES),THETA(NNODES) REAL*8 AMAT(NNODES,NNODES), BMAT(NNODES,NNODES),DELTA REAL*8 CMAT(NNODES,NNODES),DMAT(NNODES,NNODES),EVEC(NNODES) REAL*8 ZED,PHI(2),ZETA(2) REAL*8 XDELTA,YDELTA,RI REAL*8 DRDN,GI,DGDN REAL*8 LHS(2*NNODES,2*NNODES),RHS(2*NNODES) REAL*8 USOLN(NNODES),DUDN(NNODES) REAL*8 INTRADIUS,STARTRADIUS REAL*8 XINT(INTNODES),YINT(INTNODES) REAL*8 UFIN(NNODES+INTNODES),F(NNODES) REAL*8 XFIN(NNODES+INTNODES),YFIN(NNODES+INTNODES) REAL*8 UINT(INTNODES),THETACIRC(CIRCNODES) REAL*8 DGDX,DGDY,DDGDXDN,DDGDYDN REAL*8 DUINTX(INTNODES),DUINTY(INTNODES) REAL*8 DUINT(INTNODES),ANGDUINT(INTNODES) CIRCSTEP=0.1 STARTRADIUS=0.0 STARTRADIUS=0.0+0.95 INTRADIUS=INTRADIUS+STARTRADIUS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC INITIALIZE VALUES CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO II=1,INTNODES DUINTX(II)=0.0 DUINTY(II)=0.0 DUINT(II)=0.0 ENDDO DGDX=0.0 DDGDXDN=0.0 ANGDUINT=0.0 DDGDYDN=0.0 ZETA(1)=0.0 ZETA(2)=0.0 ZETA(1)= ZETA(1) - 0.5773502691 ZETA(2)= ZETA(2) + 0.5773502691 DO II=1,NUMEL IN(II,1)=0 IN(II,2)=0 ENDDO DO II=1,NNODES USOLN(II)=0.0 DUDN(II)=0.0 DO JJ=1,2 IN(II,JJ)=0 XYCOORD(II,JJ)=0.0 ENDDO ALPHA(II)=0.0 THETA(II)=0.0 EVEC(II)=0.0 F(II)=0.0 DO JJ=1,NNODES AMAT(II,JJ)=0.0 BMAT(II,JJ)=0.0 CMAT(II,JJ)=0.0 DMAT(II,JJ)=0.0 ENDDO ENDDO DO II=1,2 XL(II)=0.0 YL(II)=0.0 JL(II)=0 PHI(II)=0.0 ENDDO DO II=1,2*NNODES RHS(II)=0.0 IPIV(II)=0 DO JJ=1,2*NNODES LHS(II,JJ)=0.0 ENDDO ENDDO ZED=0.0 DELTA=0.0 DRDN=0.0 GI=0.0 DGDN=0.0 DSTEP=0.0 DO II=1,INTNODES XINT(II)=0.0 YINT(II)=0.0 UINT(II)=0.0 ENDDO DO II=1,(NNODES+INTNODES) WRITE(*,*) "CHECK POINTS",II UFIN(II)=0.0 XFIN(II)=0.0 YFIN(II)=0.0 ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CREATE CIRCLE BOUNDARY CCC CCC AND CCC CCC START UP DATA CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DSTEP=360.0/REAL(NNODES) DO II=1,NNODES XYCOORD(II,1)=RADIUS*DCOS(DSTEP*REAL(II-1)*(PI/180.0)) XYCOORD(II,2)=RADIUS*DSIN(DSTEP*REAL(II-1)*(PI/180.0)) C IN DEGREES ALPHA(II)=(180-DSTEP) C IN RADIANS THETA(II)=DSTEP*REAL(II-1)*(PI/180.0) C WRITE(*,*) "THETA",THETA(II) ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC CREATE INCIDENCE LIST CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO II=1,NUMEL-1 IN(II,1)=II IN(II,2)=II+1 C WRITE(*,*) IN(II,1),IN(II,2) ENDDO IN(NUMEL,1)=NUMEL IN(NUMEL,2)=1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC MAJOR LOOP OVER ALL NODES CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO II=1,NNODES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOOP OVER ALL ELEMENTS CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO LL=1,NUMEL XL(1)=XYCOORD(IN(LL,1),1) XL(2)=XYCOORD(IN(LL,2),1) YL(1)=XYCOORD(IN(LL,1),2) YL(2)=XYCOORD(IN(LL,2),2) JL(1)=IN(LL,1) JL(2)=IN(LL,2) DELTA=DSQRT((XL(2) - XL(1))**2 + (YL(2) - YL(1))**2) CCC PERFORM ANALYTIC INTEGRATION IF II EQUAL TO CCC NODES 1 AND 2 OF ELEMENT C IF(II.EQ.JL(1).OR.II.EQ.JL(2))THEN IF(II.EQ.JL(1))THEN BMAT(JL(1),JL(1))=BMAT(JL(1),JL(1))+(DELTA/2.0)*((3.0/2.0)- + DLOG(DELTA)) BMAT(JL(1),JL(2))=BMAT(JL(1),JL(2))+(DELTA/2.0)*((1.0/2.0)- + DLOG(DELTA)) C ENDIF ELSEIF(II.EQ.JL(2))THEN BMAT(JL(2),JL(2))=BMAT(JL(2),JL(2))+(DELTA/2.0)*((3.0/2.0)- + DLOG(DELTA)) BMAT(JL(2),JL(1))=BMAT(JL(2),JL(1))+(DELTA/2.0)*((1.0/2.0)- + DLOG(DELTA)) ELSE C ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOOP OVER ALL GAUSS PTS CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO KK=1,NUMGAUSS ZED=ZETA(KK) CCCCCCCCCCCCCCCCC DEFINE BASIS CCCCCCCCCCCCCCCCCCCCCCCCCC PHI(1) = (1.0-ZED)/2.0 PHI(2) = (1.0+ZED)/2.0 CCCCCCCCCCCCCCCC (X,Y) GAUSS POINTS CCCCCCCCCCCCCCCCCCCCC XDELTA=XL(1)*PHI(1) + XL(2)*PHI(2) YDELTA=YL(1)*PHI(1) + YL(2)*PHI(2) RI=DSQRT((XDELTA-XYCOORD(II,1))**2 + + (YDELTA-XYCOORD(II,2))**2) DRDN = ((YL(2)-YL(1))*(XDELTA-XYCOORD(II,1)) - + (XL(2)-XL(1))*(YDELTA-XYCOORD(II,2)))/ + (DELTA*RI) GI = -DLOG(RI) DGDN = -(1.0/RI)*(DRDN) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOOP OVER COLUMNS CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO JJ=1,2 AMAT(II,JL(JJ)) = AMAT(II,JL(JJ)) + PHI(JJ)* + DGDN*(DELTA/2.0)*WK BMAT(II,JL(JJ)) = BMAT(II,JL(JJ)) + PHI(JJ)* + GI*(DELTA/2.0)*WK ENDDO CCCCCCCCCCCCCCCCCC END GAUSS PTS LOOP CCCCCCCCCCCCCCCCCCCCC ENDDO ENDIF CCCCCCCCCCCCCCCCCC END ELEMENT LOOP CCCCCCCCCCCCCCCCCCCCCCC ENDDO DIST=DSQRT((PX - XYCOORD(II,1))**2 + + (PY - XYCOORD(II,2))**2) CCCC CALCULATE FORCING USING DISTANCE OF NODE II FROM SOURCE F(II) = (-DLOG(DIST)) AMAT(II,II) = AMAT(II,II) + ALPHA(II) C WRITE(*,*) BMAT(II,1) CCCCCCCCCCCCCC TYPE I BOUNDAY CONDITIONS CCCCCCCCCCCCCCCCCCCC IF(II.EQ.1.OR.II.GT.(NNODES/2).AND.II.LE.NNODES)THEN CMAT(II,II)=1.0 C WRITE(*,*) "CHECK",II EVEC(II)=0.0 CCCCCCCCCCCCCCCC TYPE II BOUNDAY CONDITIONS CCCCCCCCCCCCCCCCC ELSEIF(II.GE.2.AND.II.LE.(NNODES/2))THEN DMAT(II,II)=1.0 EVEC(II)=2.0*(DSIN((PI/2.0)-THETA(II))) C WRITE(*,*) EVEC(II) ENDIF CCCCCCCCCCCCCCCCCCCC END MAIN NODE LOOP CCCCCCCCCCCCCCCCCCCCCC ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC FILL IN FINAL GENERAL [A]X=B CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCC FILL IN LHS MATRIX CCCCCCCCCCCCCCCCCCCCCCCCCC DO II=1,NNODES DO JJ=1,NNODES LHS(II,JJ)= AMAT(II,JJ) LHS(II,JJ+NNODES)= -(BMAT(II,JJ)) LHS(II+NNODES,JJ)= CMAT(II,JJ) LHS(II+NNODES,JJ+NNODES)= DMAT(II,JJ) ENDDO CCCCCCCCCCCCCCCCCC FILL IN RHS VECTOR CCCCCCCCCCCCCCCCCCCCCCCC RHS(II)=F(II) RHS(II+NNODES)=EVEC(II) ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC SOLVE USING DGESV [A]X=B CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALL DGESV(LDLHS,1,LHS,LDLHS,IPIV,RHS,LDLHS,INFO) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC FINAL SOLUTION RHS CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO II=1,NNODES WRITE(*,*) "NODES",II,"VAL",RHS(II) USOLN(II) = RHS(II) C WRITE(*,*) "BC NODES",II,"SOLN", USOLN(II) DUDN(II) = RHS(II+NNODES) ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC SECTION FOR COMPUTING INTERIOR CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INTRADIUS=INTRADIUS + 0.900 C DSTEP=360.0/REAL(NNODES) C DO II=1,NNODES C XYCOORD(II,1)=RADIUS*DCOS(DSTEP*REAL(II-1)*(PI/180.0)) C XYCOORD(II,2)=RADIUS*DSIN(DSTEP*REAL(II-1)*(PI/180.0)) C THETA(II)=DSTEP*REAL(II-1)*(PI/180.0) DSTEP=0.0 DSTEP=360.0/REAL(CIRCNODES) DO II=1,CIRCNODES THETACIRC(II)=0.0 THETACIRC(II)=DSTEP*REAL(II-1)*(PI/180.0) C WRITE(*,*) THETACIRC(II) ENDDO C WRITE(*,*) "INTERNAL RADIUS",INTRADIUS C WRITE(*,*) CIRCSTEP DO JJ=0,LEVELS-1 DO II=1,CIRCNODES XINT(II+JJ*CIRCNODES)=INTRADIUS*DCOS(DSTEP*REAL(II-1) + *(PI/180.0)) YINT(II+JJ*CIRCNODES)=INTRADIUS*DSIN(DSTEP*REAL(II-1) + *(PI/180.0)) C XINT(II+JJ*CIRCNODES)=INTRADIUS*DCOS(THETACIRC(II)) C YINT(II+JJ*CIRCNODES)=INTRADIUS*DCOS(THETACIRC(II)) C WRITE(*,*) "YINT",YINT(II+JJ*CIRCNODES),"N",(II+JJ*CIRCNODES) ENDDO WRITE(*,*) "FIRST",INTRADIUS INTRADIUS=STARTRADIUS-(REAL(JJ+1)*REAL(CIRCSTEP)) WRITE(*,*) "INTERNAL RADIUS",INTRADIUS,"CC", + REAL((JJ+1))*REAL(CIRCSTEP) ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC MAJOR LOOP OVER ALL NODES CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO II=1,INTNODES UINT(II)=0.0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOOP OVER ALL ELEMENTS CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO LL=1,NUMEL XL(1)=XYCOORD(IN(LL,1),1) XL(2)=XYCOORD(IN(LL,2),1) YL(1)=XYCOORD(IN(LL,1),2) YL(2)=XYCOORD(IN(LL,2),2) JL(1)=IN(LL,1) JL(2)=IN(LL,2) DELTA=DSQRT((XL(2) - XL(1))**2 + (YL(2) - YL(1))**2) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC LOOP OVER ALL GAUSS PTS CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO KK=1,NUMGAUSS ZED=ZETA(KK) CCC DEFINE BASIS PHI(1) = (1.0-ZED)/2.0 PHI(2) = (1.0+ZED)/2.0 CCC (X,Y) GAUSS POINTS XDELTA=XL(1)*PHI(1) + XL(2)*PHI(2) YDELTA=YL(1)*PHI(1) + YL(2)*PHI(2) RI=DSQRT((XDELTA-XINT(II))**2 + + (YDELTA-YINT(II))**2) DRDN = ((YL(2)-YL(1))*(XDELTA-XINT(II)) - + (XL(2)-XL(1))*(YDELTA-YINT(II)))/ + (DELTA*RI) C ADDED *(1.0/(2.0*PI)) GI = -DLOG(RI) DGDN = -(1.0/RI)*(DRDN) DO JJ=1,2 UINT(II) = UINT(II) + + ((DUDN(JL(JJ))*PHI(JJ)*GI) - + (USOLN(JL(JJ))*PHI(JJ)*DGDN))*((DELTA/2.0)*WK) ENDDO CCC CALCULATE DERIVATIVES OF U CCC TOOK AWAY THE DSQRT AT THE DENOMINATOR DGDX= (1.0/(2.0*PI))*(XDELTA-XINT(II))/ + DSQRT((XDELTA-XINT(II))**2+(YDELTA-YINT(II))**2) DGDY= (1.0/(2.0*PI))*(YDELTA-YINT(II))/ + DSQRT((YDELTA-YINT(II))**2+(XDELTA-XINT(II))**2) C DDGDXDN= ((YL(1)-YL(2))/(DELTA*RI)) + C + (((YL(2)-YL(1))*(XDELTA-XINT(II)) - C + (XL(2)-XL(1))*(YDELTA-YINT(II)))/ C + (DELTA*((XDELTA-XINT(II))**2+(YDELTA-YINT(II))**2) C + **(3.0/2.0)))* C + (XDELTA-XINT(II)) c originals DDGDXDN=((XINT(II) - XDELTA)/((XDELTA-XINT(II))**2 + +(YDELTA-YINT(II))**2)**(3.0/2.0))*DRDN*(1/(2.0*PI*RI)) + + (((YL(2) - YL(1))/(DELTA*RI)) + + ((XINT(II) - XDELTA)/ + ((XDELTA-XINT(II))**2+(YDELTA-YINT(II))**2))*DRDN)* + (1.0/(2.0*PI*RI)) DDGDYDN=((YINT(II) - YDELTA)/((YDELTA-YINT(II))**2 + +(XDELTA-XINT(II))**2)**(3.0/2.0))*DRDN*(1/(2.0*PI*RI)) + + (((XL(2) - XL(1))/(DELTA*RI)) + + ((YINT(II) - YDELTA)/ + ((YDELTA-YINT(II))**2+(XDELTA-XINT(II))**2))*DRDN)* + (1.0/(2.0*PI*RI)) C DGDX=((XDELTA-XINT(II))/ C + (((XINT(II)-XDELTA)**2+(YINT(II)-YDELTA)**2)**(1.5))) C + *DRDN + C + (1.0/RI)*((YL(2)-YL(1))/(DELTA*RI)) + C + ((XDELTA-XINT(II))/(RI**2))*DRDN C DGDX=((YDELTA-INT(II))/ C + (((YINT(II)-YDELTA)**2+(XINT(II)-XDELTA)**2)**(1.5))) C + *DRDN + C + (1.0/RI)*((XL(2)-XL(1))/(DELTA*RI)) + C + ((YDELTA-YINT(II))/(RI**2))*DRDN C DDGDYDN= ((XL(1)-XL(2))/(DELTA*RI)) + C + (((XL(2)-XL(1))*(XDELTA-YINT(II)) - C + (YL(2)-YL(1))*(XDELTA-XINT(II)))/ C + (DELTA*((YDELTA-YINT(II))**2+(XDELTA-XINT(II))**2) C + **(3.0/2.0)))* C + (YDELTA-YINT(II)) DO JJ=1,2 DUINTX(II) = DUINTX(II) + + ((DUDN(JL(JJ))*PHI(JJ)*DGDX) - + (USOLN(JL(JJ))*PHI(JJ)*DDGDXDN))*((DELTA/2.0)*WK) ENDDO DO JJ=1,2 DUINTY(II) = DUINTY(II) + + ((DUDN(JL(JJ))*PHI(JJ)*DGDY) - + (USOLN(JL(JJ))*PHI(JJ)*DDGDYDN))*((DELTA/2.0)*WK) ENDDO C ANGDUINT(II)=DATAN2(DUINTY(II)/DUINTX(II)) C DUINT(II) = DUINTY(II)/(DSIN(ANGDUINT(II))) CCCCCCCCCCCCCCCCCC END GAUSS LOOP CCCCCCCCCCCCCCCCCCCCCCCCCCC ENDDO CCCCCCCCCCCCCCCCCC END ELEMENT LOOP CCCCCCCCCCCCCCCCCCCCCCCCC ENDDO CCCCCCCCCCCCCCCCCC END INTERIOR POINT LOOP CCCCCCCCCCCCCCCCCC DIST=DSQRT((PX - XINT(II))**2 + + (PY - YINT(II))**2) UINT(II) = UINT(II) + (-DLOG(DIST)) DUINTX(II) = DUINTX(II) + (PX-XINT(II))/(2*PI*DIST*DIST) DUINTY(II) = DUINTY(II) + (PY-YINT(II))/(2*PI*DIST*DIST) ENDDO CCCCCCCCCCCCCCCCCC PREPARE ALL DATA FOR VISUALIZATION CCCCCCC DO II=1,NNODES UFIN(II)= USOLN(II) XFIN(II) = XYCOORD(II,1) YFIN(II) = XYCOORD(II,2) ENDDO DO II=1,INTNODES UFIN(II+NNODES)=UINT(II) XFIN(II+NNODES) = XINT(II) YFIN(II+NNODES) = YINT(II) ENDDO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC WRITING SOME RESULTS TO FILE CCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCC WRITE [X,Y,U] OF ALL NODES TO FILE OPEN (UNIT = 7, FILE = "hw1allnodesoln.dat", STATUS="unknown") DO II=1,(NNODES+INTNODES) WRITE(7,*) XFIN(II),YFIN(II),UFIN(II) ENDDO CLOSE(7) CCC WRITE RHS SIDE SOLUTION AFTER FORWARD PROBLEM TO FILE OPEN (UNIT = 7, FILE = "hw1rhs.dat", STATUS="unknown") DO II=1,2*NNODES WRITE(7,*) RHS(II) ENDDO CLOSE(7) CCC WRITE X AND Y COORDINATES OF BOUNDARY ONLY TO FILE OPEN (UNIT = 7, FILE = "hw1coord.dat", STATUS="unknown") DO II=1,NNODES WRITE(7,*) XYCOORD(II,1),XYCOORD(II,2) ENDDO CLOSE(7) WRITE(*,*) ZETA(1),ZETA(2) CCC BCS TO FILE OPEN (UNIT = 7, FILE = "hw1belnodesoln.dat", STATUS="unknown") DO II=1,NNODES WRITE(7,*) USOLN(II) ENDDO CLOSE(7) CCC DERIVATIVE OF U OPEN (UNIT = 7, FILE = "hw1derivativeuinterior.dat", + STATUS="unknown") DO II=1,INTNODES WRITE(7,*) XINT(II),YINT(II),DUINTX(II),DUINTY(II) ENDDO CLOSE(7) END