convert fortran to matlab

6 次查看（过去 30 天）
mahdi shabani
mahdi shabani 2015-1-11
回答： Steven Lord 2016-7-11
hello Friends
I have a code based on Fortran and need to compile it to Matlab Format as soon as possible
any body can help me on this?
! *****************************************************************
! PROGRAM FOR THE SOLUTION OF 2D FRAME WITH 3 DEGREES OF FREEDOM
! *****************************************************************
PROGRAM FRAME
IMPLICIT REAL*8(A-H,O-Z)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
!
! THIS ROUTINE CONTROLS THE CALLING,IN ORDER,
! OF ALL SUBROUTINES
!
CHARACTER*6 TITLE(12)
OPEN(6,FILE='FRAME1.RES',STATUS='UNKNOWN',ACCESS='SEQUENTIAL', &
FORM='FORMATTED')
CLOSE(6,STATUS='DELETE')
OPEN(5,FILE='INPUT_FILE.DATA',STATUS='UNKNOWN',ACCESS='SEQUENTIAL'&
,FORM='FORMATTED')
OPEN(6,FILE='FRAME_RESULT.RES',STATUS='UNKNOWN',&
ACCESS='SEQUENTIAL',FORM='FORMATTED')
WRITE(*,5)
5 FORMAT (/,5X,57HPLEASE PLACE YOUR INPUT FILE WITH NAME OF &
INPUT_FILE.DATA,/,5X,37HIN THE SAME FOLDER AS YOUR SOURCE IS.,/)
CALL DATA
OPEN(1,STATUS='UNKNOWN',ACCESS='SEQUENTIAL',FORM='FORMATTED')
CALL STIFFB
CALL ASSEMB
CALL GREDUC
CALL BAKSUB
CALL FORCE
CALL RESULT
WRITE(*,6)
6 FORMAT (/,5X,47HRESULTS ARE WRITTEN IN FRAME_RESULT.RES LOCATED &
,/,5X,37HIN THE SAME FOLDER AS YOUR SOURCE
CLOSE (1)
CLOSE(5)
CLOSE(6)
STOP
END
! ****************************************************************
! THIS SUBROUTINE READS & WRITES USER INFO ABOUT STRUCTURE
! ****************************************************************
! SPECIFICATION PART
SUBROUTINE DATA
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ICODE(3),PRESC(3)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40) &
ASTIF(80,80),ASLOD(80),REACT(80)
! THESE SPECIFICATIONS ARE FOR DYNAMIC DIMENSIONING
! ALLOCATE(PROPS(NMATS,NPROP),COORD(NPOIN,NDOFN),LNODS(NELEM,NDOFN)
! ,IFPRE(6),FIXED(NSVAB),RLOAD(NPOIN,NDOFN),ELOAD(NELEM,NEVAB),
! MATNO(2),STRES(NELEM,1),XDISP(NSVAB),TDISP(NPOIN,NDOFN),
! TREAC(NPOIN,NDOFN),ASTIF(NSVAB,NSVAB),ASLOD(NSVAB),
! REACT(NSVAB))
CHARACTER*6 TITLE(12)
! READ AND WRITE THE PROBLEM TITLE
!
READ(5,915) TITLE
WRITE(6,915) TITLE
915 FORMAT(12A6)
!
! READ AND WRITE THE CONTROL DATA
!
READ(5,*) NPOIN,NELEM,NBOUN,NMATS,NPROP,NNODE,NDOFN,NDIME, &
NSTRE
900 FORMAT(16I5)
WRITE(6,905) NPOIN,NELEM,NBOUN,NMATS,NPROP,NNODE,NDOFN,NDIME, &
NSTRE
905 FORMAT(//1X,7HNPOIN =,I5,3X,7HNELEM =,I5,3X,7HNBOUN =,I5,3X, &
7HNMATS =,I5,//1X,7HNPROP =,I5,3X,7HNNODE =,I5,3X, &
7HNDOFN =,I5,3X,7HNDIME =I5,//1X,7HNSTRE =,I5)
NSVAB=NPOIN*NDOFN
NEVAB=NNODE*NDOFN
!
! READ AND WRITE THE MATERIAL PROPERTIES
!
WRITE(6,950)
950 FORMAT(/,1H0,5X,19HMATERIAL PROPERTIES)
DO 10 IMATS=1,NMATS
READ(5,*) JMATS,(PROPS(JMATS,IPROP),IPROP=1,NPROP)
10 WRITE(6,910) JMATS,(PROPS(JMATS,IPROP),IPROP=1,NPROP)
910 FORMAT(I10,4F15.5)
!
! READ AND WRITE THE ELEMENT(MEMBER) NODAL
! (JOINT) CONNECTIONS
!
WRITE(6,960)
960 FORMAT(/,1H0,2X,2HEL,3X,5HNODES,3X,4HMAT.)
DO 20 IELEM=1,NELEM
READ(5,*) JELEM,(LNODS(JELEM,INODE),INODE=1,NNODE),MATNO(JELEM)
20 WRITE(6,920) JELEM,(LNODS(JELEM,INODE),INODE=1,NNODE),MATNO(JELEM)
920 FORMAT(5I5)
!
! READ AND WRITE THE NODAL(JOINT) COORDINATES
!
WRITE(6,970)
970 FORMAT(/,1H0,5X,4HNODE,5X,6HCOORD.)
DO 30 IPOIN=1,NPOIN
READ(5,*) JPOIN,(COORD(JPOIN,IDIME),IDIME=1,NDIME)
30 WRITE(6,930) JPOIN,(COORD(JPOIN,IDIME),IDIME=1,NDIME)
930 FORMAT(I10,3F15.5)
!
! READ AND WRITE THE BOUNDARY CONDITIONS
! AND STORE IN GLOBAL VECTORS
!
DO 40 ISVAB=1,NSVAB
IFPRE(ISVAB)=0
40 FIXED(ISVAB)=0.0
WRITE(6,980)
980 FORMAT(/,1H0,1X,28HRESTRAINED NODES,FIXITY CODE, &
22H AND PRESCRIBED VALUES)
IF(NBOUN.EQ.0) GO TO 55
DO 50 IBOUN=1,NBOUN
! BOUNDRY DATA IS READ AS FOLLOWS :
! 1-NUMBER OF NODE
! 2-ICODE IN FIRST DEGREE OF FREEDOM FOR THE NODE (1 FOR FIXED
! AND 0 FOR FREE RESTRAINT.
! 3-THE MAGNITUDE OF PRESCIBED DISPLACEMENT
! THIS PROCESS IS DONE FOR EACH NODE WITH RESTRAINTS.
READ(5,*) NODFX,(ICODE(IDOFN),PRESC(IDOFN),IDOFN=1,NDOFN)
WRITE(6,940) NODFX,(ICODE(IDOFN),PRESC(IDOFN),IDOFN=1,NDOFN)
940 FORMAT(I10,3(I5,F10.5))
DO 50 IDOFN=1,NDOFN
INDEX=(NODFX-1)*NDOFN+IDOFN
IFPRE(INDEX)=ICODE(IDOFN)
50 FIXED(INDEX)=PRESC(IDOFN)
55 CONTINUE
!
! READ AND WRITE THE NODAL(JOINT) APPLIED LOADS
!
WRITE(6,990)
990 FORMAT(/,1H0,5X,4HNODE,7X,11HNODAL LOADS)
DO 60 IPOIN=1,NPOIN
DO 60 IDOFN=1,NDOFN
60 RLOAD(IPOIN,IDOFN)=0.0
70 READ(5,*) IPOIN,(RLOAD(IPOIN,IDOFN),IDOFN=1,NDOFN)
WRITE(6,930) IPOIN,(RLOAD(IPOIN,IDOFN),IDOFN=1,NDOFN)
IF(IPOIN.LT.NPOIN) GO TO 70
! READ AND WRITE ELEMENT FIXED END MOMENT & SHEAR
WRITE(6,991)
991 FORMAT(/,1H0,5X,7HELEMENT,7X,15HELEMENTAL LOADS)
DO 5 IELEM=1,NELEM
DO 5 I=1,6
5 ELOAD(IELEM,I)=0.0
88 READ(5,*) IELEM,(ELOAD(IELEM,I),I=1,6)
WRITE(6,901) IELEM,(ELOAD(IELEM,I),I=1,6)
901 FORMAT(I10,6F11.3)
IF(IELEM.LT.NELEM)GO TO 88
RETURN
END
! ****************************************************************
! THIS SUBROUTINE READS & STORES ELEMENT STIFFNESS MATRICE
! ****************************************************************
! SPECIFICATION PART
SUBROUTINE STIFFB
IMPLICIT REAL*8(A-H,O-Z)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
DIMENSION ESTIF(6,6)
!
! EVALUATION OF MEMBER STIFFNESS MATRICES FOR PIN JOINTED
!
REWIND 1
DO 20 IELEM=1,NELEM
LPROP=MATNO(IELEM)
YOUNG=PROPS(LPROP,1)
! YOUNG=E
! XAREA=A
! DINERSI=I
XAREA=PROPS(LPROP,2)
DINERSI=PROPS(LPROP,3)
NODE1=LNODS(IELEM,1)
NODE2=LNODS(IELEM,2)
XPROJ=COORD(NODE2,1)-COORD(NODE1,1)
YPROJ=COORD(NODE2,2)-COORD(NODE1,2)
ELENG=SQRT(XPROJ*XPROJ+YPROJ*YPROJ)
SINTH=YPROJ/ELENG
COSTH=XPROJ/ELENG
! SINTH=SIN(X)
1 COSTH=COS(X)
X=YOUNG*XAREA/ELENG
Y=YOUNG*DINERSI/(ELENG**3)
! X=EA/L
! Y=EI/(L^3)
ESTIF(1,1)=X*(COSTH)**2+12*Y*(SINTH)**2
ESTIF(1,2)=X*(COSTH)*(SINTH)-12*Y*(COSTH)*(SINTH)
ESTIF(1,3)=-6*(YOUNG)*(DINERSI)*(SINTH)/(ELENG)**2
ESTIF(1,4)=-ESTIF(1,1)
ESTIF(1,5)=-ESTIF(1,2)
ESTIF(1,6)=ESTIF(1,3)
ESTIF(2,1)=ESTIF(1,2)
ESTIF(2,2)=X*(SINTH)**2+12*Y*(COSTH)**2
ESTIF(2,3)=6*(YOUNG)*(DINERSI)*(COSTH)/(ELENG)**2
ESTIF(2,4)=-ESTIF(1,2)
ESTIF(2,5)=-ESTIF(2,2)
ESTIF(2,6)=ESTIF(2,3)
ESTIF(3,1)=ESTIF(1,3)
ESTIF(3,2)=ESTIF(2,3)
ESTIF(3,3)=4*(YOUNG)*(DINERSI)/(ELENG)
ESTIF(3,4)=-ESTIF(1,3)
ESTIF(3,5)=-ESTIF(2,3)
ESTIF(3,6)=2*(YOUNG)*(DINERSI)/(ELENG)
ESTIF(4,1)=-ESTIF(1,1)
ESTIF(4,2)=-ESTIF(1,2)
ESTIF(4,3)=-ESTIF(1,3)
ESTIF(4,4)=-ESTIF(1,4)
ESTIF(4,5)=-ESTIF(1,5)
ESTIF(4,6)=-ESTIF(1,6)
ESTIF(5,1)=-ESTIF(2,1)
ESTIF(5,2)=-ESTIF(2,2)
ESTIF(5,3)=-ESTIF(2,3)
ESTIF(5,4)=-ESTIF(2,4)
ESTIF(5,5)=-ESTIF(2,5)
ESTIF(5,6)=-ESTIF(2,6)
ESTIF(6,1)=ESTIF(1,3)
ESTIF(6,2)=ESTIF(3,2)
ESTIF(6,3)=ESTIF(3,6)
ESTIF(6,4)=-ESTIF(1,3)
ESTIF(6,5)=-ESTIF(2,3)
ESTIF(6,6)=ESTIF(3,3)
DO 10000 I=1,6
DO 10000 J=1,6
10000 WRITE(1,*)ESTIF(I,J)
20 CONTINUE
RETURN
END
! ****************************************************************
! THIS SUBROUTINE READ ELEMENT LOADS DATA
! ****************************************************************
SUBROUTINE ELOAD
IMPLICIT REAL*8(A-H,O-Z)
WRITE(6,*)'GIVE ME TYPE OF ELEMENT LOADS FROM 1 TO 7'
DO IELEM=1,NELEM
DO I=1,6
ELOAD(IELEM,I)=0.0
END DO
END DO
DO N=1,(2*NELEM)
READ(5,*) IELEM,IFELOAD
IF(IELEM.EQ.0) EXIT
NODE1=LNODS(IELEM,1)
NODE2=LNODS(IELEM,2)
XPROJ=COORD(NODE2,1)-COORD(NODE1,1)
YPROJ=COORD(NODE2,2)-COORD(NODE1,2)
R=SQRT(XPROJ*XPROJ+YPROJ*YPROJ)
SINTH=YPROJ/R
COSTH=XPROJ/R
IF(IFELOAD.EQ.1) THEN
READ(5,*) W
ELOAD(IELEM,1)=ELOAD(IELEM,1)-(W*R)/2*SINTH
ELOAD(IELEM,2)=ELOAD(IELEM,2)+(W*R)/2*COSTH
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(W*R*R)/12
ELOAD(IELEM,4)=ELOAD(IELEM,4)-(W*R)/2*SINTH
ELOAD(IELEM,5)=ELOAD(IELEM,5)+(W*R)/2*COSTH
ELOAD(IELEM,6)=ELOAD(IELEM,6)-(W*R*R)/12
ELSE IF(IFELOAD.EQ.2) THEN
READ(5,'(F15.5)') W
ELOAD(IELEM,1)=ELOAD(IELEM,1)+0
ELOAD(IELEM,2)=ELOAD(IELEM,2)+W*XPROJ/2
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(W*XPROJ*XPROJ)/12
ELOAD(IELEM,4)=ELOAD(IELEM,4)+0
ELOAD(IELEM,5)=ELOAD(IELEM,5)+W*XPROJ/2
ELOAD(IELEM,6)=ELOAD(IELEM,6)-(W*XPROJ*XPROJ)/12
ELSE IF(IFELOAD.EQ.3) THEN
READ(5,'(F15.5)') W
ELOAD(IELEM,1)=ELOAD(IELEM,1)+0
ELOAD(IELEM,2)=ELOAD(IELEM,2)-W*YPROJ/2
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(W*YPROJ*YPROJ)/12
ELOAD(IELEM,4)=ELOAD(IELEM,4)+0
ELOAD(IELEM,5)=ELOAD(IELEM,5)-W*YPROJ/2
ELOAD(IELEM,6)=ELOAD(IELEM,6)-(W*YPROJ*YPROJ)/12
ELSE IF(IFELOAD.EQ.4) THEN
READ(5,'(F15.5)') W
ELOAD(IELEM,1)=ELOAD(IELEM,1)-(7*W*X)/20*SINTH
ELOAD(IELEM,2)=ELOAD(IELEM,2)+(7*W*X)/20*COSTH
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(W*X*X)/20
ELOAD(IELEM,4)=ELOAD(IELEM,4)-(3*W*X)/20*SINTH
ELOAD(IELEM,5)=ELOAD(IELEM,5)+(3*W*X)/20*COSTH
ELOAD(IELEM,6)=ELOAD(IELEM,6)-(W*X*X)/30
ELSE IF(IFELOAD.EQ.5) THEN
READ(5,'(F15.5)') W
ELOAD(IELEM,1)=ELOAD(IELEM,1)-(3*W*X)/20*SINTH
ELOAD(IELEM,2)=ELOAD(IELEM,2)+(3*W*X)/20*COSTH
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(W*X*X)/30
ELOAD(IELEM,4)=ELOAD(IELEM,4)-(7*W*X)/20*SINTH
ELOAD(IELEM,5)=ELOAD(IELEM,5)+(7*W*X)/20*COSTH
ELOAD(IELEM,6)=ELOAD(IELEM,6)-(W*X*X)/20
ELSE IF(IFELOAD.EQ.6) THEN
READ(5,'(3F15.5)') P,a,b
ELOAD(IELEM,1)=ELOAD(IELEM,1)-SINTH*(P*b*b*(b+3*a)/(X*X*X))
ELOAD(IELEM,2)=ELOAD(IELEM,2)+COSTH*(P*b*b*(b+3*a)/(X*X*X))
ELOAD(IELEM,3)=ELOAD(IELEM,3)+P*a*b*b/(X*X)
ELOAD(IELEM,4)=ELOAD(IELEM,4)+(-SINTH)*(P*a*a*(a+3*b)/(X*X*X))
ELOAD(IELEM,5)=ELOAD(IELEM,5)+(+COSTH)*(P*a*a*(a+3*b)/(X*X*X))
ELOAD(IELEM,6)=ELOAD(IELEM,6)-P*b*a*a/(X*X)
ELSE IF(IFELOAD.EQ.7) THEN
READ(5,'(3F15.5)') M,a,b
ELOAD(IELEM,1)=ELOAD(IELEM,1)+(-SINTH)*(6*M*a*b)/(X*X*X)
ELOAD(IELEM,2)=ELOAD(IELEM,2)+(+COSTH)*(6*M*a*b)/(X*X*X)
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(M*b)*(b-2*a)/(X*X)
ELOAD(IELEM,4)=ELOAD(IELEM,4)+(-SINTH)*(-6*M*a*b)/(X*X*X)
ELOAD(IELEM,5)=ELOAD(IELEM,5)+(+COSTH)*(-6*M*a*b)/(X*X*X)
ELOAD(IELEM,6)=ELOAD(IELEM,6)+(M*a)*(a-2*b)/(X*X)
ELSE IF(IFELOAD.EQ.8) THEN
LPROP=MATNO(IELEM)
YOUNG=PROPS(LPROP,1)
DINERSI=PROPS(LPROP,3)
READ(5,'(4F15.5)') o,H,T1,T2
ELOAD(IELEM,1)=ELOAD(IELEM,1)+0
ELOAD(IELEM,2)=ELOAD(IELEM,2)+0
ELOAD(IELEM,3)=ELOAD(IELEM,3)+(o*YOUNG*DINERSI/H)*(T2-T1)
ELOAD(IELEM,4)=ELOAD(IELEM,4)+0
ELOAD(IELEM,5)=ELOAD(IELEM,5)+0
ELOAD(IELEM,6)=ELOAD(IELEM,6)-(o*YOUNG*DINERSI/H)*(T2-T1)
! ****************************************************************
! THIS SUBROUTINE COMPUTES 3 NODAL INTERNAL FORCES.
! ****************************************************************
SUBROUTINE FORCE
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION ESTIF(6,6),FOMEM(6)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
!
! MEMBER FORCE CALCULATIONS FOR BOTH AXIAL
! BAR AND PIN-JOINTED PLANE FRAME PROBLEMS
!
REWIND 1
DO 10 IELEM=1,NELEM
READ(1,*) ESTIF
!
! EVALUATE THE MEMBER END FORCES
!
DO 20 IEVAB=1,NEVAB
FOMEM(IEVAB)=0.0
KOUNT=0
DO 20 INODE=1,NNODE
LOCAL=LNODS(IELEM,INODE)
DO 20 IDOFN=1,NDOFN
KOUNT=KOUNT+1
20 FOMEM(IEVAB)=FOMEM(IEVAB)+ESTIF(IEVAB,KOUNT)*TDISP(LOCAL,IDOFN)
FOMEM(1)=FOMEM(1)+ELOAD(IELEM,1)
FOMEM(2)=FOMEM(2)+ELOAD(IELEM,2)
FOMEM(3)=FOMEM(3)+ELOAD(IELEM,3)
FOMEM(4)=FOMEM(4)+ELOAD(IELEM,4)
FOMEM(5)=FOMEM(5)+ELOAD(IELEM,5)
FOMEM(6)=FOMEM(6)+ELOAD(IELEM,6)
!
! EVALUATE THE AXIAL FORCE
!
LPROP=MATNO(IELEM)
YOUNG=PROPS(LPROP,1)
XAREA=PROPS(LPROP,2)
DINERSI=PROPS(LPROP,3)
NODE1=LNODS(IELEM,1)
NODE2=LNODS(IELEM,2)
XPROJ=COORD(NODE2,1)-COORD(NODE1,1)
YPROJ=COORD(NODE2,2)-COORD(NODE1,2)
ELENG=SQRT(XPROJ*XPROJ+YPROJ*YPROJ)
SINTH=YPROJ/ELENG
COSTH=XPROJ/ELENG
STRES(IELEM,1)=FOMEM(1)*COSTH+FOMEM(2)*SINTH
STRES(IELEM,2)=FOMEM(2)*COSTH-FOMEM(1)*SINTH
STRES(IELEM,3)=FOMEM(3)
STRES(IELEM,4)=FOMEM(4)*COSTH+FOMEM(5)*SINTH
STRES(IELEM,5)=FOMEM(5)*COSTH-FOMEM(4)*SINTH
10 STRES(IELEM,6)=FOMEM(6)
RETURN
END
! ****************************************************************
! THIS SUBROUTINE OUTPUTS THE RESULTS
! ****************************************************************
SUBROUTINE RESULT
IMPLICIT REAL*8(A-H,O-Z)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
!
! WRITE THE NODAL(JOINT) DISPLACEMENTS AND REACTIONS
!
WRITE(6,900)
900 FORMAT(/,1H0,5X,4HNODE,1X,30HDISPLACEMENTS (DX,DY,ROTATION))
DO 10 IPOIN=1,NPOIN
10 WRITE(6,910) IPOIN,(TDISP(IPOIN,IDOFN),IDOFN=1,NDOFN)
910 FORMAT(5X,1I3,3F15.5)
WRITE(6,991)
991 FORMAT(/,1H0,5X,4HNODE,1X,19HREACTIONS (RX,RY,M))
DO 21 IPOIN=1,NPOIN
21 WRITE(6,990)IPOIN,(TREAC(IPOIN,IDOFN),IDOFN=1,NDOFN)
990 FORMAT(5X,1I3,3F15.5)
!
! WRITE THE ELEMENT(MEMBER) STRESSES(FORCES)
!
IF(NSTRE.EQ.0) GO TO 30
WRITE(6,920)
920 FORMAT(/,1H0,5X,7HELEMENT,6X,29HNODAL INTERNAL FORCES (N,V,M))
DO 20 IELEM=1,NELEM
WRITE(6,922) IELEM,(STRES(IELEM,ISTRE),ISTRE=1,3)
20 WRITE(6,923) (STRES(IELEM,ISTRE),ISTRE=4,6)
922 FORMAT(5X,I3,6X,3F15.5)
923 FORMAT(14X,3F15.5,/)
30 CONTINUE
RETURN
END
! ****************************************************************
! THIS THIS ROUTINE ASSEMBLES APPLIED LOADS AND FIXED END SHEAR & MOMENTS
! TO FORM THE GLOBAL FORCE VECTOR & ASSEMBELS ELEMENTAL STIFFNESS MATRIXS
! TO FORM STIFFNESS MATRICE.
! NODAL APPLIED LOADS = RLOAD(I)
! GLOBAL FORCE VECTOR = ELOAD
! LOCAL STIFF. MATRIX = ESTIF
! GLOBAL STIFF. MATRIX= ASTIF
! ****************************************************************
SUBROUTINE ASSEMB
IMPLICIT REAL*8(A-H,O-Z)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
DIMENSION ESTIF(6,6)
! THIS ROUTINE ASSEMBLES THE ELEMENT (MEMBER) STIFFNESSES AND APPLIED
! LOADS TO FORM THE GLOBAL STIFFNESS MATRIX AND FORCE VECTOR
REWIND 1
DO 10 ISVAB=1,NSVAB
ASLOD(ISVAB)=0.0
DO 10 JSVAB=1,NSVAB
ASTIF(ISVAB,JSVAB)=0.0
10 CONTINUE
!
! ASSEMBLE THE ELEMENT LOADS
!
DO 15 IPOIN=1,NPOIN
DO 15 IDOFN=1,NDOFN
NROWS=(IPOIN-1)*NDOFN+IDOFN
15 ASLOD(NROWS)=ASLOD(NROWS)+RLOAD(IPOIN,IDOFN)
DO 30 IELEM=1,NELEM
DO 1 I=1,6
DO 1 J=1,6
1 READ(1,*) ESTIF(I,J)
DO 20 INODE=1,NNODE
NODEI=LNODS(IELEM,INODE)
DO 20 IDOFN=1,NDOFN
NROWS=(NODEI-1)*NDOFN + IDOFN
NROWE=(INODE-1)*NDOFN + IDOFN
! THEN ELEMENT SAGGING MOMENT IN EACH DEGREE OF FREEDOM IS ADDED TO
! GLOBAL LOAD VECTOR (ASLOD) WITH NEGATIVE SIGN
ASLOD(NROWS)=ASLOD(NROWS) - ELOAD(IELEM,NROWE)
!
! ASSEMBLE THE ELEMENT STIFFNESS MATRICES
!
DO 20 JNODE = 1,NNODE
NODEJ=LNODS(IELEM,JNODE)
DO 20 JDOFN =1,NDOFN
NCOLS=(NODEJ-1)*NDOFN + JDOFN
NCOLE=(JNODE-1)*NDOFN + JDOFN
ASTIF(NROWS,NCOLS)=ASTIF(NROWS,NCOLS)+ESTIF(NROWE,NCOLE)
20 CONTINUE
30 CONTINUE
RETURN
END
! ****************************************************************
! THIS SUBROUTINE REDUCES THE GLOBAL STIFFNESS EQUATIONS BY
! DIRECT GAUSSIAN ELIMINATION
! ****************************************************************
SUBROUTINE GREDUC
IMPLICIT REAL*8(A-H,O-Z)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
!
! THIS ROUTINE REDUCES THE GLOBAL STIFFNESS
! EQUATIONS BY DIRECT GAUSSIAN ELIMINATION
!
NEQNS=NSVAB
DO 50 IEQNS=1,NEQNS
IF(IFPRE(IEQNS).EQ.1) GO TO 30
!
! REDUCE EQUATIONS
!
PIVOT=ASTIF(IEQNS,IEQNS)
IF(ABS(PIVOT).LT.1.0E-10) GO TO 60
IF(IEQNS.EQ.NEQNS) GO TO 50
IEQN1=IEQNS+1
DO 20 IROWS=IEQN1,NEQNS
FACTR=ASTIF(IROWS,IEQNS)/PIVOT
IF(FACTR.EQ.0.0) GO TO 20
DO 10 ICOLS=IEQNS,NEQNS
ASTIF(IROWS,ICOLS)=ASTIF(IROWS,ICOLS)-FACTR*ASTIF(IEQNS,ICOLS)
10 CONTINUE
ASLOD(IROWS)=ASLOD(IROWS)-FACTR*ASLOD(IEQNS)
20 CONTINUE
GO TO 50
!
! ADJUST RHS(LOADS) FOR PRESCRIBED DISPLACEMENTS
!
30 DO 40 IROWS=IEQNS,NEQNS
ASLOD(IROWS)=ASLOD(IROWS)-ASTIF(IROWS,IEQNS)*FIXED(IEQNS)
ASTIF(IROWS,IEQNS)=0.0
40 CONTINUE
GO TO 50
60 WRITE(6,900) PIVOT,IEQNS
900 FORMAT(5X,18HINCORRECT PIVOT = ,E20.6,5X,13HEQUATION NO. ,I5)
STOP
50 CONTINUE
RETURN
END
! ****************************************************************
! THIS SUBROUTINE SOLVES THE LOWER TRINGULAR SET OF LINEAR
! EQUATIONS BY USE OF BACK SOLVE SYSTRM.
! ****************************************************************
SUBROUTINE BAKSUB
IMPLICIT REAL*8(A-H,O-Z)
COMMON/FEM1/NPOIN,NELEM,NBOUN,NPROP,NNODE,NEVAB,NSVAB,NDOFN, &
NDIME,NSTRE,NMATS
COMMON/FEM2/PROPS(5,3),COORD(40,2),LNODS(75,2),IFPRE(80), &
FIXED(80),RLOAD(40,3),ELOAD(75,6),MATNO(75), &
STRES(75,6),XDISP(80),TDISP(40,3),TREAC(40,3), &
ASTIF(80,80),ASLOD(80),REACT(80)
!
! THIS ROUTINE PERFORMS THE BACK-
! SUBSTITUTION PHASE
!
NEQNS=NSVAB
DO 5 IEQNS=1,NEQNS
REACT(IEQNS)=0.0
5 CONTINUE
NEQN1=NEQNS+1
DO 30 IEQNS=1,NEQNS
NBACK=NEQN1-IEQNS
PIVOT=ASTIF(NBACK,NBACK)
RESID=ASLOD(NBACK)
IF(NBACK.EQ.NEQNS) GO TO 20
NBAC1=NBACK+1
DO 10 ICOLS=NBAC1,NEQNS
RESID=RESID-ASTIF(NBACK,ICOLS)*XDISP(ICOLS)
10 CONTINUE
20 IF(IFPRE(NBACK).EQ.0) XDISP(NBACK)=RESID/PIVOT
IF(IFPRE(NBACK).EQ.1) XDISP(NBACK)=FIXED(NBACK)
IF(IFPRE(NBACK).EQ.1) REACT(NBACK)=-RESID
30 CONTINUE
KOUNT=0
DO 40 IPOIN=1,NPOIN
DO 40 IDOFN=1,NDOFN
KOUNT=KOUNT+1
TDISP(IPOIN,IDOFN)=XDISP(KOUNT)
40 TREAC(IPOIN,IDOFN)=REACT(KOUNT)
RETURN
END
thanks for your help
Your Friend
Mahdi shabani
Student of Isfahan uniersity

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2015-1-11
That's asking too much of me (or I think anybody to do it as a free favor to you), so if you don't want to do it yourself , then The Mathworks can do it for you. See http://www.mathworks.com/services/consulting/
  3 个评论
显示 1更早的评论
Marc
Marc 2016-7-11
I can probably do it for $99.99 an hour..... Just depends on how fast you REALLY need it.

请先登录，再进行评论。

更多回答（1 个）

Steven Lord
Steven Lord 2016-7-11
Consider creating a FORTRAN MEX-file. Search the documentation for "FORTRAN MEX" and you will find instructions.

类别

在 Help Center 和 File Exchange 中查找有关 Fortran with MATLAB 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by