! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE BT_D_B_OPTIONS (D, B, S, ROWS, COLS, IOPT, CONST) ! * * * * * * * * * * * * * * * * * * * * * * ! SPECIAL MATRIX MULTIPLICATION OPERATION ! IF IOPT=0, S = (B)T*D*B*COEFF PRODUCT ! IF IOPT=1, S = (B)T*D*B*COEFF + S NUM. INTG. ! * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: ROWS, COLS, IOPT REAL(DP), OPTIONAL :: CONST ! default = 1 REAL(DP), INTENT(IN) :: D (ROWS, ROWS), B (ROWS, COLS) REAL(DP), INTENT(OUT) :: S (COLS, COLS) ! INTEGER :: I, J, K, L ! REAL(DP) :: SUM, DBIK, COEFF REAL(DP) :: COEFF ! D(ROWS,M) = SYMMETRIC SQUARE MATRIX ! B(ROWS,COLS) = RECTANGLUAR ARRAY ! S(COLS,COLS) = RETURNED SYMMETRIC SQ MATRIX ! CONST = SCALAR COEFFICIENT IF ( PRESENT (CONST) ) THEN COEFF = CONST ELSE COEFF = 1.0_DP END IF SELECT CASE (IOPT) CASE (0) S = COEFF * MATMUL ( (MATMUL( TRANSPOSE(B),D)), B) CASE (1) S = S + COEFF * MATMUL ( (MATMUL( TRANSPOSE(B),D)), B) CASE DEFAULT ; STOP 'INVALID OPTION, BT_D_B_OPTIONS' END SELECT !9 use matmul and transpose here !DO L = 1, COLS ! DO K = 1, COLS ! SUM = 0.0_DP ! DO I = 1, ROWS ! DBIK = 0.0_DP ! DO J = 1, ROWS ! ! USE SYMMETRY OF D ! DBIK = DBIK + D (J, I) * B (J, K) ! END DO ! SUM = SUM + B (I, L) * DBIK ! END DO ! IF ( IOPT == 0 ) THEN ! S (L, K) = SUM * COEFF ! ELSE ! S (L, K) = S (L, K) + SUM * COEFF ! END IF ! END DO !END DO END SUBROUTINE BT_D_B_OPTIONS SUBROUTINE BT_DIAG_B_OPTIONS (DIAG, B, S, ROWS, COLS, IOPT, CONST) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! SPECIAL DIAGONIAL MATRIX MULTIPLICATION OPERATION ! IF IOPT=0, S = (B)T*DIAG*B*COEFF PRODUCT ONLY ! IF IOPT=1, S = (B)T*DIAG*B*COEFF + S NUM. INTEG. ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: ROWS, COLS, IOPT REAL(DP), OPTIONAL :: CONST ! default = 1 REAL(DP), INTENT(IN) :: DIAG (ROWS), B (ROWS, COLS) REAL(DP), INTENT(OUT) :: S (COLS, COLS) INTEGER :: I, K, L REAL(DP) :: SUM, COEFF ! DIAG(ROWS) = DIAGONAL MATRIX ! B(ROWS,COLS) = RECTANGLUAR ARRAY ! S(COLS,COLS) = RETURNED SYMMETRIC SQ MATRIX ! CONST = SCALAR COEFFICIENT IF ( PRESENT (CONST) ) THEN COEFF = CONST ELSE COEFF = 1.0_DP END IF !9 use matmul and transpose here DO L = 1, COLS DO K = 1, COLS SUM = 0.0_DP DO I = 1, ROWS SUM = SUM + B (I, L) * DIAG (I) * B (I, K) END DO IF ( IOPT == 0 ) THEN S (L, K) = SUM * COEFF ELSE S (L, K) = S (L, K) + SUM * COEFF END IF END DO END DO END SUBROUTINE BT_DIAG_B_OPTIONS SUBROUTINE CONDENSE_MATRICES (N_TOTAL, N_EL_FRE, S, C) ! * * * * * * * * * * * * * * * * * * * * * * * ! CONDENSATION OF ELEMENT MATRICES TO REMOVE ! INTERNAL DEGREES OF FREEDOM ! * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: N_TOTAL, N_EL_FRE REAL(DP), INTENT(INOUT) :: S (N_TOTAL, N_TOTAL), C (N_TOTAL) REAL(DP), PARAMETER :: ZERO = 0.d0 INTEGER :: I, J, K, L, M, NELIM REAL(DP) :: CK, SKK, SLKSKK ! INTERNAL DEGREES OF FREEDOM *MUST* COME LAST ! : SAA : SAB : : DA : : CA : ! :......:......: :....: = :....: ! : SBA : SBB : : DB : : CB : ! ENTER FULL S ; RETURN CONDENSED IN SAA AND CA ! DIMENSION SAA(N_EL_FRE,N_EL_FRE), CA(N_EL_FRE) ! SAA* = (SAA) - (SAB)*(SBB)I*(SAB)T ! CA* = (CA) - (SAB)*(SBB)I*(CB) ! N_TOTAL = ORIG. NO. OF D.O.F. OF ELEMENT ! N_EL_FRE = FINAL NO. OF D.O.F. OF ELEMENT ! NELIM = NUMBER OF DOF TO ELIMINATE ! S = SQUARE ELEMENT MATRIX ! C = ELEMENT COLUMN MATRIX NELIM = N_TOTAL - N_EL_FRE DO I = 1, NELIM J = N_TOTAL - I K = J + 1 SKK = S (K, K) CK = C (K) IF ( SKK /= ZERO ) THEN C (K) = CK / SKK DO L = 1, J SLKSKK = S (L, K) / SKK S (L, K) = SLKSKK DO M = L, J S (L, M) = S (L, M) - S (K, M) * SLKSKK S (M, L) = S (L, M) END DO C (L) = C (L) - CK * SLKSKK END DO END IF END DO END SUBROUTINE CONDENSE_MATRICES SUBROUTINE CONTOUR_ISOPAR_SURF (L_TYPE, MAX_PTS, PT_INOUT, V, COORD, & XYZ_PTS, KOUNT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET COORDINATES OF POINTS ON A CONTOUR ON AN ISOPARAMETRIC SURFACE ! Int. J. Num. Meth. Engr., v. 14, n. 3, 1979 ! * * * * * * ** * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP, NOD_PER_EL, N_SPACE IMPLICIT NONE INTEGER, INTENT(IN) :: L_TYPE, MAX_PTS REAL(DP), INTENT(IN) :: COORD (NOD_PER_EL, N_SPACE), V (NOD_PER_EL) REAL(DP), INTENT(INOUT) :: PT_INOUT (2) INTEGER, INTENT(OUT) :: KOUNT REAL(DP), INTENT(OUT) :: XYZ_PTS (N_SPACE, MAX_PTS) REAL(DP) :: H_L (NOD_PER_EL), DH_L (2, NOD_PER_EL), DER (2) REAL(DP) :: PT_L (2, MAX_PTS), R_NEW, S_NEW REAL(DP) :: TRY_DL, DL, DV_DR, DV_DS, GRAD INTEGER :: I, J LOGICAL :: OUTSIDE EQUIVALENCE (DER (1), DV_DR), (DER (2), DV_DS) ! COORD = GLOBAL COORD. OF NODES OF ELEMENT ! DH_L = LOCAL COORD. ELEM_DERIVATIVES OF H_L ! DL = LOCAL COORD. LENGTH OF EACH SEGMENT ! H_L = ELEMENT INTERPOLATION FUNCTIONS ! KOUNT = NO. OF PTS. ON CONTOUR CURVE ! L_TYPE = ELEM. TYPE, 0=QUADRILATERAL, 1=TRIANGLE ! MAX_PTS = MAX. NO. CONTOUR POINTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_SPACE = DIMENSION OF GLOBAL SPACE, 2 OR 3 ! PT_INOUT = LOCAL COORD OF 1ST (IN) OR LAST (OUT) POINT. ! USUALLY THE INPUT POINT IS WHERE THE CONTOUR CROSSES ! A POINT ON THE ELEMENT BOUNDARY. ! PT_L = LOCAL COORDINATES OF THE CONTOUR PATH ! V = NODAL VALUES OF QUANTITY TO BE CONTOURED ! XYZ_PTS = GLOBAL COORD. ARRAY FOR CONTOUR PTS. ! INITIALIZE KOUNT = 1 ; TRY_DL = 2.5d0 / MAX_PTS DL = TRY_DL ; OUTSIDE = .FALSE. IF ( L_TYPE == 0 ) DL = DL * 2.d0 PT_L (1:2, 1) = PT_INOUT (1:2) ! MARCH ALONG CONTOUR (ALLOWING FOR ONE DIRECTION REVERSAL) DO J = 1, MAX_PTS + 2 ! FORM SHAPE FUNCTION LOCAL ELEM_DERIVATIVES CALL SCALAR_DERIVS (PT_L (:, KOUNT), DH_L) ! FORM LOCAL ELEM_DERIVATIVE OF VARIABLE, DER = DH*V DER = MATMUL ( DH_L, V ) GRAD = SQRT ( DOT_PRODUCT (DER, DER) ) ! FIND LOCAL COORD. OF SEGMENT END R_NEW = PT_L (1, KOUNT) - DL * DV_DS / GRAD S_NEW = PT_L (2, KOUNT) + DL * DV_DR / GRAD ! IS NEXT POINT IN THE ELEMENT IF (L_TYPE == 0) THEN ! QUADRILATERAL (-1 TO +1): IS POINT OUTSIDE ? IF ( ABS (R_NEW) > 1.d0 .OR. & ABS (S_NEW) > 1.d0 ) OUTSIDE = .TRUE. ELSE ! TRIANGLE (UNIT COORDINATES): OUTSIDE ? IF ( R_NEW < 0.d0 .OR. S_NEW < 0.d0 .OR. & (R_NEW + S_NEW) > 1.d0) OUTSIDE = .TRUE. ENDIF IF ( .NOT. OUTSIDE ) THEN ! ADD POINT TO CONTOUR LIST KOUNT = KOUNT + 1 PT_L (1, KOUNT) = R_NEW ; PT_L (2, KOUNT) = S_NEW IF ( KOUNT == MAX_PTS ) EXIT ! THE LOOP ELSE ! DID IT GO OUT ON FIRST STEP ? IF ( KOUNT == 1 ) THEN DL = -1.1d0 * DL ! REVERSE DIRECTION ! IS THE CONTOUR A SINGLE POINT (TWO REVERSALS) IF ( J > 1 ) EXIT ! THE LOOP ELSE ! LINE COMPLETED IN LOCAL COORDINATES EXIT ! THE LOOP END IF END IF END DO ! OVER J TRIES ! RETURN LOCAL COORD. OF LAST POINT PT_INOUT = PT_L (1:2, KOUNT) !--> CONVERT TO GLOBAL COORDINATES DO I = 1, KOUNT ! EVALUATE SHAPE FUNCTIONS CALL SCALAR_SHAPES (PT_L (:, I), H_L) ! EVALUATE GLOBAL COORDINATES, XYZ_PTS = H*COORD XYZ_PTS (1:N_SPACE, I) = MATMUL ( H_L (1:NOD_PER_EL), COORD ) END DO !--> CONTOUR FINISHED, RETURN FOR PLOTTING END SUBROUTINE CONTOUR_ISOPAR_SURF SUBROUTINE GET_VALUE_FROM_TABLE (T, VALUE, T_TABLE, V_TABLE, NT, NV) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! INTERPOLATION FOR VALUES AT T FROM TABULATED DATA ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP IMPLICIT NONE INTEGER, INTENT(IN) :: NT, NV REAL(DP), INTENT(IN) :: T, T_TABLE (NT), V_TABLE (NT, NV) REAL(DP), INTENT(OUT) :: VALUE (NV) INTEGER :: FIRST, I, LAST REAL(DP) :: DIFF, RATIO ! NT = NO INDEPENDENT VARIABLES IN T_TABLE ! NV = NO DEPENDENT VARS FOR EACH INDEPENDENT VARIABLE ! T = GIVEN INDEPENDENT VARIABLE ! T_TABLE = INDEPENDENT VARIABLE TABLE ! V_TABLE = DEPENDENT VARIABLE TABLE ! VALUE = INTERPOLATED DEPENDENT VAR. AT T IF ( T <= T_TABLE (1) ) THEN VALUE (:) = V_TABLE (1, :) ; RETURN END IF IF ( T >= T_TABLE (NT) ) THEN VALUE (:) = V_TABLE (NT, :) ; RETURN END IF DO I = 2, NT LAST = I IF ( T <= T_TABLE (I) ) EXIT ! THIS LOOP END DO FIRST = LAST - 1 DIFF = T_TABLE (LAST) - T_TABLE (FIRST) IF ( DIFF == 0.d0 ) THEN RATIO = 0.d0 ELSE RATIO = (T - T_TABLE (FIRST) ) / DIFF END IF VALUE (:) = (1.d0 - RATIO) * V_TABLE (FIRST, :) & + RATIO * V_TABLE (LAST, :) END SUBROUTINE GET_VALUE_FROM_TABLE SUBROUTINE G_DERIV (AJ_INV, DELTA, GLOBAL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! N_SPACE GLOBAL DERIVATIVES OF LT_FREE INTERPOLATION ! FUNCTIONS AT A LOCAL POINT. ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data REAL(DP), INTENT(IN) :: AJ_INV (N_SPACE, N_SPACE) REAL(DP), INTENT(IN) :: DELTA (N_SPACE, LT_FREE) REAL(DP), INTENT(OUT) :: GLOBAL (N_SPACE, LT_FREE) ! N_SPACE = DIMENSION OF SPACE ! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! AJ_INV = INVERSE JACOBIAN MATRIX AT LOCAL POINT ! DELTA = LOCAL COORD DERIV AT POINT OF INTEREST ! GLOBAL = GLOBAL DERIVATIVES MATRIX AT LOCAL POINT ! GLOBAL = AJ_INV*DELTA GLOBAL = MATMUL (AJ_INV, DELTA) END SUBROUTINE G_DERIV SUBROUTINE GLOBAL_DERIV_OF_H (GLOBAL) !b SUBROUTINE GLOBAL_DERIV_OF_H (COORD, DLH, GLOBAL) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! N_SPACE GLOBAL DERIVATIVES OF LT_FREE INTERPOLATION ! FUNCTIONS AT A LOCAL POINT. ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for COORD (LT_N, N_SPACE), DLH (N_SPACE, LT_N) !b REAL(DP), INTENT(IN) :: COORD (LT_N, N_SPACE) !b REAL(DP), INTENT(IN) :: DLH (N_SPACE, LT_N) REAL(DP), INTENT(OUT) :: GLOBAL (N_SPACE, LT_N) REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), DET ! AJ = JACOBIAN MATRIX AT LOCAL POINT ! AJ_INV = INVERSE JACOBIAN MATRIX AT LOCAL POINT ! COORD = COORDINATES OF ALL NODES ON ELEMENT ! DLH = LOCAL COORD DERIV AT POINT OF INTEREST ! GLOBAL = GLOBAL DERIVATIVES MATRIX AT LOCAL POINT ! LT_N = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE ! N_SPACE = DIMENSION OF SPACE ! AJ = DLH*COORD, GLOBAL = AJ_INV*DLH AJ = MATMUL (DLH, COORD) CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) GLOBAL = MATMUL (AJ_INV, DLH) END SUBROUTINE GLOBAL_DERIV_OF_H !wSUBROUTINE GLOBAL_2ND_DERIV_OF_H (GLOBAL, GLOBAL_2) !w!b SUBROUTINE GLOBAL_2ND_DERIV_OF_H (COORD, DLH, D2LH, & !w!b GLOBAL, GLOBAL_2) !w! * * * * * * * * * * * * * * * * * * * * * * * * * * ! N_SPACE GLOBAL FIRST AND SECOND DERIVATIVES OF LT_N !w! INTERPOLATION FUNCTIONS AT A LOCAL POINT. !w! * * * * * * * * * * * * * * * * * * * * * * * * * * !wUse System_Constants ! for N_SPACE, N_2_DER !wUse Elem_Type_Data ! for LT_N,COORD (LT_N, N_SPACE), !w! DLH (N_SPACE, LT_N), D2LH (N_2_DER, LT_N) !w!b REAL(DP), INTENT(IN) :: COORD (LT_N, N_SPACE) !w!b REAL(DP), INTENT(IN) :: DLH (N_SPACE, LT_N) !w!b REAL(DP), INTENT(IN) :: D2LH (N_2_DER, LT_N) !w !wREAL(DP), INTENT(OUT) :: GLOBAL (N_SPACE, LT_N) !wREAL(DP), INTENT(OUT) :: GLOBAL_2 (N_2_DER, LT_N) !w !wREAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), DET !wREAL(DP) :: P_AJ (N_2_DER, N_2_DER), P_AJ_INV (N_2_DER, N_2_DER) !wINTEGER :: I_ERROR = 0 !w !w! AJ = JACOBIAN MATRIX AT LOCAL POINT !w! AJ_INV = INVERSE JACOBIAN MATRIX AT LOCAL POINT !w! DLH = LOCAL COORD DERIV AT POINT OF INTEREST !w! D1LH = LOCAL COORD SECOND DERIV AT POINT OF INTEREST !w! GLOBAL = GLOBAL DERIVATIVES MATRIX AT LOCAL POINT !w! GLOBAL_2 = GLOBAL SECOND DERIVATIVES MATRIX AT LOCAL POINT !w! LT_N = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE !w! N_SPACE = DIMENSION OF SPACE !w! N_2_DER = NUMBER OF SECOND DERIVATIVES, 1, 3, OR 6 !w! P_AJ = PRODUCT OF JACOBIAN TERMS !w! P_AJ_INV = INVERSE OF PRODUCT OF JACOBIAN TERMS !w !w! AJ = DLH*COORD, GLOBAL = AJ_INV*DLH !wAJ = MATMUL (DLH, COORD) !wCALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) !wGLOBAL = MATMUL (AJ_INV, DLH) !w !w! NOW GET SECOND DERIVATIVE: !w! GLOBAL_2 = P_AJ_INV (D2HL - HESSIAN * GLOBAL) !w! HESSIAN = D2HL * COORD !wCALL JACOBIAN_PRODUCTS (AJ, P_AJ) !wCALL INV_SMALL_MAT (N_2_DER, P_AJ, P_AJ_INV, I_ERROR) !wIF ( I_ERROR > 0 ) THEN !wPRINT *, 'WARNING: GLOBAL_2ND_DERIV_OF_H, BAD INVERSE' !wN_WARN = N_WARN + 1 !wEND IF !wGLOBAL_2 = MATMUL (P_AJ_INV, & !w(D2LH - MATMUL ( MATMUL (D2LH, COORD), GLOBAL))) !wEND SUBROUTINE GLOBAL_2ND_DERIV_OF_H FUNCTION GET_REAL_IDENTITY (N) RESULT (EYE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! COMPUTE A REAL N BY N IDENTITY MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: N ! SIZE OF MATRIX REAL(DP) :: EYE (N, N) ! IDENTITY MATRIX INTEGER :: J REAL(DP), PARAMETER :: I_2 (2, 2) & = RESHAPE ((/ 1.d0, 0.d0, & 0.d0, 1.d0 /), & (/ 2, 2/) ) REAL(DP), PARAMETER :: I_3 (3, 3) & = RESHAPE ((/ 1.d0, 0.d0, 0.d0, & 0.d0, 1.d0, 0.d0, & 0.d0, 0.d0, 1.d0 /), & (/ 3, 3/) ) SELECT CASE (N) CASE (:0) ; STOP 'INVALID SIZE, GET_REAL_IDENTITY' CASE (1) ; EYE = 1.D0 CASE (2) ; EYE = I_2 CASE (3) ; EYE = I_3 CASE DEFAULT EYE = 0.D0 DO J = 1, N ; EYE (J, J) = 1.D0 ; END DO END SELECT END FUNCTION GET_REAL_IDENTITY SUBROUTINE REAL_IDENTITY (N, EYE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! COMPUTE A REAL N BY N IDENTITY MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: N ! SIZE OF MATRIX REAL(DP), INTENT(OUT) :: EYE (N, N) ! IDENTITY MATRIX INTEGER :: J REAL(DP), PARAMETER :: I_2 (2, 2) & = RESHAPE ((/ 1.d0, 0.d0, & 0.d0, 1.d0 /), & (/ 2, 2/) ) REAL(DP), PARAMETER :: I_3 (3, 3) & = RESHAPE ((/ 1.d0, 0.d0, 0.d0, & 0.d0, 1.d0, 0.d0, & 0.d0, 0.d0, 1.d0 /), & (/ 3, 3/) ) SELECT CASE (N) CASE (:0) ; STOP 'INVALID SIZE, GET_REAL_IDENTITY' CASE (1) ; EYE = 1.D0 CASE (2) ; EYE = I_2 CASE (3) ; EYE = I_3 CASE DEFAULT EYE = 0.D0 DO J = 1, N ; EYE (J, J) = 1.D0 ; END DO END SELECT END SUBROUTINE REAL_IDENTITY SUBROUTINE INVERT_JACOBIAN (AJ, AJ_INV, DET, P_SPACE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND INVERSE AND DETERMINANT OF GEOMETRY JACOBIAN ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN ) :: P_SPACE REAL(DP), INTENT(IN ) :: AJ (P_SPACE, P_SPACE) REAL(DP), INTENT(OUT) :: AJ_INV (P_SPACE, P_SPACE), DET INTEGER :: I_ERROR = 0 ! P_SPACE = NUMBER OF SPATIAL DIMENSIONS ! AJ = JACOBIAN MATRIX AT A POINT ! AJ_INV = INVERSE OF AJ ! DET = DETERMINANT OF AJ SELECT CASE (P_SPACE) !--> 1-D CASE (1) DET = AJ(1,1) IF ( ABS(DET) <= SPACING(DET) ) THEN WRITE (N_PRT, *) 'WARNING: INVERT 1x1, NEARLY SINGULAR' WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET N_WARN = N_WARN + 1 I_ERROR = 1 !b IF ( DET == 0.d0 ) STOP 'BAD DET, INVERT_JACOBIAN' ELSE AJ_INV(1,1) = 1.d0/DET END IF !--> 2-D CASE (2) CALL INVERT_2BY2 (AJ, AJ_INV, DET, I_ERROR) !--> 3-D CASE (3) CALL INVERT_3BY3 (AJ, AJ_INV, DET, I_ERROR) !--> Error CASE DEFAULT PRINT *,'FATAL ERROR IN ELEMENT ', THIS_EL STOP 'BAD PARAMETRIC SPACE, INVERT_JACOBIAN' END SELECT IF ( DET <= 0.d0 ) THEN WRITE (N_PRT, *) 'WARNING: CORRECTED NEGATIVE JACOBIAN' WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET N_WARN = N_WARN + 1 END IF END SUBROUTINE INVERT_JACOBIAN SUBROUTINE INVERT_3BY3 (A, AINV, DET, I_ERROR) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND INVERSE AND DETERMINANT OF MATRIX A(3,3) ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: A(3,3) REAL(DP), INTENT(OUT) :: AINV(3,3), DET INTEGER, INTENT(INOUT) :: I_ERROR ! A = ORIGINAL MATRIX ! AINV = INVERSE OF MATRIX A ! DET = DETERMINANT OF A AINV(1,1) = A(2,2)*A(3,3) - A(3,2)*A(2,3) AINV(2,1) = -A(2,1)*A(3,3) + A(3,1)*A(2,3) AINV(3,1) = A(2,1)*A(3,2) - A(3,1)*A(2,2) AINV(1,2) = -A(1,2)*A(3,3) + A(3,2)*A(1,3) AINV(2,2) = A(1,1)*A(3,3) - A(3,1)*A(1,3) AINV(3,2) = -A(1,1)*A(3,2) + A(3,1)*A(1,2) AINV(1,3) = A(1,2)*A(2,3) - A(2,2)*A(1,3) AINV(2,3) = -A(1,1)*A(2,3) + A(2,1)*A(1,3) AINV(3,3) = A(1,1)*A(2,2) - A(2,1)*A(1,2) DET = A(1,1)*AINV(1,1) + A(1,2)*AINV(2,1) & + A(1,3)*AINV(3,1) IF ( ABS(DET) <= SPACING(DET) ) THEN WRITE (N_PRT, *) 'WARNING: INVERT_3BY3, NEARLY SINGULAR MATRIX' WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET N_WARN = N_WARN + 1 I_ERROR = 1 !b IF ( DET == 0.d0 ) STOP 'SINGULAR 3X3 MATRIX IN INVERT_3BY3' ELSE AINV = AINV/DET END IF END SUBROUTINE INVERT_3BY3 SUBROUTINE INVERT_2BY2 (A, AINV, DET, I_ERROR) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE THE DETERMINANT AND INVERSE OF A(2,2) ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants implicit none REAL(DP), INTENT(IN) :: A (2, 2) REAL(DP), INTENT(OUT) :: AINV (2, 2), DET INTEGER, INTENT(INOUT) :: I_ERROR ! A = ORIGINAL MATRIX ! AINV = INVERSE OF MATRIX A ! DET = DETERMINANT OF A DET = A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1) IF ( ABS(DET) <= SPACING(DET) ) THEN WRITE (N_PRT, *) 'WARNING: INVERT_2BY2, NEARLY SINGULAR MATRIX' WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET N_WARN = N_WARN + 1 I_ERROR = 1 END IF IF ( DET == 0.d0 ) THEN PRINT *, "SINGULAR MATRIX:" ; CALL RPRINT (A, 2, 2, 1) STOP 'SINGULAR 2X2 MATRIX, INVERT_2BY2' ELSE AINV (1, 1) = A (2, 2) / DET AINV (1, 2) = - A (1, 2) / DET AINV (2, 1) = - A (2, 1) / DET AINV (2, 2) = A (1, 1) / DET END IF ! singular END SUBROUTINE INVERT_2BY2 SUBROUTINE JACOBIAN_AT_PT (DELTA, COORD, AJ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE THE JACOBIAN MATRIX AT A LOCAL POINT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_SPACE, NOD_PER_EL IMPLICIT NONE REAL(DP), INTENT(IN) :: DELTA (N_SPACE, NOD_PER_EL) REAL(DP), INTENT(IN) :: COORD (NOD_PER_EL, N_SPACE) REAL(DP), INTENT(OUT) :: AJ (N_SPACE, N_SPACE) ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_SPACE = DIMENSION OF SPACE ! DELTA = LOCAL DERIVATIVES OF NOD_PER_EL INTERPOLATION ! FUNCTIONS AT POINT OF INTEREST. ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! AJ = JACOBIAN MATRIX = DELTA*COORD AJ = MATMUL (DELTA, COORD) END SUBROUTINE JACOBIAN_AT_PT FUNCTION GEOMETRIC_JACOBIAN () RESULT (AJ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET JACOBIAN MATRIX AT A POINT, LT_GEOM <= LT_N ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For N_SPACE Use Elem_Type_Data ! for: LT_GEOM, LT_PARM, ! GEOMETRY (LT_GEOM, N_SPACE), ! DLG (LT_PARM, LT_GEOM) IMPLICIT NONE REAL(DP) :: AJ (LT_PARM, LT_PARM) ! Out ! LT_GEOM = NUMBER OF NODES PER ELEMENT ! LT_PARM = DIMENSION OF PARAMETRIC SPACE ! N_SPACE = DIMENSION OF SPACE ! DLG = LOCAL DERIVATIVES OF LT_GEOM INTERPOLATION ! FUNCTIONS AT POINT OF INTEREST. ! GEOMETRY = SPATIAL COORDINATES OF GEOMETRIC NODES ! AJ = JACOBIAN MATRIX = DLG*GEOMETRY AJ = MATMUL (DLG, GEOMETRY (:, 1:LT_PARM)) END FUNCTION GEOMETRIC_JACOBIAN SUBROUTINE DIAGONALIZE_SQ_MATRIX (NL, SQ, DSQ) ! * * * * * * * * * * * * * * * * * * * * * ! Diagonialize a consistent square matrix ! * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NL REAL(DP), INTENT(IN) :: SQ (NL, NL) REAL(DP), INTENT(OUT) :: DSQ (NL) INTEGER :: I REAL(DP) :: DIAG, TOTAL, RATIO ! NL = size of array ! SQ = Original square matrix ! DSQ = Diagonal form DIAG = 0.0_DP TOTAL = SUM ( SQ ) ! sum of all terms DO I = 1, NL ! Count diagonal sum DIAG = DIAG + SQ (I, I) END DO ! Scale diagonal IF (DIAG <= 0.0_DP) STOP 'ERROR IN DIAGONALIZE_SQ_MATRIX' RATIO = TOTAL / DIAG DO I = 1, NL DSQ (I) = RATIO * SQ (I, I) END DO END SUBROUTINE DIAGONALIZE_SQ_MATRIX SUBROUTINE LINE_2D_GEOM (I_OPT, COORD, DX, DY, DS, & X_NORM, Y_NORM, X_TAN, Y_TAN) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! GEOMETRIC PROPERTIES OF A STRAIGHT LINE IN X-Y SPACE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for DP IMPLICIT NONE INTEGER, INTENT(IN) :: I_OPT REAL(DP), INTENT(IN) :: COORD (2, 2) REAL(DP), INTENT(OUT) :: DX, DY, DS, X_NORM, Y_NORM, X_TAN, Y_TAN ! COORD = GIVEN COORDINATES OF TWO DEFINING POINTS ! DX,DY,DS = X,Y COMPONENTS AND TOTAL LENGTH N ! X_NORM,Y_NORM = COMPONENTS OF UNIT NORMAL 1>>>2 V ! X_TAN,Y_TAN = COMPONENTS OF UNIT TANGENT FROM 1 TO 2 ! I_OPT = 0, CALCULATE DISTANCES ONLY ! I_OPT = 1, CALCULATE DISTANCES AND NORMAL VECTOR ! I_OPT = 2, CALCULATE DIST, NORMAL, TANGENT DX = COORD (2, 1) - COORD (1, 1) DY = COORD (2, 2) - COORD (1, 2) DS = SQRT (DX * DX + DY * DY) IF (I_OPT <= 0) RETURN IF (DS == 0.d0) THEN WRITE (6, *) 'ZERO LENGTH IN LINE_2D_GEOM' STOP 'ZERO LENGTH IN LINE_2D_GEOM' ELSE X_NORM = DY / DS Y_NORM = - DX / DS IF (I_OPT == 1) RETURN X_TAN = DX / DS Y_TAN = DY / DS END IF END SUBROUTINE LINE_2D_GEOM SUBROUTINE SQ_TO_DIAGONAL_SQ_MATRIX (NL, SQ, DSQ) ! * * * * * * * * * * * * * * * * * * * * * ! Diagonialize a consistent square matrix ! * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NL REAL(DP), INTENT(IN) :: SQ (NL, NL) REAL(DP), INTENT(OUT) :: DSQ (NL, NL) INTEGER :: I REAL(DP) :: DIAG, TOTAL, RATIO ! NL = size of array ! SQ = Original square matrix ! DSQ = Diagonal form DIAG = 0.0_DP TOTAL = SUM ( SQ ) ! sum of all terms DO I = 1, NL ! Count diagonal sum DIAG = DIAG + SQ (I, I) END DO ! Scale diagonal IF (DIAG <= 0.0_DP) STOP 'ERROR IN SQ_TO_DIAGONAL_SQ_MATRIX' RATIO = TOTAL / DIAG DO I = 1, NL DSQ (I, I) = RATIO * SQ (I, I) END DO END SUBROUTINE SQ_TO_DIAGONAL_SQ_MATRIX FUNCTION OUTER_PRODUCT (A, B) RESULT (C) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OUTER PRODUCT OF TWO VECTORS, C_IJ = A_I * B_J ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module REAL(DP), INTENT(IN) :: A (:), B (:) REAL(DP) :: C (SIZE(A), SIZE(B)) C = SPREAD (A, DIM=2, NCOPIES=SIZE(B)) & * SPREAD (B, DIM=1, NCOPIES=SIZE(A)) END FUNCTION OUTER_PRODUCT FUNCTION OUTER_DIVIDE (A, B) RESULT (C) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OUTER QUOTIENT OF TWO VECTORS, C_IJ = A_I / B_J ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module REAL(DP), INTENT(IN) :: A (:), B (:) REAL(DP) :: C (SIZE(A), SIZE(B)) C = SPREAD (A, DIM=2, NCOPIES=SIZE(B)) & / SPREAD (B, DIM=1, NCOPIES=SIZE(A)) END FUNCTION OUTER_DIVIDE FUNCTION OUTER_AND (A, B) RESULT (C) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! OUTER LOGICAL 'AND' OF TWO VECTORS, C_IJ = A_I & B_J ! * * * * * * * * * * * * * * * * * * * * * * * * * * LOGICAL, INTENT(IN) :: A (:), B (:) LOGICAL :: C (SIZE(A), SIZE(B)) C = SPREAD (A, DIM=2, NCOPIES=SIZE(B)) & .AND. SPREAD (B, DIM=1, NCOPIES=SIZE(A)) END FUNCTION OUTER_AND !SUBROUTINE NGRAND (WT, DET, V, DGV, XPT, NOD_PER_EL, N_SPACE, & ! N_EL_FRE, COL, SQ, N_FILE ) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! PROBLEM DEPENDENT INTEGRAND EVALUATION IN !! AN ISOPARAMETRIC OR SUBPARAMETRIC ELEMENT !! * * * * * * * * * * * * * * * * * * * * * * * * * * !Use Precision_Module !IMPLICIT NONE !INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_EL_FRE, N_FILE !REAL(DP), INTENT(IN) :: V (N_EL_FRE), DGV (N_SPACE, N_EL_FRE) !REAL(DP), INTENT(IN) :: XPT (N_SPACE), WT, DET !REAL(DP), INTENT(INOUT) :: COL (N_EL_FRE), SQ (N_EL_FRE, N_EL_FRE) ! !! NOD_PER_EL = NUMBER OF NODES PER ELEMENT !! N_SPACE = NUMBER OF SPATIAL DIMENSIONS !! N_EL_FRE = NUMBER OF ELEMENT DEGREES OF FREEDOM !! V = ELEMENT INTERPOLATION FUNCTIONS !! DGV = GLOBAL DERIVATIVES OF V !! XPT = GLOBAL COORDS OF THE POINT !! WT = QUADRATURE WEIGHT AT POINT !! DET = JACOBIAN DETERMINANT AT POINT !! COL = PROB DEP COLUMN MATRIX INTEGRAND !! SQ = PROB DEP SQUARE MATRIX INTEGRAND !! N_FILE = STORAGE UNIT FOR POST SOLUTION DATA !! .................................................... !! *** NGRAND PROBLEM DEPENDENT STATEMENTS FOLLOW *** !! .................................................... ! !inc include 'my_ngrand_inc' !END SUBROUTINE NGRAND !FUNCTION COSH (X) RESULT (VALUE) ! Hyperbolic cosine !Use Precision_module !IMPLICIT NONE ! REAL(DP), INTENT(IN) :: X ! REAL(DP) :: VALUE ! VALUE = (EXP (X) + EXP(-X)) * 0.5d0 !END FUNCTION COSH ! !FUNCTION SINH (X) RESULT (VALUE) ! Hyperbolic sine !Use Precision_module !IMPLICIT NONE ! REAL(DP), INTENT(IN) :: X ! REAL(DP) :: VALUE ! VALUE = (EXP (X) - EXP(-X)) * 0.5d0 !END FUNCTION SINH SUBROUTINE SYMMETRIC_INVERSE (A, A_INV, I_ERROR) ! ------------------------------------------------------ ! INVERT SYMMETRIC MATRIX A(N,N) ! ------------------------------------------------------ Use Precision_Module IMPLICIT NONE REAL (DP), INTENT (IN) :: A (:, :) REAL (DP), INTENT (OUT) :: A_INV (SIZE(A, 1), SIZE(A, 1)) INTEGER, INTENT (OUT) :: I_ERROR ! FATAL RESULT IF /= 0 INTEGER :: I, J, K, N REAL (DP) :: D ! A = ORIGINAL MATRIX ! A_INV = INVERSE OF MATRIX A ! D = A DIAGONAL ELEMENT OF A ! N = SIZE OF A I_ERROR = 0 ; N = SIZE(A, 1) ; A_INV = A ! INITIALIZE DO K = 1,N ! LOOP OVER ROWS D = A_INV (K,K) ! DIAGONAL IS PIVOT TERM !b IF ( D == 0.d0 ) THEN ! FATAL CONDITION IF ( ABS(D) <= SPACING(D) ) THEN ! FATAL CONDITION PRINT *, 'ERROR: ZERO PIVOT IN SYMMETRIC_INVERSE' I_ERROR = 1 ; RETURN ! ABORT END IF ! FATAL ERROR A_INV (K, :) = -A_INV (K, :)/D ! SCALE ROW DO I = 1, N ! LOOP OVER ROW AGAIN IF ( I /= K ) THEN ! OFF-DIAGONAL TERMS DO J = 1, N ! LOOP OVER COLUMNS IF ( J /= K ) A_INV (I,J) = & A_INV (I,J) + A_INV (I,K)*A_INV (K,J) END DO ! OVER J END IF ! OFF DIAGONAL I A_INV (I,K) = A_INV (I,K)/D END DO ! OVER ROW I A_INV (K,K) = 1.d0/D END DO ! OVER K-TH COLUMN END SUBROUTINE SYMMETRIC_INVERSE SUBROUTINE INVERT_SMALL_MATRIX (A, A_INV, I_ERROR) ! ------------------------------------------------------ ! INVERT A SMALL MATRIX A(N,N) ! ------------------------------------------------------ Use Precision_Module IMPLICIT NONE REAL (DP), INTENT (IN) :: A (:, :) REAL (DP), INTENT (OUT) :: A_INV (SIZE(A, 1), SIZE(A, 1)) INTEGER, INTENT (OUT) :: I_ERROR ! FATAL RESULT IF /= 0 INTEGER :: I, K, N REAL (DP) :: D ! A = ORIGINAL MATRIX ! A_INV = INVERSE OF MATRIX A ! D = A DIAGONAL ELEMENT OF A ! N = SIZE OF A I_ERROR = 0 ;N = SIZE(A, 1) ; A_INV = A ! INITIALIZE DO K = 1,N ! LOOP OVER DIAGONALS D = A_INV (K, K) ! PIVOT TERM !b IF ( D == 0.d0 ) THEN ! FATAL CONDITION IF ( ABS(D) <= SPACING(D) ) THEN ! FATAL CONDITION PRINT *, 'ERROR: ZERO PIVOT IN INVERT_SMALL_MATRIX' I_ERROR = 1 ; RETURN ! ABORT END IF ! FATAL ERROR A_INV (K, K) = 1.d0 A_INV (K, :) = A_INV (K, :)/D ! SCALE ROW DO I = 1, N ! LOOP OVER ROW AGAIN IF ( I == K ) CYCLE ! TO NEXT I D = A_INV (I, K) A_INV (I, K) = 0.d0 A_INV (I, :) = A_INV (I, :) - D*A_INV (K, :) END DO ! OVER ROW I END DO ! OVER K-TH COLUMN END SUBROUTINE INVERT_SMALL_MATRIX SUBROUTINE INV_SMALL_MAT (N, A, A_INV, I_ERROR) ! ------------------------------------------------------ ! INVERT A SMALL MATRIX A(N,N) ! ------------------------------------------------------ Use System_Constants IMPLICIT NONE INTEGER, INTENT (IN) :: N REAL (DP), INTENT (IN) :: A (N, N) REAL (DP), INTENT (OUT) :: A_INV (N, N) INTEGER, INTENT (OUT) :: I_ERROR ! FATAL RESULT IF /= 0 INTEGER :: I, K REAL (DP) :: D, DET ! A = ORIGINAL MATRIX ! A_INV = INVERSE OF MATRIX A ! D = A DIAGONAL ELEMENT OF A ! N = SIZE OF A I_ERROR = 0 ; A_INV = A ! INITIALIZE SELECT CASE (N) !--> 1-D CASE (1) DET = A (1,1) IF ( ABS(DET) <= SPACING(DET) ) THEN WRITE (N_PRT, *) 'WARNING: INVERT 1x1, NEARLY SINGULAR' WRITE (N_PRT, *) 'DET = ', DET N_WARN = N_WARN + 1 I_ERROR = 1 ELSE A_INV(1,1) = 1.d0/DET END IF CASE (2) !--> 2-D CALL INVERT_2BY2 (A, A_INV, DET, I_ERROR) CASE (3) !--> 3-D CALL INVERT_3BY3 (A, A_INV, DET, I_ERROR) CASE DEFAULT DO K = 1,N ! LOOP OVER DIAGONALS D = A_INV (K, K) ! PIVOT TERM IF ( ABS(D) <= SPACING(D) ) THEN ! FATAL CONDITION PRINT *, 'WARNING: ZERO PIVOT IN INV_SMALL_MAT' N_WARN = N_WARN + 1 I_ERROR = 1 EXIT ! THIS DO ON K !b END IF ! FATAL ERROR A_INV (K, K) = 1.d0 A_INV (K, :) = A_INV (K, :)/D ! SCALE ROW DO I = 1, N ! LOOP OVER ROW AGAIN IF ( I == K ) CYCLE ! TO NEXT I D = A_INV (I, K) A_INV (I, K) = 0.d0 A_INV (I, :) = A_INV (I, :) - D*A_INV (K, :) END DO ! OVER ROW I END DO ! OVER K-TH COLUMN END SELECT IF ( I_ERROR > 0 ) THEN PRINT *, 'ERROR: INV_SMALL_MAT FAILED, IDENTITY MATRIX RETURNED' CALL REAL_IDENTITY (N, A_INV) END IF ! INVERSION FAILED END SUBROUTINE INV_SMALL_MAT SUBROUTINE JACOBIAN_PRODUCTS (AJ, P_AJ) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET JACOBIAN PRODUCTS FOR SECOND DERIVATIVE CALCULATION !* * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL (DP), INTENT (IN) :: AJ (N_SPACE, N_SPACE) REAL (DP), INTENT (OUT) :: P_AJ (N_2_DER, N_2_DER) ! AJ = JACOBIAN MATRIX AT A POINT ! P_AJ = JACOBIAN PRODUCTS AT POINT IN TERMS OF THE 2ND GLOBAL ! DERIVATIVES (GENERALLY WANT P_AJ INVERSE) ! V_,AB = P_AJ_IAJB * V_,IJ (INDEX NOTATION: _ PRECEEDS SUBSCRIPTS) ! WHERE P_AJ_IAJB = X_I,A * X_J,B SELECT CASE ( N_2_DER ) ! NUMBER OF SECOND DERIVATIVES CASE ( 1 ) ! One-dimensional space ! | V,AA | = (P_AJ) | V,II | A,B IMPLY LOCAL COORD ! I,J IMPLY GLOBAL COORD P_AJ (1, 1) = AJ (1, 1) **2 CASE ( 3 ) ! Two-dimensional space ! | V,AA | | V,II | A,B IMPLY LOCAL COORD ! | V,AB | = (P_AJ) | V,IJ | I,J IMPLY GLOBAL COORD ! | V,BB | | V,JJ | P_AJ (1, 1) = AJ (1, 1) **2 P_AJ (2, 1) = AJ (1, 1) * AJ (2, 1) P_AJ (3, 1) = AJ (2, 1) **2 P_AJ (1, 2) = 2.d0 * AJ (1, 1) * AJ (1, 2) P_AJ (2, 2) = AJ (1, 1) * AJ (2, 2) + AJ (2, 1) * AJ (1, 2) P_AJ (3, 2) = 2.d0 * AJ (2, 1) * AJ (2, 2) P_AJ (1, 3) = AJ (1, 2) **2 P_AJ (2, 3) = AJ (1, 2) * AJ (2, 2) P_AJ (3, 3) = AJ (2, 2) **2 CASE ( 6 ) ! Three-dimensional space ! | V,AA | | V,II | A,B,C IMPLY LOCAL ! | V,AB | | V,IJ | I,J,K IMPLY GLOBAL ! | V,AC | = (P_AJ) | V,IK | ! | V,BB | | V,JJ | ! | V,BC | | V,JK | ! | V,CC | | V,KK | P_AJ (1, 1) = AJ (1, 1) **2 P_AJ (2, 1) = AJ (1, 2) * AJ (2, 2) P_AJ (3, 1) = AJ (1, 1) * AJ (3, 1) P_AJ (4, 1) = AJ (2, 1) **2 P_AJ (5, 1) = AJ (2, 1) * AJ (3, 1) P_AJ (6, 1) = AJ (3, 1) **2 P_AJ (1, 2) = AJ (1, 1) * AJ (1, 2) * 2.d0 P_AJ (2, 2) = AJ (1, 1) * AJ (2, 2) + AJ (1, 2) * AJ (2, 1) P_AJ (3, 2) = AJ (1, 1) * AJ (3, 2) + AJ (1, 2) * AJ (3, 1) P_AJ (4, 2) = AJ (2, 1) * AJ (2, 2) * 2.d0 P_AJ (5, 2) = AJ (2, 1) * AJ (3, 2) + AJ (2, 2) * AJ (3, 1) P_AJ (6, 2) = AJ (3, 1) * AJ (3, 2) * 2.d0 P_AJ (1, 3) = AJ (1, 1) * AJ (1, 3) * 2.d0 P_AJ (2, 3) = AJ (1, 1) * AJ (2, 3) + AJ (1, 3) * AJ (2, 1) P_AJ (3, 3) = AJ (1, 1) * AJ (3, 3) + AJ (1, 3) * AJ (3, 1) P_AJ (4, 3) = AJ (2, 1) * AJ (2, 3) * 2.d0 P_AJ (5, 3) = AJ (2, 1) * AJ (3, 3) + AJ (2, 3) * AJ (3, 1) P_AJ (6, 3) = AJ (3, 1) * AJ (3, 3) * 2.d0 P_AJ (1, 4) = AJ (1, 2) **2 P_AJ (2, 4) = AJ (1, 2) * AJ (2, 2) P_AJ (3, 4) = AJ (1, 2) * AJ (3, 2) P_AJ (4, 4) = AJ (2, 2) **2 P_AJ (5, 4) = AJ (2, 2) * AJ (3, 2) P_AJ (6, 4) = AJ (3, 2) **2 P_AJ (1, 5) = AJ (1, 2) * AJ (1, 3) * 2.d0 P_AJ (2, 5) = AJ (1, 2) * AJ (2, 3) + AJ (1, 3) * AJ (2, 2) P_AJ (3, 5) = AJ (1, 2) * AJ (3, 3) + AJ (1, 3) * AJ (3, 2) P_AJ (4, 5) = AJ (2, 2) * AJ (2, 3) * 2.d0 P_AJ (5, 5) = AJ (2, 2) * AJ (3, 3) + AJ (2, 3) * AJ (3, 2) P_AJ (6, 5) = AJ (3, 2) * AJ (3, 3) * 2.d0 P_AJ (1, 6) = AJ (1, 3) **2 P_AJ (2, 6) = AJ (1, 3) * AJ (2, 3) P_AJ (3, 6) = AJ (1, 3) * AJ (3, 3) P_AJ (4, 6) = AJ (2, 3) **2 P_AJ (5, 6) = AJ (2, 3) * AJ (3, 3) P_AJ (6, 6) = AJ (3, 3) **2 CASE DEFAULT PRINT *, 'JACOBIAN_PRODUCTS: UNKNOWN N_2_DER = ', N_2_DER STOP 'JACOBIAN_PRODUCTS: UNKNOWN N_2_DER' END SELECT END SUBROUTINE JACOBIAN_PRODUCTS subroutine Sort_Reals (lines, database) !------------------------------------------------------------ ! Bubble Sort of (changed) Real Database !------------------------------------------------------------ use System_Constants ! for DP implicit none integer, intent(in) :: lines ! number of records !b real, intent(inout) :: database (lines) ! records in database real(DP), intent(inout) :: database (lines) ! records in database integer :: swaps_Made ! number of swaps made in one pass integer :: count ! loop variable !b real :: temp ! temporary holder for making swap real(DP):: temp ! temporary holder for making swap do ! Repeat this loop forever... (until we break out of it) swaps_Made = 0 ! Initially, we've made no swaps ! Make one pass of the bubble sort algorithm do count = 1, (lines - 1) ! If item is greater than the one after it, swap them if ( database (count) > database (count + 1) ) then temp = database (count) database (count) = database (count + 1) database (count + 1) = temp swaps_Made = swaps_Made + 1 end if end do ! If we made no swaps, break out of the loop. if ( swaps_Made == 0 ) exit ! do count swaps end do end subroutine Sort_Reals subroutine Real_Vector_Order (lines, database, order) !------------------------------------------------------------ ! Ordered Bubble Sort of (Unchanged) Real Vector !------------------------------------------------------------ use System_Constants ! for DP implicit none integer, intent(in) :: lines ! number of records real(dp), intent(in) :: database (lines) ! records in database integer, intent(out) :: order (lines) ! the order array integer :: swaps_Made ! number of swaps made in one pass integer :: count ! loop variable integer :: temp ! temporary holder for making swap order = (/ (count, count = 1, lines) /) ! default order do ! Repeat this loop forever... (until we break out of it) swaps_Made = 0 ! Initially, we've made no swaps ! Make one pass of the bubble sort algorithm do count = 1, (lines - 1) ! If item is greater than the one after it, swap them if ( database (order (count)) > & database (order (count + 1)) ) then temp = order (count) order (count) = order (count + 1) order (count + 1) = temp swaps_Made = swaps_Made + 1 end if end do ! If we made no swaps, break out of the loop. if ( swaps_Made == 0 ) exit ! do count swaps end do end subroutine Real_Vector_Order subroutine Sort_string (lines, database) !------------------------------------------------------------ ! Bubble Sort of (Changed) String Database !------------------------------------------------------------ implicit none integer, intent(in) :: lines ! input size character(len=*), intent(inout) :: database (lines) ! records character (len = len(database (1))) :: temp ! swap holder integer :: swaps_Made ! number of swaps in a pass integer :: count ! loop variable interface ! to_lower function to_lower (string) result (new_string) character (len = *), intent(in) :: string character (len = len(string)) :: new_string end function to_lower end interface ! to_lower do ! Repeat this loop forever... (until we break out of it) swaps_Made = 0 ! Initially, we've made no swaps ! Make one pass of the bubble sort algorithm do count = 1, (lines - 1) ! If the element is greater than the one after it, swap them if ( LGT (to_lower (database (count )), & to_lower (database (count + 1))) ) then temp = database (count ) database (count ) = database (count + 1) database (count + 1) = temp swaps_Made = swaps_Made + 1 end if end do ! If we made no swaps, berak out of the loop. if ( swaps_Made == 0) exit ! do count swaps end do end subroutine Sort_string function to_lower (string) result (new_string) !------------------------------------------------------------ ! Convert a string or character to lower case ! (valid for ASCII or EBCDIC processors) !------------------------------------------------------------ implicit none character (len = *), intent(in) :: string ! unknown length character (len = len(string)) :: new_string ! same length character (len=26), parameter :: & UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ', & lower = 'abcdefghijklmnopqrstuvwxyz' integer :: k ! loop counter integer :: loc ! position in alphabet new_string = string ! copy everything do k = 1, len(string) ! to change letters loc = index ( UPPER, string(k:k)) if (loc /= 0 ) new_string(k:k) = lower(loc:loc) end do ! over string characters end function to_lower subroutine Integer_Sort (lines, database, order) !------------------------------------------------------------ ! Ordered Bubble Sort of (Unchanged) Integer Database !------------------------------------------------------------ implicit none integer, intent(in) :: lines ! number of records integer, intent(in) :: database (lines) ! records in database integer, intent(out) :: order (lines) ! the order array integer :: swaps_Made ! number of swaps made in one pass integer :: count ! loop variable integer :: temp ! temporary holder for making swap order = (/ (count, count = 1, lines) /) ! default order do ! Repeat this loop forever... (until we break out of it) swaps_Made = 0 ! Initially, we've made no swaps ! Make one pass of the bubble sort algorithm do count = 1, (lines - 1) ! If item is greater than the one after it, swap them if ( database (order (count)) > & database (order (count + 1)) ) then temp = order (count) order (count) = order (count + 1) order (count + 1) = temp swaps_Made = swaps_Made + 1 end if end do ! If we made no swaps, break out of the loop. if ( swaps_Made == 0 ) exit ! do count swaps end do end subroutine Integer_Sort FUNCTION GET_ABS_NORM (A) RESULT (VALUE) Use System_Constants ! for DP IMPLICIT NONE REAL(DP), INTENT (IN) :: A (:, :) REAL(DP) :: VALUE VALUE = SUM ( ABS ( A ) ) END FUNCTION GET_ABS_NORM FUNCTION GET_INF_NORM (A) RESULT (VALUE) Use System_Constants ! for DP IMPLICIT NONE REAL(DP), INTENT (IN) :: A (:, :) REAL(DP) :: VALUE, TEST INTEGER :: I VALUE = 0.d0 DO I = 1, SIZE (A, 1) TEST = SUM ( ABS ( A (I, :) ) ) IF ( TEST > VALUE ) VALUE = TEST END DO ! over rows END FUNCTION GET_INF_NORM FUNCTION GET_MAX_NORM (A) RESULT (VALUE) Use System_Constants ! for DP IMPLICIT NONE REAL(DP), INTENT (IN) :: A (:, :) REAL(DP) :: VALUE VALUE = MAXVAL ( ABS (A) ) END FUNCTION GET_MAX_NORM FUNCTION GET_SQ_NORM (A) RESULT (VALUE) Use System_Constants ! for DP IMPLICIT NONE REAL(DP), INTENT (IN) :: A (:, :) REAL(DP) :: VALUE VALUE = SQRT ( SUM ( A **2) ) / SIZE (A, 1) END FUNCTION GET_SQ_NORM FUNCTION GET_1_NORM (A) RESULT (VALUE) Use System_Constants ! for DP IMPLICIT NONE REAL(DP), INTENT (IN) :: A (:, :) REAL(DP) :: VALUE, TEST INTEGER :: I VALUE = 0.d0 DO I = 1, SIZE (A, 2) TEST = SUM ( ABS ( A (:, I) ) ) IF ( TEST > VALUE ) VALUE = TEST END DO ! over rows END FUNCTION GET_1_NORM