[文件读写] <边界元法计算界面几何性质 >的数据输入问题

!***************************************************************
!
!     CALCULATION OF THE GEOMETRICAL PROPERTIES OF SECTION (GPS)
!                   (boundary integration methon)
!                      SELINA999   2023-06-04
!***************************************************************
!
      IMPLICIT REAL*8(A-H,O-Z)
      CHARACTER*8 DATA
      CHARACTER TITLE(70)
      DIMENSION GPS(6),RES(6),ROC(101),X(101),Y(101)
!
!     NOA------------total number of part areas
!     SOA------------sign of the part area
!     NP-------------total number of nodes of the part area
!     X(101),Y(101)--node coordinates
!     NE-------------total number of boundary elements of the
!                    area
!     ROC(101)-------radiu of the boundary element (POSITIVE
!                    radiu for convex boundary, NEGATIVE radiu
!                    for concave boundary,'0'for straight line)
!     NDV------------number of divisions of each boundary element
!                    integration
!     RES(6),GPS(6)--GPS relative to global axis of the part and
!                    total area respectively,(1)A,(2)Sx,(3)Sy,
!                    (4)Ix,(5)Iy,(6)Ixy
!     XC,YC----------centriod coordinates
!     AIX,AIY,AIXY---Ix,Iy,Ixy relative to parallel centriod axis
!     PIMAX,PIMIN----centriodal principle moments of inertia
!     ANG------------declination of the principle axis relative to
!                    Imax
!
      WRITE(*,*) 'Input the file name of data please'
      READ(*,'(A)') DATA
      OPEN(5,FILE=DATA)
      OPEN(6,FILE='RESU',STATUS='NEW')
!
!        a.Input the governing data
!
      READ(5,2000) TITLE
      WRITE(6,2010) TITLE
      WRITE(*,2010) TITLE
!
      WRITE(*,2020)
      READ(*,*) NDV
      READ(5,*) NOA
      WRITE(6,2030) NOA
      WRITE(*,2030) NOA
!
      DO 10 I=1,6
   10 GPS(I)=0.0
!
      DO 90 M=1,NOA
!
!        Input part area data
!
      READ(5,*) SOA
      READ(5,*) NP
      WRITE(6,2040) M,SOA,NP
      WRITE(*,2040) M,SOA,NP
      DO 15 I=1,NP
      READ(5,*) X(I),Y(I),ROC(I)
      WRITE(6,2050) I,X(I),Y(I),I,ROC(I)
   15 WRITE(*,2050) I,X(I),Y(I),I,ROC(I)
      X(NP+1)=X(1)
      Y(NP+1)=Y(1)
!
      DO 20 I=1,6
   20 RES(I)=0.0
!
      NE=NP
      DO 70 N=1,NE
!
!          b.Determime the element constants
!
      DX=X(N+1)-X(N)
      DY=Y(N+1)-Y(N)
      CL=DSQRT(DX*DX+DY*DY)
!
      IF (ROC(N).EQ.0.) GOTO 30
!
!         For curve boundary element
!
      DCX=DY/CL
      DCY=-DX/CL
      XM=0.5*(X(N+1)+X(N))
      YM=0.5*(Y(N+1)+Y(N))
      R=DABS(ROC(N))
      S=ROC(N)/R
      D2=R*R-CL*CL/4.
      D=DSQRT(D2)
      XCC=XM-S*D*DCX
      YCC=YM-S*D*DCY
      ZETA=2.0*DASIN(0.5*CL/R)
      ALFA=DACOS(DCX)
      IF (DCY.GE.0.) GOTO 30
      ALFA=-ALFA
!
   30 X1=X(N)
        Y1=Y(N)
!
!         c.Integrate along the boundary element
!
      DO 60 I=1,NDV
      PDV=1.0/NDV
!
      IF (ROC(N).NE.0.) GOTO 40
      X2=X1+PDV*DX
      Y2=Y1+PDV*DY
      GOTO 50
!
   40 BI=I
      BETA=ALFA-S*0.5*ZETA+S*BI*PDV*ZETA
      X2=XCC+S*R*DCOS(BETA)
      Y2=YCC+S*R*DSIN(BETA)
!
!         Detetmine the integral point,length and direction cosine
!
   50 DDX=X2-X1
      DDY=Y2-Y1
      DL=DSQRT(DDX*DDX+DDY*DDY)
      DDCX=DDY/DL
      DDCY=-DDX/DL
!
      DXM=0.5*(X1+X2)
      DYM=0.5*(Y1+Y2)
      RES(1)=RES(1)+DXM*DDCX*DL
      RES(2)=RES(2)+0.5*DYM*DYM*DDCY*DL
      RES(3)=RES(3)+0.5*DXM*DXM*DDCX*DL
      RES(4)=RES(4)+DYM**3*DDCY*DL/3.
      RES(5)=RES(5)+DXM**3*DDCX*DL/3.
      RES(6)=RES(6)+0.5*DXM*DXM*DYM*DDCX*DL
      X1=X2
      Y1=Y2
   60 CONTINUE
   70 CONTINUE
!
!         Detemine GPS of total area
!
      DO 80 I=1,6
   80 GPS(I)=GPS(I)+RES(I)*SOA
   90 CONTINUE
!
      WRITE(6,2060) NDV
      WRITE(*,2060) NDV
!
!         d.Calculate GPS respect to parallel central axis
!
      A=GPS(1)
      XC=GPS(3)/A
      YC=GPS(2)/A
      AIX=GPS(4)-YC*YC*A
      AIY=GPS(5)-XC*XC*A
      AIXY=GPS(6)-XC*YC*A
!
!         e.Calculate the centriodal principle moments of inertia
!
      C1=0.5*(AIX+AIY)
      C2=0.5*(AIX-AIY)
      C3=DSQRT(C2*C2+AIXY*AIXY)
      PIMAX=C1+C3
      PIMIN=C1-C3
      IF (AIX.EQ.PIMIN) GOTO 95
      ANG=DATAN(-AIXY/(AIX-PIMIN))*180.0/3.1415926
      GOTO 100
   95 ANG=-90.0
  100 WRITE(6,2070) A,XC,YC,AIX,AIY,AIXY,PIMAX,PIMIN,ANG
      WRITE(*,2070) A,XC,YC,AIX,AIY,AIXY,PIMAX,PIMIN,ANG
!
!
2000 FORMAT(70A1)
2010 FORMAT(2X,70A1)
2020 FORMAT(/,2X,'please input the total number of divisions of &
          & each bounday element for integration  NDV')
2030 FORMAT(/,2X,'TOTAL NUMBER OF PART AREAS',6X,'NOA=',I4)
2040 FORMAT(//2X,'PART AREA NO:',I3,/,8X,'SIGN OF THE PART AREA &
     & SOA=',2X,F3.0,/,8X,'TOTAL NUMBER NODES      NP=',I4,//,&
     & 2X,'NODE',8X,'COORDINATES(MM)',8X,'ELEMENT',5X,'RADIU(MM)',/,&
     & 4X,'NO',10X,'X',9X,'Y',13X,'NO',10X,'R',/)
2050 FORMAT(3X,I3,5X,2F10.2,8X,I3,5X,F10.2)
2060 FORMAT(//,2X,'GEOMETRICAL PROPERTIES OF THE SECTION (NDV=',&
     & I4,')',//)
2070 FORMAT(4X,'TOTAL AREA',22X,'A=',E11.4,2X,'(MM**2)',//,&
     & 4X,'CENTRION COORDINATES',11X,'Xc=',E11.4,2X,'(MM)',/,&
     & 35X,'Yc=',E11.4,2X,'(MM)',//&
     & 4X,'MOMENT OF INTRIA RELATIVE TO THE PARALELL CENTRIOD AXIS',/,&
     &/,35X,'Ix=',E11.4,2X,'(MM**4)',/,35X,'Iy=',E11.4,2X,'(MM**4)',&
     &/,34X,'Ixy=',E11.4,2X,'(MM**4)',//,&
     & 4X,'CENTRIODAL PRINCIPLE MOMENTS OF INERTIA',/,&
     &/,33X,'Imax=',E11.4,2X,'(MM**4)',/,33X,'Imin=',E11.4,2X,'(MM**4)' &
     &,//,4X,'DECLINATION OF THE PRINCIPLE AXIS OF Imax',/,&
     &/,32X,'ANGLE=',E11.4,2X,'(DEG)')
      STOP
      END