! CODE1 from "Finite Elements, An Introduction, Vol. 1" ! by E. B. Becker, G. F. Carey, J. T. Oden, Prentice Hall, 1981 ! ! Fortran 90 version converted from F77, at Prof. Becker's request, ! by J. E. Akin, Rice University, Oct. 2, 1998, (akin@rice.edu) ! ! Release 2.0, Allocatable Double precision version MODULE PRECISION_MODULE INTEGER, PARAMETER :: SP = KIND (1.0) ! Single precision INTEGER, PARAMETER :: DP = KIND (1.d0) ! Double precision CONTAINS SUBROUTINE CHECK_HARDWARE_PRECISION ! Verify that the hardware supports the request IF ( KIND (SP) < 0 ) STOP 'ERROR: SP PRECISION NOT AVAILABLE' IF ( KIND (DP) < 0 ) STOP 'ERROR: DP PRECISION NOT AVAILABLE' END SUBROUTINE CHECK_HARDWARE_PRECISION END MODULE PRECISION_MODULE Module Control Use Precision_Module Integer, save :: N_NODE, N_ELEM, N_MAT, N_POINT CHARACTER(len=80) :: TITLE End Module Control Module bc_data Use Control INTEGER, PARAMETER :: MAX_BC = 2, MAX_OPT = 2 INTEGER, ALLOCATABLE :: KPOINT (:) REAL(DP), ALLOCATABLE :: POINT (:) INTEGER :: KBC (MAX_BC) REAL(DP) :: VBC (MAX_OPT, MAX_BC) contains SUBROUTINE ALLOCATE_BC_DATA ALLOCATE ( KPOINT (N_POINT) ) ; ALLOCATE ( POINT (N_POINT) ) END SUBROUTINE ALLOCATE_BC_DATA SUBROUTINE DEALLOCATE_BC_DATA DEALLOCATE ( POINT, KPOINT ) ; END SUBROUTINE DEALLOCATE_BC_DATA End Module bc_data Module elements Use Control INTEGER, PARAMETER :: NODES_PER_EL = 4 INTEGER, ALLOCATABLE :: EL_KIND (:), MAT (:), & NODES (:, :), NINT (:) contains SUBROUTINE ALLOCATE_ELEMENT_ARRAYS IF ( N_ELEM < 1 ) STOP 'INVALID N_ELEM' ALLOCATE ( EL_KIND (N_ELEM), MAT (N_ELEM) ) ALLOCATE ( NODES (NODES_PER_EL, N_ELEM), NINT (N_ELEM) ) ! INITIALIZE GLOBAL ARRAYS (may not be reqd for all systems) EL_KIND = 0 ; MAT = 0 ; NINT = 0 ; NODES = 0 END SUBROUTINE ALLOCATE_ELEMENT_ARRAYS SUBROUTINE DEALLOCATE_ELEMENT_ARRAYS ! in reverse order DEALLOCATE ( NINT, NODES, MAT, EL_KIND ) END SUBROUTINE DEALLOCATE_ELEMENT_ARRAYS End Module elements Module integration Use Precision_Module INTEGER, PARAMETER :: MAX_ORDER = 4 REAL(DP) :: XI (MAX_ORDER, MAX_ORDER), W (MAX_ORDER, MAX_ORDER) End Module integration Module materials Use Control INTEGER, PARAMETER :: MAX_PROP = 4 REAL(DP), ALLOCATABLE :: PROP (:, :) contains SUBROUTINE ALLOCATE_MATERIALS ALLOCATE ( PROP (MAX_PROP, N_MAT) ) ; END SUBROUTINE SUBROUTINE DEALLOCATE_MATERIALS DEALLOCATE ( PROP ) ; END SUBROUTINE DEALLOCATE_MATERIALS End Module materials Module matrix Use Control REAL(DP), ALLOCATABLE :: GK (:, :), GF (:) contains SUBROUTINE ALLOCATE_MATRICES ALLOCATE ( GK (N_NODE, N_NODE), GF (N_NODE) ) ; END SUBROUTINE SUBROUTINE DEALLOCATE_MATRICES DEALLOCATE ( GF, GK ) ; END SUBROUTINE DEALLOCATE_MATRICES End Module matrix Module nodes_data Use Control REAL(DP), ALLOCATABLE :: X (:), U (:) contains SUBROUTINE ALLOCATE_NODES IF ( N_NODE < 2 ) STOP 'INVALID N_NODE' ALLOCATE ( X (N_NODE), U (N_NODE) ) ; END SUBROUTINE SUBROUTINE DEALLOCATE_NODES DEALLOCATE ( U, X ) ; END SUBROUTINE DEALLOCATE_NODES End Module nodes_data PROGRAM CODE1 Use bc_data ; Use elements ; Use materials Use nodes_data ; Use matrix OPEN (5, FILE = 'CODE1.dat') OPEN (6, FILE = 'Code1.out', STATUS = 'UNKNOWN') CALL PRE_PROCESS CALL PROCESS CALL POST_PROCESS CALL DEALLOCATE_MATRICES ; CALL DEALLOCATE_NODES CALL DEALLOCATE_MATERIALS ; CALL DEALLOCATE_ELEMENT_ARRAYS CALL DEALLOCATE_BC_DATA CLOSE (5) ; CLOSE (6) END PROGRAM CODE1 SUBROUTINE PRE_PROCESS Use control ; Use bc_data ; Use elements Use materials ; Use nodes_data ; Use matrix Use integration CALL CHECK_HARDWARE_PRECISION READ (5, * ) TITLE ; WRITE (6, * ) TITLE CALL READ_CONTROLS CALL ALLOCATE_BC_DATA ; CALL ALLOCATE_ELEMENT_ARRAYS CALL ALLOCATE_MATERIALS ; CALL ALLOCATE_NODES CALL ALLOCATE_MATRICES CALL READ_NODES CALL READ_ELEMENTS CALL READ_MATERIALS CALL READ_BOUNDARY_CONDITIONS CALL SET_INTEGRATION_DATA END SUBROUTINE PRE_PROCESS SUBROUTINE READ_CONTROLS Use control READ (5, * ) N_NODE, N_ELEM, N_MAT, N_POINT WRITE (6, 101) N_NODE, N_ELEM, N_MAT, N_POINT 101 FORMAT(/,' N_NODE =',I3/,' N_ELEM =',I3/ & ' N_MAT =',I3/ ' N_POINT=',I3 /) END SUBROUTINE READ_CONTROLS SUBROUTINE READ_NODES Use control ; Use nodes_data IMPLICIT NONE INTEGER :: NREC, NR, N1, N2, DN, I REAL(DP) :: X1, X2, DX, XX REAL(DP), PARAMETER :: LARGE = HUGE(1.d0)/2 X = LARGE ! initalize to invalid number READ (5, * ) NREC DO NR = 1, NREC READ (5, * ) N1, N2, X1, X2 IF ( N1 < 1 .OR. N2 > N_NODE ) & STOP 'ERROR IN NODAL POINT NUMBER DATA' IF ( N2 == 0 ) THEN ! single point X (N1) = X1 ELSE ! create points DN = N2 - N1 ! number of spaces DX = (X2 - X1) / DN ! size of space XX = X1 - DX ! initialize DO I = N1, N2 XX = XX + DX ! increment X (I) = XX END DO END IF END DO WRITE (6, 101) (I, X (I), I = 1, N_NODE) 101 FORMAT(' NODE NO.',I6,5X,'X-COORDINATE',F11.3) DO I = 1, N_NODE IF ( X (I) >= LARGE ) STOP & 'ERROR IN NODAL POINT COORDINATE DATA' END DO END SUBROUTINE READ_NODES SUBROUTINE READ_ELEMENTS Use control ; Use elements IMPLICIT NONE INTEGER :: NREC, NR, NEL1, EL_KIND1, MAT1, NINT1, NODE1 INTEGER :: NUMBER, NN, NEL2, NOD, NEL, N, I REAL(DP) :: X1, X2, DX, XX READ (5, * ) NREC NOD = 1 DO NR = 1, NREC READ (5, * ) NEL1, EL_KIND1, MAT1, NINT1, NODE1, NUMBER NN = EL_KIND1 + 1 NEL2 = NEL1 + NUMBER - 1 NOD = NODE1 DO NEL = NEL1, NEL2 NOD = NOD-1 EL_KIND (NEL) = EL_KIND1 MAT (NEL) = MAT1 NINT (NEL) = NINT1 DO N = 1, NN NOD = NOD+1 NODES (N, NEL) = NOD END DO END DO END DO WRITE (6, "(//,'EL NO EL_KIND MAT NINT NODES')") DO NEL = 1, N_ELEM WRITE (6, "(8I6)") NEL, EL_KIND (NEL), MAT (NEL), NINT (NEL), & (NODES (I, NEL), I = 1, NODES_PER_EL) END DO END SUBROUTINE READ_ELEMENTS SUBROUTINE READ_MATERIALS Use control ; Use materials IMPLICIT NONE INTEGER :: I, J READ (5, * ) ( (PROP (I, J), I = 1, MAX_PROP), J = 1, N_MAT) WRITE (6, "(//,'MAT NO',10X,'K',10X,'C',10X,'B',10X,'F')") WRITE (6, "(I7,(4(1PE11.3)))") (J, (PROP (I, J), & I = 1, MAX_PROP), J = 1, N_MAT) END SUBROUTINE READ_MATERIALS SUBROUTINE READ_BOUNDARY_CONDITIONS Use control ; Use bc_data IMPLICIT NONE INTEGER :: I READ (5, * ) (KBC (I), VBC (1, I), VBC (2, I), I = 1, MAX_BC) WRITE (6, 101) (KBC (I), VBC (1, I), VBC (2, I), I = 1, MAX_BC) 101 FORMAT(//,' BOUNDARY CONDITION DATA', & /,10X,'TYPE',7X,'V1',10X,'V2', & /,' AT X=0:',I5,2(1PE12.3), & /,' AT X=L:',I5,2(1PE12.3)) IF ( N_POINT == 0 ) RETURN READ (5, * ) (KPOINT (I), POINT (I), I = 1, N_POINT) WRITE (6, 103) (KPOINT (I), POINT (I), I = 1, N_POINT) 103 FORMAT(//,' POINT LOAD DATA:',/,(I5,1PE12.3)) END SUBROUTINE READ_BOUNDARY_CONDITIONS SUBROUTINE PROCESS Use control CALL SET_INTEGRATION_DATA CALL FORM_K_AND_F CALL APPLY_BOUNDARY_CONDITIONS CALL SOLVE END SUBROUTINE PROCESS SUBROUTINE SET_INTEGRATION_DATA Use control ; Use integration !.....GAUSS QUADRATURE ORDER 1 XI (1, 1) = 0.d0 W (1, 1) = 2.d0 !.....GAUSS QUADRATURE ORDER 2 XI (1, 2) = - 1.d0 / SQRT (3.d0) XI (2, 2) = - XI (1, 2) W (1, 2) = 1.0d0 W (2, 2) = W (1, 2) !....GAUSS QUADRATURE ORDER 3 XI (1, 3) = - SQRT (3.d0 / 5.d0) XI (2, 3) = 0.d0 XI (3, 3) = - XI (1, 3) W (1, 3) = 5.d0 / 9.d0 W (2, 3) = 8.d0 / 9.d0 W (3, 3) = W (1, 3) !.....GAUSS QUADRATURE ORDER 4 XI (1, 4) = - 0.8611363115d0 XI (2, 4) = - 0.3399810435d0 XI (3, 4) = - XI (2, 4) XI (4, 4) = - XI (1, 4) W (1, 4) = 0.3478548451d0 W (2, 4) = 0.6521451549d0 W (3, 4) = W (2, 4) W (4, 4) = W (1, 4) END SUBROUTINE SET_INTEGRATION_DATA SUBROUTINE FORM_K_AND_F Use control ; Use elements ; Use matrix Use integration ; Use nodes_data IMPLICIT NONE INTEGER :: NG, NEL, N, I1, I2, I3 REAL(DP) :: EK (NODES_PER_EL, NODES_PER_EL), EF (NODES_PER_EL) GK = 0.d0 ! zero global stiffness GF = 0.d0 ! zero global source DO NEL = 1, N_ELEM N = EL_KIND (NEL) + 1 I1 = NODES (1, NEL) I2 = NODES (N, NEL) I3 = NINT (NEL) CALL ELEM (X (I1), X (I2), N, EK, EF, I3, XI (1, I3), & W (1, I3), MAT (NEL) ) CALL ASSEMBLE (EK, EF, N, NODES (1, NEL), GK, GF, N_NODE) END DO END SUBROUTINE FORM_K_AND_F SUBROUTINE ELEM (X1, X2, N, EK, EF, NL, XI, W, MAT) Use control IMPLICIT NONE INTEGER, INTENT(IN) :: N, NL, MAT REAL(DP), INTENT(IN) :: X1, X2 REAL(DP), INTENT(INOUT) :: EK (N, N), EF (N), XI (NL), W (NL) REAL(DP) :: PSI (N), D_PSI (N), DX, X, XK, XC, XB, XF INTEGER :: L, I, J DX = (X2 - X1) / 2. ! jacobian !.....INITIALIZE ELEMENT ARRAYS EK = 0.d0 ! zero element stiffness EF = 0.d0 ! zero element source !.....BEGIN INTEGRATION POINT LOOP DO L = 1, NL ! over quafrature points X = X1 + (1. + XI (L)) * DX ! integration coordinate CALL GET_MAT (MAT, X, XK, XC, XB, XF) CALL SHAPE (XI (L), N, PSI, D_PSI) DO I = 1, N ! over element rows EF (I) = EF (I) + PSI (I) * XF * W (L) * DX DO J = 1, N ! over element columns EK (I, J) = EK (I, J) & + (XK * D_PSI (I) * D_PSI (J) / DX / DX & + XC * PSI (I) * D_PSI (J) / DX & + XB * PSI (I) * PSI (J) ) * W (L) * DX END DO ! over element columns END DO ! over element rows END DO ! over quafrature points END SUBROUTINE ELEM SUBROUTINE SHAPE (XI, N, PSI, D_PSI) Use control IMPLICIT NONE INTEGER, INTENT(IN) :: N REAL(DP), INTENT(IN) :: XI REAL(DP), INTENT(OUT) :: PSI (N), D_PSI (N) SELECT CASE (N) ! the number of nodes CASE (2) !.....LINEAR SHAPE FUNCTIONS PSI (1) = 0.5d0 * (1. - XI) PSI (2) = 0.5d0 * (1. + XI) D_PSI (1) = - 0.5d0 D_PSI (2) = 0.5d0 CASE (3) !.....QUADRATIC ELEMENTS PSI (1) = XI * (XI - 1.d0) * 0.5d0 PSI (2) = 1.d0 - XI**2 PSI (3) = XI * (XI + 1.d0) * 0.5d0 D_PSI (1) = XI - 0.5d0 D_PSI (2) = - 2.d0 * XI D_PSI (3) = XI + 0.5d0 CASE (4) !.....CUBIC ELEMENTS PSI (1) = 9.d0 / 2.d0 * (1.d0 / 9.d0 - XI**2) * (1.d0 - XI) PSI (2) = 16.d0 / 27.d0 * (1.d0 - XI**2) * (1.d0 / 3.d0 - XI) PSI (3) = 16.d0 / 27.d0 * (1.d0 - XI**2) * (1.d0 / 3.d0 + XI) PSI (4) = 9.d0 / 2.d0 * (1.d0 / 9.d0 - XI**2) * (1.d0 + XI) D_PSI (1) = 9.d0 / 2.d0 * (3.d0 * XI**2 & - 2.d0 * XI - 1.d0 / 9.d0) D_PSI (2) = 16.d0 / 27.d0 * (3.d0 * XI**2 & - 2.d0 / 3.d0 * XI - 1.d0) D_PSI (3) = 16.d0 / 27.d0 * ( - 3.d0 * XI**2 & - 2.d0 / 3.d0 * XI + 1.d0) D_PSI (4) = 9.d0 / 2.d0 * ( - 3.d0 * XI**2 - 2.d0 * XI & + 1.d0 / 9.d0) CASE DEFAULT WRITE (6, "(' ERROR IN CALL TO SHAPE, N=',I3)") STOP ' ERROR IN CALL TO SHAPE' END SELECT END SUBROUTINE SHAPE SUBROUTINE GET_MAT (MAT, X, XK, XC, XB, XF) Use control ; Use materials IMPLICIT NONE INTEGER, INTENT(IN) :: MAT REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: XK, XC, XB, XF XK = PROP (1, MAT) XC = PROP (2, MAT) XB = PROP (3, MAT) XF = X END SUBROUTINE GET_MAT SUBROUTINE ASSEMBLE (EK, EF, N, NODE, GK, GF, NG) Use control IMPLICIT NONE INTEGER, INTENT(IN) :: N, NG, NODE (N) REAL(DP), INTENT(IN) :: EK (N, N), EF (N) REAL(DP), INTENT(INOUT) :: GK (NG, NG), GF (NG) INTEGER :: I, IG, J, JG DO I = 1, N IG = NODE (I) GF (IG) = GF (IG) + EF (I) DO J = 1, N JG = NODE (J) GK (IG, JG) = GK (IG, JG) + EK (I, J) END DO END DO END SUBROUTINE ASSEMBLE SUBROUTINE APPLY_BOUNDARY_CONDITIONS Use control ; Use matrix ; Use bc_data IMPLICIT NONE !.....BOUNDARY CONDITION AT NODE 1 IF ( KBC (1) < 1 .OR. KBC (1) > 3 ) GO TO 99 GO TO (10, 11, 12) KBC (1) 10 CALL DIRCHLET (VBC (1, 1), 1, N_NODE) GO TO 15 11 GF (1) = GF (1) + VBC (1, 1) GO TO 15 12 GK (1, 1) = GF (1) + VBC (1, 1) GF (1) = GF (1) + VBC (1, 1) * VBC (2, 1) !.....BOUNDARY CONDITION AT NODE N_NODE 15 IF ( KBC (2) < 1 .OR. KBC (2) > 3 ) GO TO 99 GO TO (20, 21, 22) KBC (2) 20 CALL DIRCHLET (VBC (1, 2), N_NODE, N_NODE) RETURN 21 GF (N_NODE) = GF (N_NODE) + VBC (1, 2) RETURN 22 GK (N_NODE, N_NODE) = GK (N_NODE, N_NODE) + VBC (1, 2) GF (N_NODE) = GF (N_NODE) + VBC (1, 2) * VBC (2, 2) RETURN 99 WRITE (6, "('ERROR IN APPLY_BOUNDARY_CONDITIONS KBC= ',I3)") KBC STOP 'ERROR IN APPLY_BOUNDARY_CONDITIONS' END SUBROUTINE APPLY_BOUNDARY_CONDITIONS SUBROUTINE DIRCHLET (VAL, N, NG) Use control ; Use matrix IMPLICIT NONE INTEGER, INTENT(IN) :: N, NG REAL(DP), INTENT(IN) :: VAL INTEGER :: I DO I = 1, NG GF (I) = GF (I) - GK (I, N) * VAL GK (I, N) = 0.d0 GK (N, I) = 0.d0 END DO GK (N, N) = 1.d0 GF (N) = VAL END SUBROUTINE DIRCHLET SUBROUTINE SOLVE Use control ; Use matrix ; Use nodes_data IMPLICIT NONE CALL TRI (GK, N_NODE, N_NODE) CALL BACK (GK, U, GF, N_NODE, N_NODE) END SUBROUTINE SOLVE SUBROUTINE TRI (A, N, M) Use control IMPLICIT NONE INTEGER, INTENT(IN) :: N, M REAL(DP), INTENT(INOUT) :: A (N, N) INTEGER :: I, J, K REAL(DP) :: FAC REAL(DP), PARAMETER :: NIL = 2 * TINY (1.d0) DO I = 1, M - 1 IF ( ABS (A (I, I) ) < NIL ) THEN ! nearly singular WRITE (6, 100) I, A (I, I) 100 FORMAT('ELIMINATION FAILED DUE TO SMALL PIVOT',/ & ' EQUATION NO.',I5,' PIVOT',1PE10.3) STOP 'ELIMINATION FAILED DUE TO SMALL PIVOT' END IF DO J = I + 1, M IF ( A (J, I) == 0.d0 ) CYCLE ! to next J FAC = A (J, I) / A (I, I) DO K = I + 1, M A (J, K) = A (J, K) - A (I, K) * FAC END DO END DO END DO END SUBROUTINE TRI SUBROUTINE BACK (A, X, B, N, M) Use control IMPLICIT NONE INTEGER, INTENT(IN) :: N, M REAL(DP), INTENT(INOUT) :: A (N, N), X (M), B (M) INTEGER :: M1, I, J1, J, IB M1 = M - 1 DO I = 1, M1 J1 = I + 1 DO J = J1, M B (J) = B (J) - B (I) * A (J, I) / A (I, I) END DO END DO X (M) = B (M) / A (M, M) DO I = 1, M1 IB = M - I J1 = IB + 1 DO J = J1, M B (IB) = B (IB) - A (IB, J) * X (J) END DO X (IB) = B (IB) / A (IB, IB) END DO END SUBROUTINE BACK SUBROUTINE POST_PROCESS Use control IMPLICIT NONE INTEGER :: NEL1, NEL2, I, NEL, NPTS, NPT REAL(DP) :: OUTPTS (10), XX, UH, UX, SIGH, SIGX, DU, DSIG CHARACTER (len = 5) :: CMMND !ea 1 READ (5, 100) CMMND 100 FORMAT(A5,2I5) WRITE (6, * ) cmmnd IF ( CMMND == ' SOLN' ) THEN !.... OUTPUT SOLUTION AT SELECTED POINTS READ (5, * ) NEL1, NEL2 IF ( NEL1 + NEL2 == 0 ) THEN NEL1 = 1 NEL2 = N_ELEM END IF READ (5, * ) NPTS, (OUTPTS (I), I = 1, NPTS) WRITE (6, 102) NPTS, (OUTPTS (I), I = 1, NPTS) 102 FORMAT(//,'EVALUATE SOLUTION AT',I3, & ' OUTPUT POINTS PER ELEMENT',/,' XI= ',(10(F6.3,','))) WRITE (6, '(//,A,//)') TITLE WRITE (6, 103) 103 FORMAT(4X, 'X', 8X, 'U', 10X, 'U', 10X, 'U', 7X, & 'K DU/DX', 4X, 'K DU/DX', 4X, 'K DU/DX',/, 12X, 'F E', 7X, & 'EXACT', 6X, 'ERROR', 7X, 'F E', 7X, 'EXACT', 6X, 'ERROR') DO NEL = NEL1, NEL2 WRITE (6, "('ELEMENT NO.',I3)") NEL DO NPT = 1, NPTS CALL EVAL (NEL, OUTPTS (NPT), XX, UH, UX, SIGH, SIGX) DU = UX - UH DSIG = SIGX - SIGH WRITE (6, "(F8.4,6(1PE11.3))") XX, UH, UX, DU, & SIGH, SIGX, DSIG END DO END DO GO TO 1 ELSE IF ( CMMND == ' NORM' ) THEN !..... CALL ERROR NORM ROUTINE (IF EXACT SOLUTION IS AVAILABLE) 20 CALL ERROR_NORMS GO TO 1 END IF TITLE (1:5) = CMMND END SUBROUTINE POST_PROCESS SUBROUTINE EVAL (NEL, XI, XX, UH, UX, SIGH, SIGX) Use control ; Use nodes_data ; Use elements IMPLICIT NONE INTEGER, INTENT(IN) :: NEL REAL(DP), INTENT(IN) :: XI REAL(DP), INTENT(OUT) :: XX, UH, UX, SIGH, SIGX INTEGER :: N, I, I1 REAL(DP) :: X1, X2, DX, XK, XC, XB, XF REAL(DP) :: PSI (NODES_PER_EL), D_PSI (NODES_PER_EL) N = EL_KIND (NEL) + 1 CALL SHAPE (XI, N, PSI, D_PSI) !.... CALCULATE UH AND DUH/DX FROM SHAPE FUNCTIONS UH = 0. SIGH = 0. DO I = 1, N I1 = NODES (I, NEL) IF ( I == 1 ) X1 = X (I1) IF ( I == N ) X2 = X (I1) UH = UH + PSI (I) * U (I1) SIGH = SIGH + D_PSI (I) * U (I1) END DO DX = (X2 - X1) * 0.5d0 XX = X1 + (1.d0 + XI) * DX CALL GET_MAT (MAT (NEL), XX, XK, XC, XB, XF) SIGH = SIGH * XK / DX CALL EXACT (XX, UX, SIGX) END SUBROUTINE EVAL SUBROUTINE EXACT (X, U, SIG) Use control IMPLICIT NONE REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: U, SIG REAL(DP) :: E, EX, SINH_1, SINH_X, COSH_X !.... Note: COSH, SINH are in F90, but F77 style kept here !.... MODEL PROBLEM: -D2U/DX2 + U = X !.... U(X) = X - SINH(X)/SINH(1) !.... SIG(X) = 1 - COSH(X)/SINH(1) E = EXP (1.d0) EX = EXP (X) SINH_1 = E - 1.d0 / E SINH_X = EX - 1.d0 / EX COSH_X = EX + 1.d0 / EX U = X - SINH_X / SINH_1 SIG = 1.d0 - COSH_X / SINH_1 END SUBROUTINE EXACT SUBROUTINE ERROR_NORMS Use control ; Use nodes_data ; Use elements IMPLICIT NONE INTEGER :: I, N1, N, N2, K REAL(DP) :: SQ_NORM, E_NORM, WT, XI, XX, UH, UX, SIGH, SIGX REAL(DP) :: XK, XC, XB, XF, W INTEGER, PARAMETER :: TRAPEZOID = 50 ! quadrature intervals REAL(DP), PARAMETER :: DXI = 2.0d0 / TRAPEZOID ! spacings SQ_NORM = 0.d0 E_NORM = 0.d0 ! DXI = 1.d0 / 25.d0 !ea ?? DO I = 1, N_ELEM N1 = NODES (1, I) N = EL_KIND (I) + 1 N2 = NODES (N, I) ! WT = (X (N2) - X (N1) ) / 50.d0 WT = (X (N2) - X (N1) ) / TRAPEZOID XI = - 1.d0 DO K = 1, TRAPEZOID + 1 XX = X (N1) + (1.d0 + XI) * (X (N2) - X (N1) ) / 2.d0 CALL EVAL (I, XI, XX, UH, UX, SIGH, SIGX) CALL GET_MAT (MAT (I), XX, XK, XC, XB, XF) IF ( K == 1 .OR. K == TRAPEZOID + 1 ) THEN W = WT / 2.d0 ELSE W = WT END IF ! XI = XI + 0.02d0 !ea error ?? xxxxxxxxxxxxx XI = XI + DXI SQ_NORM = SQ_NORM + (UX - UH) **2 * W E_NORM = E_NORM + (SIGX - SIGH) **2 * W / XK & + (UX - UH) **2 * W * XB END DO END DO SQ_NORM = SQRT (SQ_NORM) E_NORM = SQRT (E_NORM) WRITE (6, "(//,A,/)") TITLE WRITE (6, 100) SQ_NORM, E_NORM 100 FORMAT('THE H0 NORM OF THE ERROR IS ',1PE13.3/ & 'THE ENERGY NORM OF THE ERROR IS',1PE12.3) END SUBROUTINE ERROR_NORMS