Demo entry 6633110

123

   

Submitted by anonymous on Jul 31, 2017 at 16:56
Language: Fortran. Code size: 7.6 kB.

    PROGRAM FEM
!
!   === VARIABLE DECLARATION ===
!
    INTEGER NE,NP,NP2
    INTEGER ELEM,NODE
    REAL X(1000),Y(1000)
    INTEGER NBCX(1000),NBCY(1000),NBC(1000),NDBC
    INTEGER IA(4,1000)
    REAL TT(1000),EE(1000),POI(1000)
    INTEGER NL,NLOAD
    REAL LOAD,TRF(1000)
!
    REAL D(3,3)
    REAL XI(4),ETA(4)
    REAL XX(4),YY(4),JAC(2,2)
    REAL DETJ,INVJ(2,2)
    REAL W(4,4,4),DNDXI(4),DNDETA(4)
    REAL BB(3,2,4),B(3,8)
    REAL KK(8,8),KE(8,8),KG(1000,1000)
    INTEGER ORD(8)
!
    REAL AA,U(1000)
!
!   === READ DATA ===
!
    OPEN (10,FILE='data.dat')
!
!       +++ INITIAL COODINATES & BOUNDARY CONDITION +++
!
        READ (10,'(2I5)') NE,NP
        NP2=2*NP
        DO 100 I=1,NP
            READ (10,'(I5,2F10.3,2I5)') NODE,X(I),Y(I),NBCX(I),NBCY(I)
            NBC(2*I-1)=NBCX(I)
            NBC(2*I)=NBCY(I)
100     CONTINUE
        NDBC=0
        DO 110 I=1,NP2
            IF (NBC(I) .EQ. 0) THEN
                NDBC=NDBC+1
            END IF
110     CONTINUE
!
!       +++ ELEMENT CONDITION +++
!
        DO 120 I=1,NE
            READ (10,'(5I5,F10.3,F10.1,F5.2)') ELEM,(IA(J,I),J=1,4),TT(I),EE(I),POI(I)
120     CONTINUE
!
!       +++ LOAD DATA +++
!
        READ (10,'(I5)') NL
        DO 130 I=1,NL
            TRF(I)=0
130     CONTINUE
        DO 140 I=1,NL
            READ (10,'(I5,F20.5)') NLOAD,LOAD
            TRF(NLOAD)=LOAD
140     CONTINUE
!
    CLOSE (10)
!
!   === CREATE STIFFNESS MATRIX ===
!
!   +++ GAUSSIAN POINT +++
!
    XI(1)=1.0D0/SQRT(3.0D0)
    XI(2)=1.0D0/SQRT(3.0D0)
    XI(3)=-1.0D0/SQRT(3.0D0)
    XI(4)=-1.0D0/SQRT(3.0D0)
!
    ETA(1)=1.0D0/SQRT(3.0D0)
    ETA(2)=-1.0D0/SQRT(3.0D0)
    ETA(3)=1.0D0/SQRT(3.0D0)
    ETA(4)=-1.0D0/SQRT(3.0D0)
!
!   +++ ELEMENT STIFFNESS MATRIX +++
!
    DO 700 I=1,NP2
        DO 710 J=1,NP2
            KG(I,J)=0.0D0
710     CONTINUE
700 CONTINUE
!
    DO 200 N=1,NE
!
!   +++ D MATRIX +++
!
        D(1,1)=1.0D0
        D(1,2)=POI(N)
        D(1,3)=0.0D0
        D(2,1)=POI(N)
        D(2,2)=1.0D0
        D(2,3)=0.0D0
        D(3,1)=0.0D0
        D(3,2)=0.0D0
        D(3,3)=(1.0D0-POI(N))/2.0D0
!
        DO 300 I=1,3
            DO 310 J=1,3
                D(I,J)=D(I,J)*EE(N)/(1.0D0-POI(N)**2.0D0)
310         CONTINUE
300     CONTINUE
!
        DO 210 M=1,4
            XX(M)=X(IA(M,N))
            YY(M)=Y(IA(M,N))
210     CONTINUE
!
        DO 400 M=1,4
!
!       *** J MATRIX ***
!
            JAC(1,1)=1.0D0/4.0D0*((XX(2)-XX(1))*(1-ETA(M))+(XX(3)-XX(4))*(1+ETA(M)))
            JAC(1,2)=1.0D0/4.0D0*((YY(2)-YY(1))*(1-ETA(M))+(YY(3)-YY(4))*(1+ETA(M)))
            JAC(2,1)=1.0D0/4.0D0*((XX(4)-XX(1))*(1-XI(M))+(XX(3)-XX(2))*(1+XI(M)))
            JAC(2,2)=1.0D0/4.0D0*((YY(4)-YY(1))*(1-XI(M))+(YY(3)-YY(2))*(1+XI(M)))
!
!       *** DETERMINANT AND INVERSE OF J MATRIX ***
!
            DETJ=JAC(1,1)*JAC(2,2)-JAC(1,2)*JAC(2,1)
!
            INVJ(1,1)=JAC(2,2)/DETJ
            INVJ(1,2)=-JAC(1,2)/DETJ
            INVJ(2,1)=-JAC(2,1)/DETJ
            INVJ(2,2)=JAC(1,1)/DETJ
!
!       *** COMPONENTS OF B MATRIX ***
!
            DO 410 I=1,4
                DO 420 J=1,4
                    DO 430 K=1,4
                        W(I,J,K)=0.0D0
430                 CONTINUE
420             CONTINUE
410         CONTINUE
!
            DNDXI(1)=-1.0D0/4.0D0*(1.0D0-ETA(M))
            DNDXI(2)=1.0D0/4.0D0*(1.0D0-ETA(M))
            DNDXI(3)=1.0D0/4.0D0*(1.0D0+ETA(M))
            DNDXI(4)=-1.0D0/4.0D0*(1.0D0+ETA(M))
!
            DNDETA(1)=-1.0D0/4.0D0*(1.0D0-XI(M))
            DNDETA(2)=-1.0D0/4.0D0*(1.0D0+XI(M))
            DNDETA(3)=1.0D0/4.0D0*(1.0D0+XI(M))
            DNDETA(4)=1.0D0/4.0D0*(1.0D0-XI(M))
!
            W(1,1,1)=1.0D0
            W(2,4,1)=1.0D0
            W(3,2,1)=1.0D0
            W(3,3,1)=1.0D0
!
            DO 440 I=1,2
                DO 450 J=1,2
                    W(I,J,2)=INVJ(I,J)
                    W(I+2,J+2,2)=INVJ(I,J)
450             CONTINUE
440         CONTINUE
            DO 460 I=1,3
                DO 470 J=1,4
                    DO 480 K=1,4
                        W(I,J,4)=W(I,J,4)+W(I,K,1)*W(K,J,2)
480                 CONTINUE
470             CONTINUE
460         CONTINUE
!
            DO 490 I=1,3
                DO 500 J=1,2
                    DO 510 K=1,4
                        BB(I,J,K)=0.0D0
510                 CONTINUE
500             CONTINUE
490         CONTINUE
!
            DO 520 L=1,4
                W(1,1,3)=DNDXI(L)
                W(2,1,3)=DNDETA(L)
                W(3,2,3)=DNDXI(L)
                W(4,2,3)=DNDETA(L)
                DO 530 I=1,3
                    DO 540 J=1,2
                        DO 550 K=1,4
                            BB(I,J,L)=BB(I,J,L)+W(I,K,4)*W(K,J,3)
550                     CONTINUE
540                 CONTINUE
530             CONTINUE
520         CONTINUE
!
!       *** B MATRIX ***
!
            DO 560 I=1,3
                DO 570 J=1,8
                    B(I,J)=0.0D0
570             CONTINUE
560         CONTINUE
!
            DO 580 L=1,4
                DO 590 I=1,3
                    DO 600 J=1,2
                        B(I,J+2*L-2)=BB(I,J,L)
600                 CONTINUE
590             CONTINUE
580         CONTINUE
!
!       *** STIFFNESS MATRIX AT GAUSSIAN POINT ***
!
            DO 610 I=1,8
                DO 620 J=1,8
                    KK(I,J)=0.0D0
                    IF (M .EQ. 1) THEN
                        KE(I,J)=0.0D0
                    END IF
620             CONTINUE
610         CONTINUE
!
            DO 630 I=1,8
                DO 640 J=1,8
                    DO 650 K=1,3
                        DO 660 L=1,3
                            KK(I,J)=KK(I,J)+B(K,I)*D(K,L)*B(L,J)
660                     CONTINUE
650                 CONTINUE
                    KE(I,J)=KE(I,J)+TT(N)*KK(I,J)*DETJ
640             CONTINUE
630         CONTINUE
!
400     CONTINUE
!
!   +++ GLOBAL STIFFNESS MATRIX +++
!
        DO 750 I=1,4
            ORD(2*I-1)=2*IA(I,N)-1
            ORD(2*I)=2*IA(I,N)
750     CONTINUE
!
        DO 730 I=1,8
            DO 740 J=1,8
                KG(ORD(I),ORD(J))=KG(ORD(I),ORD(J))+KE(I,J)
740         CONTINUE
730     CONTINUE
!
200 CONTINUE
!
!   === CONSIDER BOUNDARY CONDITION ===
!
    DO 800 I=1,NP2
        IF (NBC(I) .EQ. 0) THEN
            DO 810 J=1,NP2
                KG(I,J)=0.0D0
                KG(J,I)=0.0D0
810         CONTINUE
            KG(I,I)=1.0D0
            TRF(I)=0.0D0
        END IF
800 CONTINUE
!
!   === SOLVE GLOBAL STIFFNESS EQUATION  ===
!
	DO 900 I=1,NP2
		AA=KG(I,I)
		DO 910 J=1,NP2
			KG(I,J)=KG(I,J)/AA
910     CONTINUE
		TRF(I)=TRF(I)/AA
		DO 920 J=1,NP2
			IF(I .NE. J) THEN
				AA=KG(J,I)
				DO 930 K=1,NP2
					KG(J,K)=KG(J,K)-KG(I,K)*AA
930             CONTINUE
				TRF(J)=TRF(J)-TRF(I)*AA
			END IF
920     CONTINUE
900	CONTINUE
!
	DO 950 I=1,NP2
		U(I)=TRF(I)
950 CONTINUE
!
!   === WRITE DATA ===
!
    OPEN (20,FILE='result.dat')
        WRITE (20,'(A25)') '========================='
        WRITE (20,'(A25)') '    i      u(i)      v(i)     '
        WRITE (20,'(A25)') '-------------------------'
        DO 999 I=1,NP
            WRITE (20,'(I5,2F10.3)') I,U(2*I-1),U(2*I) 
999     CONTINUE
        WRITE (20,'(A25)') '========================='
    CLOSE (20)
!
    WRITE (6,*) '           '
    WRITE (6,*) '*** END ***'
    WRITE (6,*) '           '
!
    PAUSE
    STOP
!
    END

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).