! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE ASSEMBLE_SKYLINE_EQS (I_DIAG, NODES, S_SKY, CC, X, DD_OLD, & S_SKY_LOWER, ITER, M_SKY, M_SKY_LOWER) !b S_SKY_LOWER, ITER, M_SKY) !b X, DD_OLD, S_SKY_LOWER, ITER, I_BC) !b !b X, DD_OLD, S_SKY_LOWER, ITER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ASSEMBLE SYSTEM EQUATIONS AND STORE POST ! SOLUTION ELEMENT DATA ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_FREE), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_FREE, LT_QP), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! EL_M (LT_FREE, LT_FREE), DIAG_M (LT_FREE), ! D (LT_FREE), ELEM_NODES (LT_N) Use Sys_Properties_Data Use Interface_Header IMPLICIT NONE ! ALWAYS USED INTEGER, INTENT(IN) :: ITER ! Current iteration number INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) !b INTEGER, INTENT(IN) :: I_BC (MAX_NP) !b REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) !b REAL(DP), INTENT(IN) :: DD_OLD (N_D_FRE) !b REAL(DP), INTENT(INOUT) :: CC (N_D_FRE) REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF) REAL(DP), INTENT(INOUT) :: S_SKY_LOWER (NOT_SYM) REAL(DP), INTENT(INOUT) :: M_SKY (M_COEFF) REAL(DP), INTENT(INOUT) :: M_SKY_LOWER (NOT_SYM) ! effective unsym mass ! Locals INTEGER :: IE, J, LT REAL(DP) :: COL_SUMS ! Error checking ! Work items for time integration REAL(DP) :: K (LT_FREE, LT_FREE), M (LT_FREE, LT_FREE) REAL(DP) :: F (LT_FREE) ! VARIABLES: ! CC = SYSTEM EQUATIONS COLUMN MATRIX ! COL_SUMS + SUM OF ALL COLUMN MATRICES ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD_OLD = SYSTEM NODAL PARAMETERS FROM LAST ITERATION ! DELETED = TRUE IF THE ELEMENT HAS BEEN DELETED ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! DIAG_M = DIAGONALIZED ELEMENT MASS MATRIX (LT_FREE) ! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION ! DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = CONSTITUTIVE MATRIX ! EB = PRODUCT OF E*B ! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT ! EL_M = ELEMENT MASS MATRIX ! FLT_EL = REAL PROPERTIES OF ELEMENTS ! FLT_MISC = MISCELLANEOUS REAL PROPERTIES OF SYSTEM ! FLT_NP = REAL PROPERTIES OF SYSTEM NODES ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! ITER = CURRENT ITERATION NUMBER (not used) !b ! I_BC = NODAL BOUNDARY RESTRAINT INDICATOR ARRAY !b ! I_DIAG = DIAGONAL LOCATION IN SKYLINE VECTOR ! INDEX = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT ! L_HOMO = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS ! L_S_TOT = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS ! LP_FIX = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES ! L_PROP = ARRAY INTEGER ELEMENT PROPERTIES ! LP_TEST > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED ! L_SHAPE = SHAPE FLAG FOR QUADRATURE RULE SELECTION ! LT_DATA = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE ! LT = ELEMENT TYPE NUMBER (IF USED) ! LT_QP = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE ! LT_USER = OPTIONAL USER CODE ASSOCIATED WITH TYPE ! L_TYPE = ELEMENT TYPE NUMBER ! MAX_NP = NUMBER OF SYSTEM NODES ! MODE = STORAGE FLAG, 0=SYMMETRIC, OR 1=UNSYMMETRIC SKYLINE ! MISC_FIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_COEFF = TOTAL NUMBER OF TERMS IN SS ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_ELEMS = NUMBER OF ELEMENTS ! N_EL_FRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! N_GEOM = NUMBER OF GEOMETRY NODES ! N_HOMO = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS ! N_ITER = NO. OF ITERATIONS TO BE RUN (USUALLY 1) ! N_MAT = NUMBER OF MATERIAL TYPES ! NODES = NODAL INCIDENCES OF ALL ELEMENTS ! N_PARM = DIMENSION OF PARAMETRIC SPACE ! NP_FIX = INTEGER PROPERTIES AT ALL NODES ! N_QP = NUMBER OF QUADRATURE POINTS, >= LT_QP ! N_R_B = NUMBER OF ROWS IN B AND E MATRICES ! N_SPACE = DIMENSION OF SPACE ! N_FILE1 = UNIT FOR POST SOLUTION MATRICES STORAGE ! N_FILE2 = OPTIONAL UNITS FOR USER (USED WHEN > 0) ! N_FILE3,4 = OPTIONAL UNITS FOR USER (USED WHEN > 0) ! NUL_COL = 1, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO ! PT = QUADRATURE COORDINATES ! S = ELEMENT SQUARE MATRIX ! S_SKY = SYSTEM EQUATIONS SQUARE MATRIX ! STRAIN = STRAIN OR GRADIENT VECTOR ! STRAIN_0 = INITIAL STRAIN OR GRADIENT VECTOR ! STRESS = STRESS VECTOR ! WT = QUADRATURE WEIGHTS ! X = COORDINATES OF SYSTEM NODES ! XYZ = SPACE COORDINATES AT A POINT ! GENERATE ELEMENT EQUATIONS & POST SOLUTION MATRICES WRITE (N_PRT, "(/, 'BEGINNING ASSEMBLY')") IF ( I_BUG > 0 .AND. N_ITER > 1 ) PRINT *, 'ITERATION ', ITER IF ( U_FLUX > 0 ) REWIND U_FLUX IF ( N_FILE1 > 0 ) REWIND N_FILE1 IF ( N_FILE2 > 0 ) REWIND N_FILE2 COL_SUMS = 0.d0 LT = 1 ; LAST_LT = 0 ! INITIALIZE ELEMENT TYPE ! LOOP OVER ALL STANDARD ELEMENTS FIRST, ! THEN ANY OF THEIR BOUNDARY SEGMENT ELEMENTS WRITE (N_PRT, "('BEGINNING STANDARD ELEMENT ASSEMBLY')") DO IE = 1, L_S_TOT ! for elements and any boundary segments !b call at(125) IF ( DELETED (IE) ) CYCLE ! to skip deleted element CALL SET_THIS_ELEMENT_NUMBER (IE) IF ( N_MIXED > 0 .AND. IE == (N_ELEMS + 1) ) WRITE (N_PRT, & "('BEGINNING MIXED_BC SEGMENTS ASSEMBLY')") IF ( N_F_SEG > 0 .AND. IE == (N_ELEMS + N_MIXED + 1) ) & WRITE (N_PRT, "('BEGINNING FLUX BOUNDARY SEGMENTS ASSEMBLY')") !--> GET ELEMENT TYPE NUMBER IF ( N_L_TYPE > 1 ) LT = L_TYPE (IE) ! SAME AS LAST TYPE ? IF ( LT /= LAST_LT ) THEN ! this is a new type, or first element !b call at(137) ! GET CONTROLS FOR THIS TYPE, AND ALLOCATE ITS ARRAYS CALL SET_ELEM_TYPE_INFO (LT) ! including LAST_LT END IF ! a new element type ! INITIALIZE AND CHECK IF ACTIVE !b call list_type_alloc_status !b S = 0.d0 ; EL_M = 0.d0 ; C = 0.d0 !--> EXTRACT ELEMENT NODE NUMBERS ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES) !--> CALCULATE DEGREE OF FREEDOM NUMBERS INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES) if (debug_el_sq .AND. (N_ITER > 0 .OR. DYNAMIC .OR. TRANSIENT)) then !b call at(155) ; print *,'DD_OLD ' ; call rprint(DD_OLD, 1, N_D_FRE, 1) end if !--> GENERATE ELEMENT PROBLEM DEPENDENT MATRICES IF ( EXAMPLE > 0 ) THEN ! Select an old example CALL GENERATE_EL_MATRICES_NEW (IE, X, DD_OLD) IF ( DEBUG_EL_COL .OR. DEBUG_EL_SQ ) THEN PRINT *, 'Exiting GENERATE_EL_MATRICES_NEW eith C, S:' IF ( IE <= 3 ) CALL RPRINT (C, 1, LT_FREE, 1) !b IF ( IE <= 3 ) CALL RPRINT (S, LT_FREE, LT_FREE, 1) !b END IF ELSE CALL GENERATE_EL_MATRICES (IE, X, DD_OLD) END IF ! Select an old example !b call at(153) !--> SAVE ELEMENT MATRICES FOR ELEMENT REACTIONS IF ( L_REACT > 0 ) WRITE (L_REACT) S, C !--> STORE THE MATRICES IN SYSTEM EQUATIONS IF ( STATIC ) THEN IF ( I_BUG > 0 ) PRINT *,'BEGIN STATIC ANALYSIS' CALL SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, S_SKY, & S_SKY_LOWER) IF ( NUL_COL == 0) THEN ! SOME NON_ZERO ELEMENT SOURCES CALL STORE_COLUMN (LT_FREE, INDEX, C, CC) COL_SUMS = COL_SUMS + SUM ( C ) ! Error checking IF ( DEBUG_EL_COL .AND. IE <= 3 ) CALL RPRINT (C, 1, LT_FREE, 1) !b ELSE ! validate key words IF ( IE == 1 .AND. ANY (C /= 0.d0 ) ) THEN N_WARN = N_WARN + 1 PRINT *, 'WARNING: KEYWORD el_no_col IS WRONG' END IF END IF ! ELEMENT SOURCES ELSE IF ( TRANSIENT ) THEN ! BRANCH ON TIME INTEGRATION METHOD ! So element matrix S * D (t) + EL_M * D_dot (t) = C (t) ! Converts to algebraic K * D (t) = F (t) with ! K a function of S, EL_M, STEP_SIZE ! M a function of S, EL_M, STEP_SIZE, so that ! F is a function of M * D and C IF ( I_BUG > 0 ) PRINT *,'BEGIN TRANSIENT ANALYSIS' ! Convert element level arrays ! S (LT_FREE, LT_FREE), EL_M (LT_FREE, LT_FREE), ! C (LT_FREE), or DIAG_M (LT_FREE) with ! D (LT_FREE) to form ! One step transient forms ! K (LT_FREE, LT_FREE), M (LT_FREE, LT_FREE), ! F (LT_FREE) for assembly into ! K(S, EL_M, Dt) * D(t) = F(C, S, EL_M, Dt, D) CALL ONE_STEP_TRANSIENT_AT_EL (K, M, F) ! ASSEMBLE, ASSUMING CONSTANT C W.R.T. TIME CALL SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, K, S_SKY, & !b S_SKY_LOWER) CALL SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, M, M_SKY, & !b M_SKY_LOWER) CALL STORE_COLUMN (LT_FREE, INDEX, F, CC) !b ! ASSEMBLE REACTION ROWS FILE IF TRANSIENT REACTIONS DESIRED !b like SAVE_SKY_ROWS_GEN for use in GET_REACTIONS IF ( DEBUG_ASSEMBLY ) THEN print *,'K_SKY at 214 for element ', THIS_EL CALL SKY_FUL_SYM (S_SKY, I_DIAG, CC) print *,'M_SKY at 216' CALL SKY_FUL_SYM (M_SKY, I_DIAG, CC) END IF IF ( DEBUG_EL_COL ) THEN PRINT *,'EL, F, CC ', THIS_EL CALL RPRINT (F, 1, LT_FREE, 1) !b CALL RPRINT (CC, 1, N_D_FRE, 1) !b END IF ELSE IF ( DYNAMIC ) THEN IF ( I_BUG > 0 ) PRINT *,'BEGIN DYNAMIC ANALYSIS' ! ASSEMBLE, ASSUMING CONSTANT C W.R.T. TIME CALL SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, S_SKY, & !b S_SKY_LOWER) ! Symmetric mass matrix only CALL SKY_STORE_SQ_SYM (LT_FREE, INDEX, I_DIAG, EL_M, M_SKY) ! Symmetric element damping matrix only IF ( EL_DAMPING ) CALL SKY_STORE_SQ_SYM (LT_FREE, INDEX, & I_DIAG, DAMP, M_SKY_LOWER) CALL STORE_COLUMN (LT_FREE, INDEX, C, CC) !b ELSE IF ( EIGEN ) THEN !b PRINT *, 'BEGINNING EIGEN ANALYSIS' ! ASSEMBLE STIFFNESS AND MASS MATRICES FROM ELEMENT !b CALL SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, S_SKY, & !b S_SKY_LOWER) !b CALL SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, EL_M, M_SKY, & !b M_SKY_LOWER) STOP 'ASSEMBLE_SKYLINE_EQS: EIGEN ASSEMBLY UNDEFINED' ELSE ! UNKNOWN DEFAULT STOP 'ASSEMBLE_SKYLINE_EQS: ELEMENT ASSEMBLY CLASS UNDEFINED' END IF ! ANALYSIS CLASS !b IF ( .NOT. STATIC ) THEN ! NEED MASS ASSEMBLY !b IF ( DIAGONAL_MASS ) THEN !b CALL DIAGONALIZE_SQ_MATRIX (LT_FREE, EL_M, DIAG_M) !b CALL STORE_COLUMN (LT_FREE, INDEX, DIAG_M, M_SKY) !b ELSE !b CALL SKY_STORE_SQ_SYM (LT_FREE, INDEX, I_DIAG, EL_M, M_SKY) !b IF ( .NOT. SYMMETRIC ) THEN !b WRITE (N_PRT, *) 'WARNING: ' !b N_WARN = N_WARN + 1 !b END IF !b END IF ! DIAGONAL !b END IF ! NOT STATIC IF ( IE == N_ELEMS ) WRITE (N_PRT, & "('COMPLETED STANDARD ELEMENT ASSEMBLY')") IF ( N_MIXED > 0 .AND. IE == (N_ELEMS + N_MIXED) ) & WRITE (N_PRT, "('COMPLETED MIXED_BC SEGMENTS ASSEMBLY')") IF ( N_F_SEG > 0 .AND. & IE == (N_ELEMS + N_MIXED + N_F_SEG) ) WRITE (N_PRT, & "('COMPLETED FLUX BOUNDARY SEGMENTS ASSEMBLY')") IF ( IE == L_S_TOT .AND. (N_MIXED > 0 .OR. N_F_SEG > 0) ) & WRITE (N_PRT, "('ASSEMBLY COMPLETE.', /)") END DO ! over all elements and boundary segments IF ( NUL_COL == 0 .AND. COL_SUMS == 0.d0 ) THEN PRINT *, 'WARNING: ALL ELEMENT SOURCE VECTORS ARE NULL,' PRINT *, 'PUT el_no_col IN KEYWORD CONTROL' N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF IF ( DEBUG_EL_COL .AND. IE <= 3 ) CALL RPRINT (CC, 1, N_D_FRE, 1) !b WRITE (N_PRT, "(/, 'ASSEMBLY COMPLETED')") END SUBROUTINE ASSEMBLE_SKYLINE_EQS SUBROUTINE CEQ_SKYLINE (N_RES_EQ, NDEX_C_EQ, I_DOF_HI, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND SKYLINE STORAGE REQUIRED BY CONSTRAINT EQUATION ! MODIFICATION PROCEDURES. IT ALSO UPDATES I_DIAG AND N_COEFF ! AND I_DOF_HI. SOME OF THE CALCULATIONS ARE CONSERVATIVE ! (MORE STORAGE PROVIDED THAN NEEDED) ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE INTEGER, INTENT(INOUT) :: I_DOF_HI (N_D_FRE), I_DIAG (N_D_FRE) INTEGER, INTENT(IN) :: NDEX_C_EQ (MAX_ACT, N_CEQ) INTEGER, INTENT(IN) :: N_RES_EQ (MAX_ACT) INTEGER :: IC, IEQ, II, IIC, IICTOP, IJ, IJTOP, JJ INTEGER :: MINDEP, NR, N_R_BOT, NRTOP, NTEST, ILOOP ! MAX_ACT = NO ACTIVE CONSTR TYPES ! N_CEQ = TOTAL NO CONSTR EQS ! NDEX_C_EQ = (I,J) CONSTR DOF NO I OF EQ J ! N_D_FRE = TOTAL NO OF SYSTEM DEGREES OF FREEDOM IF ( I_BUG == 1 ) THEN WRITE (N_BUG, * ) "SKYLINE, I_DOF_HI: ", & (I_DOF_HI (II), II = 1, N_D_FRE) WRITE (N_BUG, * ) "I_DIAG: ", & (I_DIAG (II), II = 1, N_D_FRE) END IF IEQ = N_RES_EQ (1) !--> LOOP OVER NON DIAGONAL CONSTRAINTS DO IC = 2, MAX_ACT NTEST = N_RES_EQ (IC) IF ( NTEST > 0 ) THEN !--> * LOOP OVER TYPE IC EQUATIONS DO ILOOP = 1, NTEST IEQ = IEQ + 1 ! ! * THE MODIFICATIONS OF STORAGE ARE DONE FOR THE CURRENT ! * EQUATION: DJ + SUM ( AK * DK ) = E ! * TWO STEPS ARE NEEDED FOR THE STORAGE MODIFICATION ! ! * FIND THE NR DOF AND LOWEST "MINDEP" OF THE DEPENDENT DOFS NR = NDEX_C_EQ (1, IEQ) MINDEP = N_D_FRE DO JJ = 2, IC MINDEP = MIN (NDEX_C_EQ (JJ, IEQ), MINDEP) END DO ! ! * FIND THE TOP AND THE BOTTOM OF THE SKYLINE OF COLUMN NR CALL SKY_TOP_BOT (I_DIAG, NR, NRTOP, N_R_BOT) ! ! * STEP 1: HORIZONTAL FILL-IN (READ THEORY) ! DO IJ = NRTOP, N_R_BOT IJTOP = IJ - I_DOF_HI (IJ) + 1 ! * SEARCH FOR "PROBABLE NON ZEROS" IN COL NR & MODIFY TOP IF ( IJTOP <= NR .OR. IJ < NR ) THEN IJTOP = MIN (MINDEP, IJTOP) I_DOF_HI (IJ) = IJ - IJTOP + 1 END IF END DO ! ! * STEP 2: MODIFY THE CONSTRAINT DOFS COLS ! * THE TOP WILL BE THE MIN OF THREE NUMBERS ! DO II = 1, IC IIC = NDEX_C_EQ (II, IEQ) IICTOP = IIC - I_DOF_HI (IIC) + 1 IICTOP = MIN (IICTOP, MINDEP) IICTOP = MIN (IICTOP, NRTOP) I_DOF_HI (IIC) = IIC - IICTOP + 1 END DO ! ! * END OF MODIFICATIONS OF THE SKYLINE STORAGE ! * UPDATE I_DIAG AND N_COEFF ! I_DIAG = GET_SKY_DIAG (I_DOF_HI) N_COEFF = I_DIAG (N_D_FRE) ! IF ( I_BUG == 1 ) THEN WRITE (N_BUG, * ) "SKYLINE, I_DOF_HI: ", & (I_DOF_HI (II), II = 1, N_D_FRE) WRITE (N_BUG, * ) "I_DIAG: ", & (I_DIAG (II), II = 1, N_D_FRE) END IF END DO END IF END DO END SUBROUTINE CEQ_SKYLINE SUBROUTINE MESH_SKYLINE (NODES, I_DOF_HI, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! USE MESH COLUMN HEIGHTS TO FIND DIAGONAL COEFFICIENTS ! FOR SYMMETRIC SKYLINE STORAGE MODE, ASSUMING SYMMETRIC ! COLUMNS STORED FROM TOP DOWN ! * * * * * * * * * * * * * * * * * * * * * * * * * ! TOTAL NUMBER OF SQ MATRIX TERMS = I_DIAG(N_D_FRE) Use System_Constants Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER, INTENT(OUT) :: I_DIAG (N_D_FRE), I_DOF_HI (N_D_FRE) INTEGER :: J, COUNT, N_PT ! LOOP, INACTIVE COUNTER ! I_DOF_HI(I) = COLUMN HEIGHT OF EQ I, WITH DIAG ! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ NUMBER ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ! N_ELEMS = NUMBER OF ELEMENTS ! NODES = NODAL INCIDENCES OF ALL ELEMENTS CALL SKY_HI (NODES, I_DOF_HI) !b IF ( I_BUG > 0 ) PRINT *, 'I_DOF_HI, ', I_DOF_HI ! TEST FOR INACTIVE EQUATIONS TO BE CORRECTED COUNT = 0 DO J = 1, N_D_FRE IF ( I_DOF_HI (J) < 0 ) STOP 'NEGATIVE HEIGHTS, MESH_SKYLINE' IF ( I_DOF_HI (J) == 0 ) THEN ! INACTIVE EQ WITH NO ENTRIES IF ( COUNT == 0 ) THEN PRINT *,'WARNING: NULL EQUATION AT EQUATION NUMBER = ', J CALL GET_NODE_AT_DOF (J, N_G_DOF, N_PT) PRINT *,'NOTE: FIRST NULL EQUATION IS AT NODE ', N_PT N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF COUNT = COUNT + 1 I_DOF_HI (J) = 1 END IF END DO IF ( COUNT > 0 ) THEN PRINT *, 'WARNING: TOTAL NUMBER OF NULL EQS = ', COUNT N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF ! CONVERT TO STORAGE LOCATIONS I_DIAG = GET_SKY_DIAG (I_DOF_HI) END SUBROUTINE MESH_SKYLINE SUBROUTINE SKY_BC_GEN (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, S_SKY, & CC, I_DIAG, S_SKY_LOWER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! APPLY BOUNDARY CONSTRAINT EQUATIONS TO GENERAL SKYLINE SYSTEM ! OF EQS, S*D = CC ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: COEF_C_EQ (MAX_ACT, N_CEQ) INTEGER, INTENT(IN) :: N_RES_EQ (MAX_ACT), & NDEX_C_EQ (MAX_ACT, N_CEQ), & I_DIAG (N_D_FRE) REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF), CC (N_D_FRE) REAL(DP), INTENT(INOUT) :: S_SKY_LOWER (NOT_SYM) INTEGER :: IC, IEQ, NEQ, NTEST ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_COEFF = NUMBER OF NON-ZERO COEFFICIENTS IN S_SKY ! NOT_SYM = NUMBER OF NON-ZERO COEFF IN S_SKY_LOWER, ! = 0 OR N_COEFF ! S_SKY = SYS. EQ. SQ. MATRIX ! CC = SYS. EQ. COL. MATRIX ! N_G_DOF = NO. PARAMETERS PER NODE ! N_RES_EQ = NO OF CONSTRAINT EQS OF EACH TYPE ! COEF_C_EQ (I,J) = CONSTRAINT EQS COEFF I FOR EQ J ! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J ! N_CEQ = NUMBER OF CONSTRAINT EQUATIONS IF ( I_BUG == 1 ) THEN PRINT *, 'Entering SKY_BC_GEN' END IF ! WORK BACKWARD FROM LARGEST CONSTRAINT TYPE (DO TYPE 1 LAST) IF ( SUM(N_RES_EQ) < 1) THEN PRINT *, 'WARNING: NO RESTRAINTS APPLIED' N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF IEQ = SUM ( N_RES_EQ ) + 1 !b XXX if MAX_ACT > 1, split into 2 passes: LCE1, LCE2 DO IC = MAX_ACT, 1, -1 NTEST = N_RES_EQ (IC) IF ( NTEST > 0 ) THEN ! Apply constraint of type IC IF ( IC == 1 ) THEN !--> TYPE 1 D(L1) = C1 DO NEQ = 1, NTEST IEQ = IEQ - 1 CALL SKY_TYPE_1_GEN (NDEX_C_EQ (1, IEQ), & COEF_C_EQ (1, IEQ), S_SKY, CC, I_DIAG, S_SKY_LOWER) END DO ELSE IF ( .NOT. SYMMETRIC ) & STOP 'NOT SYMMETRIC, SKY_BC_GEN for TYPE > 1' !--> TYPE 2 D(L1) + C1*D(L2) = C2 !--> TYPE 3 D(L1) + C1*D(L2) + C2*D(L3) = C3, ETC. DO NEQ = 1, NTEST IEQ = IEQ - 1 !b CALL SKY_LCE (S_SKY, CC, IC, NDEX_C_EQ (1:IC, IEQ), & !b COEF_C_EQ (1:IC,IEQ), I_DIAG) CALL SKY_LCE_SYM (S_SKY, CC, IC, NDEX_C_EQ (1:IC, IEQ), & COEF_C_EQ (1:IC,IEQ), I_DIAG) END DO END IF END IF END DO END SUBROUTINE SKY_BC_GEN SUBROUTINE CHECK_SKY_SYS (S, C, I_DIAG, S_LOWER) !jea ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! WARNING, this code is not fully checked. No known bugs. ! CHECK SKYLINE SYSTEM FOR INVALID EQUATIONS & WARN ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE) REAL(DP), INTENT(INOUT) :: S_LOWER (NOT_SYM) !jea INTEGER :: I, J, K REAL(DP) :: SMAX, TEST, ZERO = 0.d0 LOGICAL :: WARN ! C = SYSTEM COLUMN MATRIX ! I_DIAG = POINTER TO DIAGONAL COEFFICIENT IN S ! S = SYSTEM SQUARE MATRIX IN SKY MODE (upper) ! S_LOWER = SYSTEM SQUARE MATRIX IN SKY MODE (lower) !jea ! N_COEFF = NUMBER OF COEFFICIENTS IN S ! NOT_SYM = NUMBER OF COEFFICIENTS IN S_LOWER !jea ! N_D_FRE = NUMBER OF EQUATIONS SMAX = ZERO; WARN = .FALSE. ! INITALIZE DO I = 1, N_D_FRE !b TEST = ABS (S (I_DIAG (I) ) ) TEST = S (I_DIAG (I) ) IF ( ABS (TEST) > SMAX) SMAX = TEST IF ( TEST <= 0.d0 ) WARN = .TRUE. END DO IF ( WARN ) THEN WRITE (N_PRT, "(/, 'WARNINGS FROM CHECK_SKY_SYS:')") N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF IF ( SMAX <= ZERO ) STOP & 'ALL ELEMENT STIFFNESSES ZERO, CHECK_SKY_SYS' K = 0 DO I = 1, MAX_NP DO J = 1, N_G_DOF K = K + 1 TEST = S (I_DIAG (K) ) IF ( TEST <= ZERO ) THEN IF ( TEST == ZERO ) THEN WRITE (N_PRT, 200) I, J 200 FORMAT ('WARNING, NODE ',I5,' DOF',I3,' RESTRAINED') N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF IF ( TEST < ZERO ) THEN WRITE (N_PRT, 300) I, J 300 FORMAT ('ERROR, NODE ',I5,' DOF',I3,' RESTRAINED') N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF ! SET DOF K TO ZERO IN SKYLINE MODE CALL SKY_TYPE_1_GEN (K, ZERO, S, C, I_DIAG, S_LOWER) !jea END IF END DO END DO END SUBROUTINE CHECK_SKY_SYS SUBROUTINE SKY_DIAG (I_DOF_HI, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! USE COLUMN HEIGHTS TO FIND DIAGONAL LOCATIONS ! FOR SKYLINE STORAGE MODE ! * * * * * * * * * * * * * * * * * * * * * * * * * ! ASSUMING SYMMETRIC COLUMNS STORED FROM TOP DOWN Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DOF_HI (N_D_FRE) INTEGER, INTENT(OUT) :: I_DIAG (N_D_FRE) INTEGER :: I, IPOINT ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ! I_DOF_HI(I) = COLUMN HEIGHT OF EQ I, WITH DIAG ! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ IN UPPER TRIANGLE ! TOTAL NUMBER OF SQ MATRIX TERMS = I_DIAG(N_D_FRE) IPOINT = 0 DO I = 1, N_D_FRE IPOINT = IPOINT + I_DOF_HI (I) I_DIAG (I) = IPOINT END DO END SUBROUTINE SKY_DIAG FUNCTION GET_SKY_DIAG (I_DOF_HI) RESULT (I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! USE COLUMN HEIGHTS TO FIND DIAGONAL LOCATIONS ! FOR SKYLINE STORAGE MODE ! * * * * * * * * * * * * * * * * * * * * * * * * * ! ASSUMING SYMMETRIC COLUMNS STORED FROM TOP DOWN Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DOF_HI (N_D_FRE) INTEGER :: I_DIAG (N_D_FRE) INTEGER :: I, IPOINT ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ! I_DOF_HI(I) = COLUMN HEIGHT OF EQ I, WITH DIAG ! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ IN UPPER TRIANGLE ! TOTAL NUMBER OF SQ MATRIX TERMS = I_DIAG(N_D_FRE) IPOINT = 0 DO I = 1, N_D_FRE IPOINT = IPOINT + I_DOF_HI (I) I_DIAG (I) = IPOINT END DO END FUNCTION GET_SKY_DIAG SUBROUTINE SKY_HI (NODES, I_DOF_HI) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND COLUMN HEIGHTS OF SYSTEM EQUATIONS IN SKYLINE ! STORAGE MODE (ZERO HEIGHT IS INACTIVE EQ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for ELEM_NODES (LT_N), INDEX (LT_FREE) Use Interface_Header Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER, INTENT(OUT) :: I_DOF_HI (N_D_FRE) ! Automatic Arrays INTEGER :: LT_NODES (NOD_PER_EL), LT_INDEX (N_EL_FRE) INTEGER :: L_HIGH (N_EL_FRE) INTEGER :: IE, J, NDX ! DELETED = LOGICAL ELEM FLAG, TRUE IF DELETED ! ELEM_NODES = ELEMENT NODAL INCIDENCES ! INDEX(I) = SYS DOF NUMBER OF ELEMENT DOF I ! I_DOF_HI(I) = COL HEIGHT OF SYS DOF I ! N_D_FRE = TOTAL NO OF SYSTEM DOF ! N_EL_FRE = MAX NUMBER OF ELEMENT PARAMETERS (DOF) ! N_ELEMS = NUMBER OF ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! NODES = NODAL INCIDENCES OF ALL ELEMENTS ! ZERO HEIGHTS I_DOF_HI = 0 ! LOOP OVER ELEMENTS DO IE = 1, N_ELEMS ! (boundary segments use same nodes) IF ( DELETED (IE) ) CYCLE ! to next element LT_NODES = 0; LT_INDEX = 0 !b ! EXTRACT NODES, FIND DOF NOS LT_N = LT_DATA (1, L_TYPE(IE)) !b CALL ... !b LT_NODES = GET_ELEM_NODES (IE, LT_N, NODES) LT_NODES (1:LT_N) = GET_ELEM_NODES (IE, LT_N, NODES) ! FIND ELEMENT COLUMN HEIGHTS LT_FREE = LT_N*N_G_DOF !b LT_INDEX = GET_ELEM_INDEX (LT_N, LT_NODES) LT_INDEX (1:LT_FREE) = GET_ELEM_INDEX (LT_N, LT_NODES) CALL EL_HIGH (LT_INDEX, L_HIGH) ! COMPARE WITH CURRENT MAXIMUMS DO J = 1, LT_FREE NDX = LT_INDEX (J) IF ( NDX > 0 ) THEN IF ( I_DOF_HI (NDX) < L_HIGH (J) ) & I_DOF_HI (NDX) = L_HIGH (J) END IF END DO ! for element dof END DO ! for all elements END SUBROUTINE SKY_HI FUNCTION GET_SKY_HEIGHTS (NODES) RESULT (I_DOF_HI) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND COLUMN HEIGHTS OF SYSTEM EQUATIONS IN SKYLINE ! STORAGE MODE (ZERO HEIGHT IS INACTIVE EQ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for ELEM_NODES (LT_N), INDEX (LT_FREE) Use Interface_Header Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER :: I_DOF_HI (N_D_FRE) ! Automatic Arrays INTEGER :: LT_NODES (NOD_PER_EL), LT_INDEX (N_EL_FRE) INTEGER :: L_HIGH (N_EL_FRE) INTEGER :: IE, J, NDX ! DELETED = LOGICAL ELEM FLAG, TRUE IF DELETED ! LT_NODES = ELEMENT NODAL INCIDENCES ! LT_INDEX(I) = SYS DOF NUMBER OF ELEMENT DOF I ! I_DOF_HI(I) = COL HEIGHT OF SYS DOF I ! N_D_FRE = TOTAL NO OF SYSTEM DOF ! N_EL_FRE = MAX NUMBER OF ELEMENT PARAMETERS (DOF) ! N_ELEMS = NUMBER OF ELEMENTS ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF PARAMETERS PER NODE ! NODES = NODAL INCIDENCES OF ALL ELEMENTS ! ZERO HEIGHTS I_DOF_HI = 0 ! LOOP OVER ELEMENTS DO IE = 1, N_ELEMS ! (boundary segments use same nodes) IF ( DELETED (IE) ) CYCLE ! to next element ! EXTRACT NODES, FIND DOF NOS LT_N = LT_DATA (1, L_TYPE(IE)) LT_NODES = GET_ELEM_NODES (IE, LT_N, NODES) ! FIND ELEMENT COLUMN HEIGHTS LT_FREE = LT_N*N_G_DOF LT_INDEX = GET_ELEM_INDEX (LT_N, LT_NODES) CALL EL_HIGH (LT_INDEX, L_HIGH) ! COMPARE WITH CURRENT MAXIMUMS DO J = 1, LT_FREE NDX = LT_INDEX (J) IF ( NDX > 0 ) THEN IF ( I_DOF_HI (NDX) < L_HIGH (J) ) & I_DOF_HI (NDX) = L_HIGH (J) END IF END DO ! for element dof END DO ! for all elements END FUNCTION GET_SKY_HEIGHTS SUBROUTINE SKY_LCE (S, C, NCD, NDX, A, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! APPLY LINEAR CONSTRAINTS TO SKYLINE SYMMETRIC EQUATIONS ! S*D = C, WITH S IN VECTOR MODE AND ! D(NDX(1))+A(1)*D(NDX(2))+...A(NCD-1)*D(NDX(NCD))=A(NCD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), PARAMETER :: ZERO = 0.d0, ONE = 1.d0 INTEGER, INTENT(IN) :: NCD, I_DIAG (N_D_FRE), NDX (NCD) REAL(DP), INTENT(IN) :: A (NCD) REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE) LOGICAL :: SKIP INTEGER :: I, IBOT, I_J_V, I_NR_V, ITOP, J, K, L, LL, NR REAL(DP) :: CR, E, ESC, SRR, SRRP1 ! ! N_COEFF = NUMBER OF TERMS FOR SQ MATRIX IN VECTOR MODE ! NCD = TOTAL NUMBER OF DOF IN CONSTRAINT EQUATION ! N_D_FRE = TOTAL NUMBER OF DEGREES OF FREEDOM ! NDX(I) = SYS DOF NOS OF CONSTRAINT TERM I ! NR = REDUNDANT DEGREE OF FREEDOM = NDX(1) ! A(J) = NORMALIZED COEFF OF (J+1) TERM, A0 = 1.0 ! S*D=C SYSTEM PARTITIONED ! C = SYSTEM COLUMN MATRIX | SII | SIR | SID | |DI| |CI| ! S = SYSTEM SQUARE MATRIX |-----|-----|-----| |--| |--| ! I-INDEPENDENT | SRI | SRR | SRD |*|DR|=|CR| ! R-REDUNDANT |-----|-----|-----| |--| |--| ! D-DEPENDENT | SDI | SDR | SDD | |DD| |CD| IF ( I_BUG > 0 ) PRINT *, 'Entering SKY_LCE' IF ( ANY (NDX < 1) ) STOP 'SKY_LCE: BAD CONSTRAINT INDEX' IF ( ANY (NDX > N_D_FRE) ) STOP 'SKY_LCE: BAD CONSTRAINT INDEX' E = A (NCD) NR = NDX (1) SRR = S (I_DIAG (NR) ) SRRP1 = SRR + ONE CR = C (NR) ESC = E * SRRP1 - CR ! FIND TOP AND BOTTOM OF THE SKYLINE OF THE NR COLUMN CALL SKY_TOP_BOT (I_DIAG, NR, ITOP, IBOT) ! FORM MODIFIED COLUMN MATRIX, CX = CX - E*SXR IF ( E /= ZERO ) THEN DO I = ITOP, IBOT !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V /= 0 ) C (I) = C (I) - E * S (I_NR_V) END DO END IF ! ADDITIONAL COLUMN CHANGES FOR CD AND CR IF ( NCD > 1 .AND. ESC /= ZERO ) THEN DO K = 2, NCD C (NDX (K) ) = C (NDX (K) ) + A (K - 1) * ESC END DO END IF C (NR) = E ! *** SQUARE MATRIX COLUMN MODIFICATIONS *** IF ( NCD > 1) THEN ! FORM SID, BEGIN SDD DO K = 2, NCD J = NDX (K) DO I = ITOP, IBOT !b IF ( I_BUG > 0 ) PRINT *, 'I, J ', I, J !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V) CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V /= 0 ) THEN IF ( S (I_NR_V) /= ZERO ) THEN IF ( I_J_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE' ! DO NOT MODIFY NR COLUMN, MAY NEED IT LATER SKIP = .FALSE. DO LL = K + 1, NCD IF ( I == NDX (LL) ) SKIP = .TRUE. END DO IF ( (.NOT. SKIP) .AND. I /= NR ) & S (I_J_V) = S (I_J_V) - S (I_NR_V) * A (K - 1) END IF ! S (I_NR_V) END IF ! I_NR_V END DO END DO ! COMPLETE SDD DO K = 2, NCD I = NDX (K) DO L = 2, K J = NDX (L) !b IF ( I_BUG > 0 ) PRINT *, 'I, J ', I, J !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V) CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_J_V == 0 .OR. I_NR_V == 0 ) & STOP 'INVALID ADDRESS, SKY_LCE' S (I_J_V) = S (I_J_V) + SRRP1 * A (K - 1) * A (L - 1) & - A (L - 1) * S (I_NR_V) END DO END DO END IF !b XXX end LCE1 begin LCE2, split into 2 routines ! *** INSERT CONSTRAINT EQUATION *** ! WARNING: NEXT LOOP NOT VALID FOR COUPLED CONSTRAINTS DO I = ITOP, IBOT !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V /= 0 ) S (I_NR_V) = ZERO END DO S (I_DIAG (NR)) = 1.d0 IF ( NCD > 1 ) THEN DO K = 2, NCD I = NDX (K) !b IF ( I_BUG > 0 ) PRINT *, 'I, NR ', I, NR CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE' S (I_NR_V) = A (K - 1) END DO ! *** MODIFICATIONS COMPLETED, CHECK DIAGONAL *** DO K = 2, NCD I = NDX (K) IF ( S (I_DIAG (I)) <= ZERO ) THEN ; WRITE (6, * ) & 'WARNING: NEGATIVE DIAGONAL FOR CONSTRAINT SET', NDX N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END DO END IF IF ( I_BUG > 0 ) PRINT *, 'Exiting SKY_LCE' END SUBROUTINE SKY_LCE SUBROUTINE SKY_SOLVE_SYM (A, B, X, I_DIAG, FACT, BACK) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PERFORM A =(U)T*D*U FACTORIZATION AND/OR BACKSUBSTITUTION ! OF SYMMETRIC POSITIVE DEFINITE SYSTEM OF EQUATIONS A*X = B ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_COEFF, N_D_FRE IMPLICIT NONE LOGICAL, INTENT(IN) :: FACT, BACK INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) REAL(DP), INTENT(INOUT) :: A (N_COEFF), B (N_D_FRE) REAL(DP), INTENT(OUT) :: X (N_D_FRE) INTEGER :: I, ID, IE, IH, IR, IS, J, JD, JH, JR, K REAL(DP) :: D, ENERGY, ZERO = 0.d0 ! ! A(N_COEFF) = UPPER TRIANGULAR COEFFICIENT MATRIX STORED IN ! COLUMN FORM (HOLDS D & U ON EXIT) ! B(N_D_FRE) = RIGHT SIDE VECTOR (HOLDS X ON EXIT) ! X(N_D_FRE) = SOLUTION VECTOR ! I_DIAG(N_D_FRE) = ADDRESSES OF DIAGONAL TERMS IN A(N_COEFF) ! N_D_FRE = NUMBER OF EQUATIONS ! FACT = .TRUE., FACTOR A(N_COEFF) ! .FALSE., DO NOT FACTOR A(N_COEFF) ! BACK = .TRUE., FORWARD REDUCE B(N_D_FRE) & BACKSUBSTITUTE ! .FALSE., DO NOT FORWARD REDUCE & BACKSUBSTITUTE ! !.... FACTOR A, REDUCE B ! JR = 0 DO J = 1, N_D_FRE JD = I_DIAG (J) JH = JD-JR IS = J - JH + 2 IF ( JH < 2 ) GO TO 600 IF ( JH == 2 ) GO TO 300 100 IF ( .NOT. FACT ) GO TO 500 IE = J - 1 K = JR + 2 ID = I_DIAG (IS - 1) ! !.... REDUCE ALL EQUATIONS EXCEPT DIAGONAL ! DO I = IS, IE IR = ID ID = I_DIAG (I) IH = MIN0 (ID-IR - 1, I - IS + 1) IF ( IH > 0 ) A (K) = A (K) & - DOT_PRODUCT (A (K-IH:K-1), A (ID-IH:ID-1)) K = K + 1 END DO ! on I ! !.... REDUCE DIAGONAL TERM ! 300 IF ( .NOT. FACT ) GO TO 500 IR = JR + 1 IE = JD-1 K = J - JD DO I = IR, IE ID = I_DIAG (K + I) IF ( A (ID) == ZERO ) THEN STOP 'ZERO DIAGONAL IN SKY_SOLVE_SYM' ELSE D = A (I) A (I) = A (I) / A (ID) A (JD) = A (JD) - D * A (I) END IF END DO ! on I ! !.... REDUCE RHS ! 500 IF ( BACK ) B (J) = B (J) & - DOT_PRODUCT(A (JR+1:JR+JH-1), B (IS-1:IS+JH-3)) 600 CONTINUE JR = JD END DO ! on J IF ( .NOT. BACK ) RETURN ! !.... DIVIDE BY DIAGONAL PIVOTS ! ENERGY = 0.d0 !b DO I = 1, N_D_FRE D = B (I) ID = I_DIAG (I) IF ( A (ID) /= ZERO ) B (I) = B (I) / A (ID) ENERGY = ENERGY + D * B (I) END DO !b PRINT *, ' ' !b PRINT *, 'EQUATION STRAIN ENERGY VALUE = ', ENERGY * 0.5d0 IF ( ENERGY < 0.d0 ) THEN PRINT *, 'WARNING, SKY_SOLVE_SYM: NEGATIVE ENERGY IN SYSTEM' PRINT *, 'WARNING: IS unsymmetric NEEDED AS KEYWORD INPUT ?' N_WARN = N_WARN + 2 PRINT *, 'ERROR, SKY_SOLVE_SYM: SYMMETRIC SOLVER FATAL ERROR' STOP 'SKY_SOLVE_SYM: FATAL ERROR, NEGATIVE ENERGY' ELSE !b PRINT *, 'EQUATION STRAIN ENERGY NORM = ', SQRT (ENERGY) END IF ! INVALID SYSTEM ! !.... BACKSUBSTITUTE ! J = N_D_FRE JD = I_DIAG (J) 800 D = B (J) !9 replace with while on J J = J - 1 IF ( J <= 0 ) RETURN JR = I_DIAG (J) IF ( JD-JR > 1 ) THEN IS = J - JD+JR + 2 K = JR - IS + 1 DO I = IS, J B (I) = B (I) - A (I + K) * D END DO X = B END IF JD = JR GO TO 800 END SUBROUTINE SKY_SOLVE_SYM SUBROUTINE SKY_FORWARD_SYM (A, B, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PREFORM FORWARD SUBSTITUTION, ONLY, FOR PREVIOUSLY ! FACTORED SYMMETRIC SKYLINE MATRIX SYSTEM, A*X = B VIA ! SKY_SOLVE_SYM (A, B, X, I_DIAG, .TRUE., .FALSE.) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_COEFF, N_D_FRE IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) REAL(DP), INTENT(IN) :: A (N_COEFF) REAL(DP), INTENT(INOUT) :: B (N_D_FRE) !b REAL(DP), INTENT(IN) :: X (N_D_FRE) INTEGER :: IS, J, JD, JH, JR ! A(N_COEFF) = UPPER TRIANGULAR COEFFICIENT MATRIX STORED IN ! COLUMN FORM (HOLDS D & U ON EXIT) ! B(N_D_FRE) = RIGHT SIDE VECTOR (HOLDS X ON EXIT) ! X(N_D_FRE) = SOLUTION VECTOR ! I_DIAG(N_D_FRE) = ADDRESSES OF DIAGONAL TERMS IN A(N_COEFF) ! N_D_FRE = NUMBER OF EQUATIONS JR = 0 DO J = 1, N_D_FRE JD = I_DIAG (J) JH = JD-JR IS = J - JH + 2 !b IF ( JH < 2 ) GO TO 600 IF ( JH > 1 ) THEN B (J) = B (J) - DOT_PRODUCT(A (JR+1:JR+JH-1), B (IS-1:IS+JH-3)) END IF !b 600 CONTINUE JR = JD END DO ! on J END SUBROUTINE SKY_FORWARD_SYM SUBROUTINE SKY_SOLVE_UNSYM (A, C, B, X, I_DIAG, FACT, BACK) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PROGRAM TO PERFORM (C\A)=L*D*U FACTORIZATION AND/OR BACKSUBSTITUTION ! OF AN UNSYMMETRIC SYSTEM OF EQUATIONS, (C\A)*X = B ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! for N_COEFF, N_D_FRE IMPLICIT NONE LOGICAL, INTENT(IN) :: FACT, BACK INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) REAL(DP), INTENT(INOUT) :: A (N_COEFF), C (N_COEFF), B (N_D_FRE) REAL(DP), INTENT(OUT) :: X (N_D_FRE) INTEGER :: I, ID, IR, J, JD, JH, JR, K, M, N_END, N_START, NT REAL(DP) :: D ! ! A(N_COEFF) = UPPER TRIANGULAR COEFFICIENT MATRIX STORED IN ! COLUMN FORM (HOLDS D & U ON EXIT) ! C(N_COEFF) = LOWER TRIANGULAR COEFFICIENT MATRIX STORED IN ! COLUMN FORM (HOLDS L ON EXIT), DIAGONAL IGNORED ! B(N_D_FRE) = RIGHT SIDE VECTOR (COULD HOLD X ON EXIT) ! X(N_D_FRE) = SOLUTION VECTOR ! I_DIAG(N_D_FRE) = ADDRESSES OF DIAGONAL TERMS IN A(N_COEFF) AND C ! N_COEFF = NUMBER OF STORED TERMS IN A AND C, = I_DIAG(N_D_FRE) ! N_D_FRE = NUMBER OF EQUATIONS ! FACT = .TRUE., FACTOR A(N_COEFF) AND C(N_COEFF) ! .FALSE., DO NOT FACTOR A(N_COEFF) AND C(N_COEFF) ! BACK = .TRUE., FORWARD REDUCE B(N_D_FRE) & BACKSUBSTITUTE ! .FALSE., DO NOT FORWARD REDUCE & BACKSUBSTITUTE ! JR = 0 DO 15 J = 1, N_D_FRE JD = I_DIAG (J) JH = JD-JR IF ( JH <= 1 ) GO TO 10 N_START = J + 1 - JH N_END = J - 1 IF ( .NOT. FACT ) GO TO 20 K = JR + 1 ID = 0 ! !.... REDUCE ALL EQUATIONS EXCEPT DIAGONAL ! DO 30 I = N_START, N_END IR = ID ID = I_DIAG (I) NT = MIN0 (ID-IR - 1, I - N_START) IF ( NT == 0 ) GO TO 40 A (K) = A (K) - DOT_PRODUCT (A (K-NT:K-1), C (ID-NT:ID-1)) C (K) = C (K) - DOT_PRODUCT (C (K-NT:K-1), A (ID-NT:ID-1)) 40 IF ( A (ID) /= 0. ) THEN C (K) = C (K) / A (ID) ELSE STOP 'ZERO DIAGONAL IN SKY_SOLVE_UNSYM' END IF K = K + 1 30 END DO ! !.... REDUCE DIAGONAL TERM ! A (JD) = A (JD) & - DOT_PRODUCT (A (JR+1:JR+JH-1), C (JR+1:JR+JH-1)) ! !.... FORWARD REDUCE THE R.H.S. ! 20 IF ( BACK ) B (J) = B (J) & - DOT_PRODUCT (C (JR+1:JR+JH-1), B (N_START:N_START+JH-2)) 10 JR = JD 15 END DO IF ( .NOT. BACK ) RETURN ! !.... BACKSUBSTITUTION ! J = N_D_FRE JD = I_DIAG (J) IF ( A (JD) /= 0. ) B (J) = B (J) / A (JD) IF ( N_D_FRE == 1 ) THEN X = B ; RETURN END IF 50 D = B (J) J = J - 1 JR = I_DIAG (J) IF ( JD-JR <= 1 ) GO TO 60 M = J - JD+JR + 2 K = JR - M + 1 DO 70 I = M, J B (I) = B (I) - A (I + K) * D 70 END DO 60 JD = JR IF ( A (JD) /= 0. ) B (J) = B (J) / A (JD) IF ( J > 1 ) GO TO 50 ! MOVE INTO X X = B END SUBROUTINE SKY_SOLVE_UNSYM SUBROUTINE SKY_STORE_SQ (LT_FREE, INDEX, I_DIAG, S, SS, SS_LOWER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! STORE ELEMENT SQUARE MATRIX TO SYSTEM SQUARE MATRIX STORED ! IN SYMMETRIC OR NON-SYMMETRIC SKYLINE MODE (Revised 5/9/98) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: LT_FREE INTEGER, INTENT(IN) :: INDEX (LT_FREE), I_DIAG (N_D_FRE) REAL(DP), INTENT(INOUT) :: SS (N_COEFF), SS_LOWER (NOT_SYM) REAL(DP), INTENT(IN) :: S (LT_FREE, LT_FREE) INTEGER :: I, J, NDX_I, NDX_J, NDX_V ! N_COEFF = NO COEFF IN UPPER HALF OF SQ MATRIX = I_DIAG(N_D_FRE) ! NOT_SYM = NO COEFF IN LOWER HALF OF SQ MATRIX = I_DIAG(N_D_FRE) ! N_D_FRE = TOTAL NO OF DOF IN SYSTEM ! LT_FREE = NUMBER OF ELEMENT DEGREES OF FREEDOM ! INDEX(I) = SYS DOF NO OF ELEMENT DOF I ! I_DIAG(I) = LOCATION OF DIAGONAL OF I-TH EQ ! S = ELEMENT SQUARE MATRIX ! SS = UPPER HALF OF SYS SQ MATRIX IN SKYLINE VECTOR MODE ! SS_LOWER = LOWER HALF OF SYS SQ MATRIX IN SKYLINE VECTOR MODE ! SYMMETRIC = .true. if symmetric so: ss_lower = 0, NOT_SYM = 0 ! LOOP OVER ELEMENT COEFFICIENTS DO J = 1, LT_FREE ! element columns NDX_J = INDEX (J) ! system column ! ALLOW FOR OMITTED NODES IF ( NDX_J > 0 ) THEN ! equation exists DO I = 1, LT_FREE ! element rows NDX_I = INDEX (I) ! system row IF ( NDX_I > 0 ) THEN ! equation exists IF ( I_BUG > 0 .AND. SYMMETRIC .AND. I == J & .AND. S (I, J) < 0.d0) THEN PRINT *,'BAD DIAGONAL ELEMENT TERM = ', S (I, J) PRINT *,'ELEMENT ROW = ', I PRINT *,'SYSTEM ROW = ', NDX_I STOP 'NEGATIVE ELEMENT DIAGONAL, SKY_STORE_SQ' END IF ! FIND SYSTEM COEFF IN VECTOR SS (AND SS_LOWER) NDX_V = I_DIAG ( MAX0 (NDX_I, NDX_J) ) - IABS (NDX_J-NDX_I) IF ( NDX_I <= NDX_J ) THEN ! Upper half SS (NDX_V) = SS (NDX_V) + S (I, J) ELSE ! Lower half IF ( .NOT. SYMMETRIC ) & SS_LOWER (NDX_V) = SS_LOWER (NDX_V) + S (I, J) END IF END IF END DO END IF END DO END SUBROUTINE SKY_STORE_SQ SUBROUTINE SKY_TOP_BOT (I_DIAG, J_COL, J_TOP, J_BOT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND THE TOP (J_TOP) AS WELL AS (MAXIMUM) BOTTOM ! (J_BOT) INDEX OF THE SKYLINE OF COL "J_COL" ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! ASSUMING SYMM EQS, COLS STORED FROM TOP DOWN Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE), J_COL INTEGER, INTENT(OUT) :: J_TOP, J_BOT INTEGER :: I, MINI LOGICAL :: DEBUG = .FALSE. ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ ! ! * FINDING J_TOP IF ( J_COL > 1 ) THEN J_TOP = J_COL - (I_DIAG (J_COL) - I_DIAG (J_COL - 1) ) + 1 ELSE J_TOP = 1 END IF ! * FINDING J_BOT J_BOT = J_COL DO I = J_COL + 1, N_D_FRE MINI = I - (I_DIAG (I) - I_DIAG (I - 1) ) + 1 IF ( MINI <= J_COL ) J_BOT = I END DO IF ( DEBUG ) WRITE (N_BUG, *) "SKY_TOP_BOT: ", J_COL, J_TOP, J_BOT END SUBROUTINE SKY_TOP_BOT SUBROUTINE SUB_IN_SKY (I_DIAG, I, J, I_J_V) ! * * * * * * * * * * * * * * * * * * * * * * * ! CONVERT (I,J) FULL SYMMETRIC MATRIX SUBSCRIPTS TO ! (I_J_V) SUBSCRIPT OF VECTOR SKYLINE STORAGE MODE. ! I_J_V = 0 IF OUTSIDE SKYLINE AND THUS HAS 0 VALUE. ! * * * * * * * * * * * * * * * * * * * * * * * ! ASSUMING SYMM EQS, COLS STORED FROM TOP DOWN ! REVISED FOR I_J_V = 0, 9/97 Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE), I, J INTEGER, INTENT(OUT) :: I_J_V INTEGER :: ID ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ ID = MAX0 (I,J) !b if ( i < 1 .or. I > N_D_FRE ) stop 'invalid diag. SUB_IN_SKY' !b if ( j < 1 .or. j > N_D_FRE ) stop 'invalid diag. SUB_IN_SKY' I_J_V = I_DIAG (ID) - IABS(I-J) ! test for out of bounds if (I_J_V < 1 .or. I_J_V > n_coeff ) then print *, 'out of bounds, SUB_IN_SKY: I, J, V = ', I, J, I_J_V CALL SKY_FUL_LOCATIONS (I_DIAG) ! locations stop 'out of bounds, SUB_IN_SKY' end if ! TEST FOR ZERO ENTRY OUTSIDE SKYLINE IF ( ID > 1 ) THEN IF ( I_J_V <= I_DIAG (ID - 1) ) I_J_V = 0 END IF END SUBROUTINE SUB_IN_SKY FUNCTION GET_SKY_SUBSCRIPT (I_DIAG, I, J) RESULT (I_J_V) ! * * * * * * * * * * * * * * * * * * * * * * * ! CONVERT (I,J) FULL SYMMETRIC MATRIX SUBSCRIPTS TO ! (I_J_V) SUBSCRIPT OF VECTOR SKYLINE STORAGE MODE. ! I_J_V = 0 IF OUTSIDE SKYLINE AND THUS HAS 0 VALUE. ! * * * * * * * * * * * * * * * * * * * * * * * ! ASSUMING SYMM EQS, COLS STORED FROM TOP DOWN ! REVISED FOR I_J_V = 0, 9/97 Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE), I, J INTEGER :: I_J_V INTEGER :: ID ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ! I_DIAG(I) = LOCATION OF DIAG OF I-TH EQ ID = MAX0 (I,J) IF ( I < 1 .OR. I > N_D_FRE ) STOP 'INVALID INPUT GET_SKY_SUBSCRIPT' IF ( J < 1 .OR. J > N_D_FRE ) STOP 'INVALID INPUT GET_SKY_SUBSCRIPT' I_J_V = I_DIAG (ID) - IABS(I-J) ! test for out of bounds if (I_J_V < 1 .or. I_J_V > n_coeff ) then print *, 'out of bounds, GET_SKY_SUBSCRIPT: I, J, V = ', I, J, I_J_V CALL SKY_FUL_LOCATIONS (I_DIAG) ! locations stop 'out of bounds, GET_SKY_SUBSCRIPT' end if ! TEST FOR ZERO ENTRY OUTSIDE SKYLINE IF ( ID > 1 ) THEN IF ( I_J_V <= I_DIAG (ID - 1) ) I_J_V = 0 END IF END FUNCTION GET_SKY_SUBSCRIPT SUBROUTINE CHECK_SKY_SYS_SYM (S, C, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! WARNING, this code is not fully checked. No known bugs. ! CHECK SKYLINE SYSTEM FOR INVALID EQUATIONS & WARN ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE) INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) INTEGER :: I, J, K REAL(DP) :: SMAX, TEST, ZERO = 0.d0 ! C = SYSTEM COLUMN MATRIX ! I_DIAG = POINTER TO DIAGONAL COEFFICIENT IN S ! S = SYSTEM SQUARE MATRIX IN SKYLINE MODE ! N_COEFF = NUMBER OF COEFFICIENTS IN S ! N_D_FRE = NUMBER OF EQUATIONS SMAX = ZERO DO I = 1, N_D_FRE TEST = ABS (S (I_DIAG (I) ) ) IF ( TEST > SMAX) SMAX = TEST END DO IF ( SMAX <= ZERO ) STOP & 'ALL ELEMENT STIFFNESSES ZERO, CHECK_SKY_SYS' K = 0 DO I = 1, MAX_NP DO J = 1, N_G_DOF K = K + 1 TEST = S (I_DIAG (K) ) IF ( TEST <= ZERO ) THEN IF ( TEST == ZERO ) THEN ; WRITE (N_PRT, 2) I, J 2 FORMAT ('WARNING, NODE ',I5,' DOF',I3,' WAS RESTRAINED') N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF IF ( TEST < ZERO ) THEN WRITE (N_PRT, 3) I, J 3 FORMAT ('ERROR, NODE ',I5,' DOF',I3,' WAS RESTRAINED') STOP 'ERROR, NEGATIVE DIAGONAL. CHECK_SKY_SYS_SYM' END IF ! negative diagonal ! SET DOF K TO ZERO IN SKYLINE MODE CALL SKY_TYPE_1_SYM (K, ZERO, S, C, I_DIAG) ENDIF END DO END DO END SUBROUTINE CHECK_SKY_SYS_SYM SUBROUTINE SAVE_SKY_ROWS_SYM (I_BC, S_SKY, CC, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * !--> SAVE INDEPENDENT REACTION EQUATIONS FROM SKYLINE MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: S_SKY (N_COEFF), CC (N_D_FRE) INTEGER, INTENT(IN) :: I_BC (MAX_NP), I_DIAG (N_D_FRE) ! Automatic Arrays INTEGER :: INDEX (N_G_DOF), KODES (N_G_DOF) INTEGER :: J, JKV, J1, J2, IG, INDX, K ! N_G_DOF = NUMBER OF DOF PER NODE ! N_REACT = SEQUENTIAL UNIT TO STORE REACTION DATA ! N_D_FRE = TOTAL NUMBER OF EQUATIONS IF ( N_REACT > 0 ) THEN ; REWIND (N_REACT) ELSE ; STOP 'NO UNIT FOR REACTION STORAGE, SAVE_SKY_ROWS_SYM' END IF IF ( I_BUG > 0 ) PRINT *, 'Entering SAVE_SKY_ROWS_SYM' ! LOOP OVER EQUATIONS FOR ESSENTIAL BC FLAG DO J = 1, MAX_NP IF ( I_BC (J) > 0 ) THEN CALL PT_CODE (J, I_BC (J), KODES) INDEX = GET_INDEX_AT_PT (J) DO IG = 1, N_G_DOF IF ( KODES (IG) > 0 ) THEN INDX = INDEX (IG) ! COLUMN NUMBERS J1 = INDX J2 = INDX IF (INDX > 1) J1 = J2 - (I_DIAG(INDX) - I_DIAG(INDX-1)) + 1 IF ( INDX < N_D_FRE ) THEN ! FIND LAST COLUMN DO K = N_D_FRE, INDX, -1 CALL SUB_IN_SKY (I_DIAG, INDX, K, JKV) IF ( JKV > 0 ) THEN ! FOUND NON-ZERO ENTRY J2 = K ; EXIT ! THIS LOOP END IF END DO END IF ! WRITE NODE, PARAMETER, RANGE OF NON-ZEROS if ( j1 < 1 ) then !b j1 = 1 !b N_WARN = N_WARN + 1 print *, 'WARNING: fixed bounds error, SAVE_SKY_ROWS_SYM' !b end if !b if ( j2 > N_D_FRE ) then !b j2 = N_D_FRE !b N_WARN = N_WARN + 1 print *, 'fWARNING: ixed bounds error, SAVE_SKY_ROWS_SYM' !b end if !b WRITE (N_REACT) J, IG, J1, J2 ! SAVE ROW OF EQUILIBRIUM EQ DO K = J1, J2 CALL SUB_IN_SKY (I_DIAG, INDX, K, JKV) IF ( JKV > 0 ) THEN WRITE (N_REACT) S_SKY (JKV) ELSE WRITE (N_REACT) 0.d0 ENDIF END DO WRITE (N_REACT) CC (INDX) ENDIF END DO ENDIF END DO END SUBROUTINE SAVE_SKY_ROWS_SYM SUBROUTINE SKY_BC_SYM (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, S_SKY, & CC, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! APPLY BOUNDARY CONSTRAINT EQUATIONS TO SYMMETRIC ! SKYLINE SYSTEM OF EQS, S_SKY*D = CC ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: N_RES_EQ (MAX_ACT) INTEGER, INTENT(IN) :: NDEX_C_EQ (MAX_ACT, N_CEQ) INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) REAL(DP), INTENT(IN) :: COEF_C_EQ (MAX_ACT, N_CEQ) REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF), CC (N_D_FRE) INTEGER :: IC, IEQ, NEQ, NTEST ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_COEFF = NUMBER OF NON-ZERO COEFFICIENTS IN S_SKY ! S_SKY = SYS. EQ. SQ. MATRIX ! CC = SYS. EQ. COL. MATRIX ! N_G_DOF = NO. PARAMETERS PER NODE ! N_RES_EQ = NO OF CONSTRAINT EQS OF EACH TYPE ! COEF_C_EQ (I,J) = CONSTRAINT EQS COEFF I FOR EQ J ! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J ! N_CEQ = NUMBER OF CONSTRAINT EQUATIONS IF ( I_BUG > 0 ) PRINT *, 'Entering SKY_BC_SYM' IEQ = 0 !b XXX if MAX_ACT > 1, split into 2 passes: LCE1, LCE2 ! DO TYPE ONE LAST DO IC = MAX_ACT, 1, -1 NTEST = N_RES_EQ (IC) IF ( NTEST > 0 ) THEN IF ( IC == 1 ) THEN !--> TYPE 1 D(L1) = C1 DO NEQ = 1, NTEST IEQ = IEQ + 1 CALL SKY_TYPE_1_SYM (NDEX_C_EQ (1, IEQ), & COEF_C_EQ (1, IEQ), S_SKY, CC, I_DIAG) END DO ELSE !--> TYPE 2 D(L1) + C1*D(L2) = C2 !--> TYPE 3 D(L1) + C1*D(L2) + C2*D(L3) = C3, etc DO NEQ = 1, NTEST IEQ = IEQ + 1 CALL SKY_LCE_SYM (S_SKY, CC, IC, NDEX_C_EQ (1:IC, IEQ), & COEF_C_EQ (1:IC, IEQ), I_DIAG) END DO ENDIF ! SPECIAL TYPE 1, OR GENERAL TYPE ENDIF ! TYPE ANY TYPE IC ARE ACTIVE END DO END SUBROUTINE SKY_BC_SYM SUBROUTINE SKY_LCE_SYM (S, C, NCD, NDX, A, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! APPLY LINEAR CONSTRAINTS TO SKYLINE SYMMETRIC EQUATIONS ! S*D = C, WITH S IN VECTOR MODE AND ! D(NDX(1))+A(1)*D(NDX(2))+...A(NCD-1)*D(NDX(NCD))=A(NCD) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), PARAMETER :: ZERO = 0.d0, ONE = 1.d0 INTEGER, INTENT(IN) :: NCD, I_DIAG (N_D_FRE), NDX (NCD) REAL(DP), INTENT(IN) :: A (NCD) REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE) LOGICAL :: SKIP INTEGER :: I, IBOT, I_J_V, I_NR_V, ITOP, J, K, L, LL, NR REAL(DP) :: CR, E, ESC, SRR, SRRP1 ! N_COEFF = NUMBER OF TERMS FOR SQ MATRIX IN VECTOR MODE ! NCD = TOTAL NUMBER OF DOF IN CONSTRAINT EQUATION ! N_D_FRE = TOTAL NUMBER OF DEGREES OF FREEDOM ! NDX(I) = SYS DOF NOS OF CONSTRAINT TERM I ! NR = REDUNDANT DEGREE OF FREEDOM = NDX(1) ! A(J) = NORMALIZED COEFF OF (J+1) TERM, A0 = 1.0 ! S*D=C SYSTEM PARTITIONED ! C = SYSTEM COLUMN MATRIX | SII | SIR | SID | |DI| |CI| ! S = SYSTEM SQUARE MATRIX |-----|-----|-----| |--| |--| ! I-INDEPENDENT | SRI | SRR | SRD |*|DR|=|CR| ! R-REDUNDANT |-----|-----|-----| |--| |--| ! D-DEPENDENT | SDI | SDR | SDD | |DD| |CD| IF ( I_BUG > 0 ) PRINT *, 'Entering SKY_LCE_SYM' IF ( ANY (NDX < 1) ) STOP 'SKY_LCE_SYM: BAD CONSTRAINT INDEX' IF ( ANY (NDX > N_D_FRE) ) STOP 'SKY_LCE_SYM: BAD CONSTRAINT INDEX' E = A (NCD) ! INITIALIZE NR = NDX (1) SRR = S (I_DIAG (NR) ) SRRP1 = SRR + ONE CR = C (NR) ESC = E * SRRP1 - CR ! FIND TOP AND BOTTOM OF THE SKYLINE OF THE NR COLUMN CALL SKY_TOP_BOT (I_DIAG, NR, ITOP, IBOT) ! FORM MODIFIED COLUMN MATRIX, CX = CX - E*SXR IF ( E /= ZERO ) THEN DO I = ITOP, IBOT CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V /= 0 ) C (I) = C (I) - E * S (I_NR_V) END DO ENDIF ! ADDITIONAL COLUMN CHANGES FOR CD AND CR IF ( NCD > 1 .AND. ESC /= ZERO ) THEN DO K = 2, NCD C (NDX (K) ) = C (NDX (K) ) + A (K - 1) * ESC END DO ENDIF C (NR) = E !b add averaging data ! *** SQUARE MATRIX COLUMN MODIFICATIONS *** IF ( NCD > 1) THEN ! FORM SID, BEGIN SDD DO K = 2, NCD J = NDX (K) DO I = ITOP, IBOT CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V) CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V /= 0 ) THEN IF ( S (I_NR_V) /= ZERO ) THEN IF ( I_J_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE_SYM' ! DO NOT MODIFY NR COLUMN, MAY NEED IT LATER SKIP = .FALSE. DO LL = K + 1, NCD IF ( I == NDX (LL) ) SKIP = .TRUE. END DO IF ( (.NOT. SKIP) .AND. I /= NR ) & S (I_J_V) = S (I_J_V) - S (I_NR_V) * A (K - 1) END IF ! S (I_NR_V) END IF ! I_NR_V END DO END DO ! COMPLETE SDD DO K = 2, NCD I = NDX (K) DO L = 2, K J = NDX (L) CALL SUB_IN_SKY (I_DIAG, I, J, I_J_V) CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_J_V == 0 .OR. I_NR_V == 0 ) & STOP 'INVALID ADDRESS, SKY_LCE_SYM' S (I_J_V) = S (I_J_V) + SRRP1 * A (K - 1) * A (L - 1) & - A (L - 1) * S (I_NR_V) END DO END DO ENDIF !b XXX end LCE1 begin LCE2, split into 2 routines ! *** INSERT CONSTRAINT EQUATION *** ! WARNING: NEXT LOOP NOT VALID FOR COUPLED CONSTRAINTS DO I = ITOP, IBOT CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V /= 0 ) S (I_NR_V) = ZERO END DO S (I_DIAG (NR)) = ONE !b add averaging data IF ( NCD > 1 ) THEN DO K = 2, NCD I = NDX (K) CALL SUB_IN_SKY (I_DIAG, I, NR, I_NR_V) IF ( I_NR_V == 0 ) STOP 'INVALID ADDRESS, SKY_LCE_SYM' S (I_NR_V) = A (K - 1) !b add averaging data END DO ! *** MODIFICATIONS COMPLETED, CHECK DIAGONAL *** DO K = 2, NCD I = NDX (K) IF ( S (I_DIAG (I)) <= ZERO ) THEN ; WRITE (6, * ) & 'WARNING: NEGATIVE DIAGONAL FOR CONSTRAINT SET', NDX N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END DO ENDIF IF ( I_BUG > 0 ) PRINT *, 'Leaving SKY_LCE_SYM' END SUBROUTINE SKY_LCE_SYM SUBROUTINE SKY_TYPE_1_SYM (LOC, VALUE, S, C, I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! APPLY TYPE 1 MODIFICATION TO SYMMERTIC SKYLINE EQS ! S*D = C, D(LOC) = VALUE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE), LOC REAL(DP), INTENT(IN) :: VALUE REAL(DP), INTENT(INOUT) :: S (N_COEFF), C (N_D_FRE) INTEGER :: I, ID, INV, ITOP ! N_D_FRE = DFREE NUMBER OF EQUATIONS ! N_COEFF = NUMBER OF NON-ZERO COEFFICIENTS IN S ! LOC = DOF NUMBER OF CONSTRAINED PARAMETER ! VALUE = GIVEN VALUE OF DOF NUMBER LOC ! S = SYSTEM SQUARE MATRIX IN SKYLINE STORAGE MODE ! C = FULL COLUMN MATRIX ! I_DIAG = POINTER TO DIAGONAL COEFFICIENT IN S ! IF ( LOC < 1 .OR. LOC > N_D_FRE ) STOP 'BAD LOC, SKY_TYPE_1_SYM' ! SUBTRACT COLUMN*VALUE FROM RHS DO I = 1, N_D_FRE ! FIND S(I, LOC) IN S VECTOR ID = MAX0 (I, LOC) INV = I_DIAG (ID) - IABS (I - LOC) ITOP = 1 IF (ID > 1) ITOP = I_DIAG (ID-1) + 1 ! IS IT OUTSIDE SKYLINE AND THUS ZERO? IF (INV >= ITOP) THEN C (I) = C (I) - VALUE * S (INV) S (INV) = 0.d0 ENDIF END DO ! RESET THE EQUATION ROW S (I_DIAG (LOC)) = 1.d0 C (LOC) = VALUE END SUBROUTINE SKY_TYPE_1_SYM SUBROUTINE SKY_FUL_SYM (S, I_DIAG, C) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! COPY SKYLINE S INTO SYMMETRIC FULL S AND PRINT WITH C ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: S (N_COEFF), C (N_D_FRE) INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! Automatic Arrays REAL(DP) :: FULL (N_D_FRE, N_D_FRE) INTEGER :: I, ID, IJV, ITOP, J ! N_COEFF = NO COEFF IN SQ MATRIX = I_DIAG(N_D_FRE) ! N_D_FRE = TOTAL NO OF DOF IN SYSTEM FULL = 0.d0 DO I = 1, N_D_FRE DO J = I, N_D_FRE !b FULL (I, J) = 0.d0 ! FIND S(I, J) IN S VECTOR ID = MAX0 (I, J) IJV = I_DIAG (ID) - IABS (I - J) ITOP = 1 IF ( ID > 1 ) ITOP = I_DIAG (ID-1) + 1 ! IS IT OUTSIDE SKYLINE AND THUS ZERO? IF ( IJV >= ITOP ) THEN FULL (I, J) = S (IJV) FULL (J, I) = FULL (I, J) END IF END DO END DO PRINT *,'SKYLINE COEFF VALUES FOLLOW' CALL RPRINT (FULL, N_D_FRE, N_D_FRE, 1) PRINT *,'COLUMN COEFF VALUES FOLLOW' CALL RPRINT (C, N_D_FRE, 1, 1) END SUBROUTINE SKY_FUL_SYM SUBROUTINE SKY_STORE_SQ_SYM (LTFREE, INDEX, I_DIAG, S, SS) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! STORE ELEMENT SQUARE MATRIX TO SYSTEM SQUARE ! MATRIX STORED IN SYMMETRIC SKYLINE MODE ! * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE INTEGER, INTENT(IN) :: LTFREE INTEGER, INTENT(IN) :: INDEX (LTFREE), I_DIAG (N_D_FRE) REAL(DP), INTENT(IN) :: S (LTFREE, LTFREE) REAL(DP), INTENT(INOUT) :: SS (N_COEFF) INTEGER :: I, J, JTEMP, NDX_I, NDX_J, NDX_V ! N_COEFF = NO COEFF IN SQ MATRIX = I_DIAG(N_D_FRE) ! N_D_FRE = TOTAL NO OF DOF IN SYSTEM ! LTFREE = NUMBER OF ELEMENT DEGREES OF FREEDOM ! INDEX(I) = SYS DOF NO OF ELEMENT DOF I ! I_DIAG(I) = LOCATION OF DIAGONAL OF I-TH EQ ! S = ELEMENT SQUARE MATRIX ! SS = SYS SQ MATRIX IN SKYLINE VECTOR MODE ! ! LOOP OVER ELEMENT COEFFICIENTS DO J = 1, LTFREE NDX_J = INDEX (J) ! ALLOW FOR OMITTED NODES IF ( NDX_J > 0 ) THEN JTEMP = I_DIAG (NDX_J) - NDX_J DO I = 1, LTFREE NDX_I = INDEX (I) IF ( NDX_I <= NDX_J .AND. NDX_I > 0 ) THEN IF ( (I == J) .AND. (S (I, J) < 0.d0)) THEN PRINT *,'BAD DIAGONAL ELEMENT TERM = ', S (I, J) PRINT *,'ELEMENT ROW = ', I PRINT *,'SYSTEM ROW = ', NDX_I STOP 'NEGATIVE ELEMENT DIAGONAL, SKY_STORE_SQ_SYM' END IF ! FIND SYSTEM COEFF IN VECTOR S NDX_V = JTEMP + NDX_I ! NDX_V = I_DIAG ( AMAX0 (NDX_I, NDX_J) ) ! - IABS (NDX_J-NDX_I) SS (NDX_V) = SS (NDX_V) + S (I, J) ENDIF END DO ENDIF END DO END SUBROUTINE SKY_STORE_SQ_SYM SUBROUTINE SAVE_SKY_ROWS_GEN (I_BC, S_SKY, CC, I_DIAG, S_SKY_LOWER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! SAVE INDEPENDENT REACTION EQUATIONS FROM A GENERAL SKYLINE ROW ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: S_SKY (N_COEFF), CC (N_D_FRE) ! always used REAL(DP), INTENT(IN) :: S_SKY_LOWER (NOT_SYM) ! may be used INTEGER, INTENT(IN) :: I_BC (MAX_NP), I_DIAG (N_D_FRE) ! always used ! Automatic Arrays INTEGER :: INDEX (N_G_DOF) ! equation numbers at a node INTEGER :: KODES (N_G_DOF) ! bc flags at a node INTEGER :: J ! node loop INTEGER :: IG ! dof loop at node INTEGER :: ROW ! equation number of dof at node INTEGER :: IN_SKY ! location in skyline, 0 if absent INTEGER :: COL, COL_START, COL_STOP ! column search at row ! N_COEFF = NUMBER OF TERMS IN UPPER STIFFNESS ! N_G_DOF = NUMBER OF DOF PER NODE ! N_REACT = SEQUENTIAL UNIT TO STORE REACTION DATA ! N_D_FRE = TOTAL NUMBER OF EQUATIONS ! NOT_SYM = NUMBER OF TERMS IN LOWER STIFFNESS, 0 OR N_COEFF IF ( N_REACT > 0 ) THEN ; REWIND (N_REACT) ELSE ; STOP 'NO UNIT FOR REACTION STORAGE' END IF ! LOOP OVER EQUATIONS FOR ESSENTIAL BC FLAG DO J = 1, MAX_NP ! node loop IF ( I_BC (J) > 0 ) THEN ! bc applied CALL PT_CODE (J, I_BC (J), KODES) ! get kodes INDEX = GET_INDEX_AT_PT (J) ! get dof at node DO IG = 1, N_G_DOF ! dof loop IF ( KODES (IG) < 0 ) THEN ! component IG has bc ! XXX > ROW = INDEX (IG) ! extract row number ! FIND COLUMN NUMBER BOUNDS FOR NON-ZERO COEFFICIENTS IN ROW COL_START = FIRST_COL_AT_SKY_ROW (ROW, I_DIAG) COL_STOP = LAST_COL_AT_SKY_ROW (ROW, I_DIAG) ! CORRECT FOR ANY INITIALLY NULL ROWS (GRAB DIAGONAL) IF ( COL_START > ROW ) COL_START = ROW IF ( COL_STOP < ROW ) COL_STOP = ROW ! WRITE NODE, PARAMETER, RANGE OF NON-ZEROS TO STORAGE WRITE (N_REACT) J, IG, COL_START, COL_STOP ! SAVE ROW OF NON_ZEROS IN SYSTEM EQUATIONS DO COL = COL_START, COL_STOP ! column loop IN_SKY = LOCATION_IN_SKY (ROW, COL, I_DIAG) ! locate IF ( IN_SKY > 0 ) THEN ! non-zero ? IF ( ROW <= COL ) THEN ! in upper WRITE (N_REACT) S_SKY (IN_SKY) ! upper value ELSE ! ROW > COL, ! in lower IF ( SYMMETRIC ) THEN ! symmetric ? WRITE (N_REACT) S_SKY (IN_SKY) ! upper value ELSE ! non_sym WRITE (N_REACT) S_SKY_LOWER (IN_SKY) ! lower value END IF END IF ELSE WRITE (N_REACT) 0.d0 ! zero term END IF END DO WRITE (N_REACT) CC (ROW) ! save row END IF END DO END IF END DO END SUBROUTINE SAVE_SKY_ROWS_GEN FUNCTION FIRST_COL_AT_SKY_ROW (ROW, I_DIAG) RESULT (COLUMN) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND FIRST COLUMN IN A ROW IN A SKYLINE STORAGE ! WARNING: CAN OCCUR AFTER A DIAGONAL ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! (ALSO FIRST ROW IN A COLUMN IN A SKYLINE STORAGE) Use system_constants ! N_D_FRE IMPLICIT NONE INTEGER, INTENT(IN) :: ROW ! row (or column) number INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! skyline diagonals INTEGER :: COLUMN ! first column ( or row) INTEGER :: COL ! loop INTEGER :: LOCATION_IN_SKY ! N_D_FRE = TOTAL NUMBER OF EQUATIONS COLUMN = 0 ! may be null DO COL = 1, N_D_FRE ! forward column loop IF ( LOCATION_IN_SKY (ROW, COL, I_DIAG) > 0 ) THEN ! non-zero COLUMN = COL ; EXIT ! search over END IF END DO IF ( COLUMN == 0 ) THEN PRINT *,'WARNING, NULL FIRST_COL_AT_SKY_ROW' N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END FUNCTION FIRST_COL_AT_SKY_ROW FUNCTION LAST_COL_AT_SKY_ROW (ROW, I_DIAG) RESULT (COLUMN) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND LAST COLUMN IN A ROW IN A SKYLINE STORAGE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! (ALSO LAST ROW IN A COLUMN IN A SKYLINE STORAGE) Use system_constants ! N_D_FRE ! Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: ROW ! row (or column) number INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! skyline diagonals INTEGER :: COLUMN ! first column ( or row) INTEGER :: COL ! loop INTEGER :: LOCATION_IN_SKY ! N_D_FRE = TOTAL NUMBER OF EQUATIONS COLUMN = 0 ! maybe a null row DO COL = N_D_FRE, 1, -1 ! backward column loop IF ( LOCATION_IN_SKY (ROW, COL, I_DIAG) > 0 ) THEN ! non-zero COLUMN = COL ; EXIT ! search over END IF END DO IF ( COLUMN == 0 ) THEN PRINT *,'WARNING, NULL LAST_COL_AT_SKY_ROW' N_WARN = N_WARN + 1 ! INCREMENT WARNING END IF END FUNCTION LAST_COL_AT_SKY_ROW FUNCTION FIRST_IN_SKY_ROW_OR_COL (NUMBER, I_DIAG) RESULT (IN_SKY) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND VECTOR LOCATION OF FIRST NON-ZERO COEFFICIENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use system_constants ! N_D_FRE Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: NUMBER ! row or column number INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! skyline diagonals INTEGER :: IN_SKY ! location in skyline INTEGER :: COL, LOCATION ! loop ! N_D_FRE = TOTAL NUMBER OF EQUATIONS IN_SKY = I_DIAG(1) ! IF NUMBER = 1, maybe 0 IF ( NUMBER > 1 ) THEN ! general case IN_SKY = 0 ! maybe null IF ( I_DIAG(NUMBER) > I_DIAG(NUMBER-1) ) THEN ! not null IN_SKY = I_DIAG(NUMBER-1) + 1 ! search over ELSE ! zero diagonal DO COL = NUMBER, N_D_FRE ! check columns LOCATION = LOCATION_IN_SKY (NUMBER, COL, I_DIAG) IF ( LOCATION > 0 ) THEN ! possible non-zero IN_SKY = LOCATION ; EXIT ! search over END IF END DO END IF END IF END FUNCTION FIRST_IN_SKY_ROW_OR_COL FUNCTION LAST_IN_SKY_ROW_OR_COL (NUMBER, I_DIAG) RESULT (IN_SKY) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND VECTOR LOCATION OF LAST NON-ZERO COEFFICIENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use system_constants ! N_D_FRE Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: NUMBER ! row or column number INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! skyline diagonals INTEGER :: IN_SKY ! location in skyline INTEGER :: COLUMN, LOCATION ! looping data ! N_D_FRE = TOTAL NUMBER OF EQUATIONS IN_SKY = I_DIAG(NUMBER) ! upper limit, maybe 0 DO COLUMN = N_D_FRE, 1, -1 ! backward column loop LOCATION = LOCATION_IN_SKY (NUMBER, COLUMN, I_DIAG) IF ( LOCATION > 0 ) THEN ! non-zero possible IN_SKY = LOCATION ; EXIT ! search over END IF END DO END FUNCTION LAST_IN_SKY_ROW_OR_COL FUNCTION LOCATION_IN_SKY (I, J, I_DIAG) RESULT (IN_SKY) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! CONVERT (I,J) FULL SYMMETRIC MATRIX SUBSCRIPTS TO ! (IN_SKY) SUBSCRIPT OF VECTOR SKYLINE STORAGE MODE. ! IN_SKY = 0 IF OUTSIDE SKYLINE AND THUS HAS 0 VALUE. ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ASSUMING COLUMNS STORED FROM TOP DOWN Use system_constants ! N_D_FRE IMPLICIT NONE INTEGER, INTENT(IN) :: I, J ! full subscripts INTEGER, INTENT(IN) :: I_DIAG(N_D_FRE) ! diagonal locations INTEGER :: IN_SKY ! location in skyline INTEGER :: ID ! max equation ! N_D_FRE = TOTAL NO OF SYSTEM EQUATIONS ID = MAX0 (I,J) ! max equation IF ( ID < 1 .OR. ID > N_D_FRE) STOP 'BAD ID, LOCATION_IN_SKY' IN_SKY = I_DIAG(ID) - IABS(I-J) ! possible non-zero ! TEST FOR ZERO ENTRY OUTSIDE SKYLINE IF ( ID > 1 ) THEN IF ( IN_SKY <= I_DIAG(ID - 1) ) IN_SKY = 0 END IF END FUNCTION LOCATION_IN_SKY SUBROUTINE SKY_FUL_GEN (S, I_DIAG, C, S_LOWER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! COPY SKYLINE VALUES INTO FULL FORM AND PRINT WITH C ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, SYMMETRIC Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: S (N_COEFF), C (N_D_FRE) REAL(DP), INTENT(IN) :: S_LOWER (NOT_SYM) INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! Automatic Arrays REAL(DP) :: FULL (N_D_FRE, N_D_FRE) INTEGER :: I, IJV, J ! N_D_FRE = TOTAL NUMBER OF EQUATIONS FULL = 0. DO I = 1, N_D_FRE ! row loop DO J = I, N_D_FRE ! (upper) column loop ! FIND S(I, J) IN S VECTOR IJV = LOCATION_IN_SKY (I, J, I_DIAG) ! position ! IS IT OUTSIDE SKYLINE AND THUS ZERO? IF ( IJV > 0 ) THEN ! stored ? FULL (I, J) = S (IJV) ! copy upper value IF ( I /= J ) THEN ! in lower column IF ( SYMMETRIC ) THEN ! symmetric ? FULL (J, I) = S (IJV) ! copy upper value ELSE FULL (J, I) = S_LOWER (IJV) ! copy lower value END IF END IF END IF END DO END DO PRINT *,'SKYLINE COEFF VALUES FOLLOW' CALL RPRINT (FULL, N_D_FRE, N_D_FRE, 1) ! print full system PRINT *,'COLUMN COEFF VALUES FOLLOW' CALL RPRINT (C, N_D_FRE, 1, 1) END SUBROUTINE SKY_FUL_GEN SUBROUTINE SKY_FUL_LOCATIONS (I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! PRINT SKYLINE STORAGE LOCATIONS IN FULL MATRIX FORM ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! Automatic Arrays INTEGER :: FULL (N_D_FRE, N_D_FRE) INTEGER :: I, IJV, J ! N_D_FRE = TOTAL NUMBER OF EQUATIONS DO I = 1, N_D_FRE ! row loop DO J = I, N_D_FRE ! upper col loop ! FIND S(I, J) LOCATION IN S VECTOR(S) IJV = LOCATION_IN_SKY (I, J, I_DIAG) ! may be 0 FULL (I, J) = IJV ; FULL (J, I) = IJV ! set location IF ( I == J ) FULL (J, J) = -IJV ! flag diagonal END DO END DO PRINT *,'SKYLINE COEFF LOCATIONS FOLLOW (NEG ON DIAGONAL)' CALL IPRINT (FULL, N_D_FRE, N_D_FRE) ! print locations END SUBROUTINE SKY_FUL_LOCATIONS SUBROUTINE SHOW_SKY_FUL (I_DIAG) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! SHOW SKYLINE AS A FULL CHARACTER MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! Automatic Arrays CHARACTER(LEN=1) :: FULL (N_D_FRE, N_D_FRE) ! work area CHARACTER(LEN=1) :: TEN(10) ! work area INTEGER :: I, I_J_V, J ! loops INTEGER :: FAKE ! avoid warnings ! N_D_FRE = TOTAL NUMBER OF EQUATIONS FAKE = N_D_FRE ! avoid warnings TEN = (/ '1', '2', '3', '4', '5', '6', '7', '8', '9', ':' /) PRINT *, 'SKYLINE IN FULL MATRIX:' select case (N_D_FRE) case (1:9); print *, ' ', ten(1:FAKE) case (10:19); print *, ' ', ten, ten(1:(FAKE-10)) case (20:29); print *, ' ', ten, ten, ten(1:(FAKE-20)) case (30:39); print *, ' ', ten, ten, ten, ten(1:(FAKE-30)) case default; print *, ' ', ten, ten, ten, ten, ten end select FULL = 'o' ! zero symbol DO I = 1, N_D_FRE ! row loop DO J = I, N_D_FRE ! upper column in row ! FIND S(I, J) IN S VECTOR I_J_V = LOCATION_IN_SKY (I, J, I_DIAG) ! may be 0 ! IS IT OUTSIDE SKYLINE AND THUS ZERO? IF ( I_J_V > 0 ) THEN ! is in skyline FULL (I, J) = 'U' ; FULL (J, I) = 'L' ! upper or lower END IF IF ( I == J ) FULL (J,J) = '\\' ! single \ most systems END DO END DO DO I = 1, N_D_FRE ! row J = I ; IF ( J > 9 ) J = MOD(J,10) ! row id PRINT *, J, ' ', FULL(I,:) ! id and columns in row END DO END SUBROUTINE SHOW_SKY_FUL SUBROUTINE SKY_TYPE_1_GEN (LOC, VALUE, S, C, I_DIAG, S_LOWER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! APPLY TYPE 1 MODIFICATION TO SYMMERTIC OR UNSYMMETRIC ! SKYLINE EQS, S*D = C, D(LOC) = VALUE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, SYMMETRIC Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE), LOC REAL(DP), INTENT(IN) :: VALUE REAL(DP), INTENT(INOUT) :: S (N_COEFF), S_LOWER (NOT_SYM) REAL(DP), INTENT(INOUT) :: C (N_D_FRE) REAL(DP) :: AVE INTEGER :: I, INV ! N_D_FRE = DFREE NUMBER OF EQUATIONS ! N_COEFF = NUMBER OF NON-ZERO COEFFICIENTS IN S ! NOT_SYM = NUMBER OF NON-ZERO COEFFICIENTS IN S_LOWER ! LOC = DOF NUMBER OF CONSTRAINED PARAMETER ! VALUE = GIVEN VALUE OF DOF NUMBER LOC ! S = UPPER SYSTEM MATRIX IN SKYLINE STORAGE MODE ! S_LOWER = LOWER SYSTEM MATRIX IN SKYLINE STORAGE MODE ! C = FULL COLUMN MATRIX ! I_DIAG = POINTER TO DIAGONAL COEFFICIENT IN S and S_LOWER ! SYMMETRIC = true if symmetric ! AVE = AVERAGE DIAGONAL VALUE IN CONSTRAINT ! IF ( LOC < 1 .OR. LOC > N_D_FRE ) THEN PRINT *,'WARNING: IMPOSSIBLE (FATAL) CONSTRAINT DOF = ', LOC PRINT *,'ASSOCIATED CONSTRAINT VALUE WAS = ', VALUE N_WARN = N_WARN + 1 STOP 'BAD LOC, SKY_TYPE_1_GEN' END IF AVE = S (I_DIAG(LOC)) ! scaling if ( ave < 1.e-6) ave = 1.d0 !b ! SUBTRACT THE LOC_TH COLUMN*VALUE FROM RHS ! AND ZERO THE LOC_TH ROW DO I = 1, N_D_FRE ! run down the rows ! FIND S(I, LOC) IN S VECTOR INV = LOCATION_IN_SKY (I, LOC, I_DIAG) ! IS IT OUTSIDE SKYLINE AND THUS ZERO? IF ( INV > 0 ) THEN ! may not be zero, get it IF ( SYMMETRIC ) THEN ! ALL IN UPPER C (I) = C (I) - VALUE * S (INV) ! subtract from rhs S (INV) = 0.d0 ! zero row and column ELSE ! may be in upper or lower IF ( I <= LOC ) THEN ! use upper C (I) = C (I) - VALUE * S (INV) ! subtract from rhs S (INV) = 0.d0 ! zero column S_LOWER (INV) = 0.d0 ! zero row ELSE ! use lower C (I) = C (I) - VALUE * S_LOWER (INV) ! subtract from rhs S_LOWER (INV) = 0.d0 ! zero column S (INV) = 0.d0 ! zero row END IF END IF END IF END DO ! RESET THE EQUATION ROW, scaled IF ( .NOT. SYMMETRIC ) S_LOWER (I_DIAG(LOC)) = 0.d0 S (I_DIAG (LOC)) = AVE C (LOC) = VALUE*AVE END SUBROUTINE SKY_TYPE_1_GEN SUBROUTINE EL_HIGH (LT_INDEX, L_HIGH) ! * * * * * * * * * * * * * * * * * * * * * * * * ! FIND SYSTEM COLUMN HEIGHTS OF AN ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * Use Elem_Type_Data ! for LT_INDEX (LT_FREE) IMPLICIT NONE INTEGER, INTENT(IN) :: LT_INDEX (LT_FREE) INTEGER, INTENT(OUT) :: L_HIGH (LT_FREE) INTEGER :: I, MIN, NDX ! LT_FREE = NO OF DEGREES OF FREEDOM OF ELEMENT ! LT_INDEX = SYSTEM DOF NUMBERS OF ELEMENT PARAMETERS ! L_HIGH(I) = COLUMN HEIGHT FOR EQUATION LT_INDEX(I) MIN = LT_INDEX (1) L_HIGH = 0 ! FIND MINIMUM LT_INDEX DO I = 1, LT_FREE NDX = LT_INDEX (I) ! ALLOW FOR OMITTED NODES IF ( NDX > 0 .AND. NDX < MIN ) MIN = NDX END DO ! CONVERT TO COLUMN HEIGHTS MIN = MIN - 1 DO I = 1, LT_FREE NDX = LT_INDEX (I) IF ( NDX > 0 ) L_HIGH (I) = NDX - MIN END DO END SUBROUTINE EL_HIGH SUBROUTINE VECTOR_MULT_SKY_GEN (S_SKY, VECTOR, I_DIAG, S_SKY_LOWER, & RESULT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! MULTIPLY A GENERAL SKYLINE ARRAY BY A VECTOR, SKY*VECTOR = RESULT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: S_SKY (N_COEFF) ! always used REAL(DP), INTENT(IN) :: VECTOR (N_D_FRE) ! always used REAL(DP), INTENT(IN) :: S_SKY_LOWER (NOT_SYM) ! may be used INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! always used REAL(DP), INTENT(OUT) :: RESULT (N_D_FRE) ! always used ! Automatic Arrays INTEGER :: ROW ! equation number of dof at node INTEGER :: IN_SKY ! location in skyline, 0 if absent INTEGER :: COL, COL_START, COL_STOP ! column search at row REAL (DP) :: DOT ! row dot product ! I_DIAG = DIAGONAL LOCATION IN SKYLINE VECTOR ! N_COEFF = NUMBER OF TERMS IN UPPER STIFFNESS ! N_G_DOF = NUMBER OF DOF PER NODE ! N_D_FRE = TOTAL NUMBER OF EQUATIONS ! NOT_SYM = NUMBER OF TERMS IN LOWER STIFFNESS, 0 OR N_COEFF ! S_SKY = UPPER HALF OF SS MATRIX IN SKYLINE MODE ! S_SKY_LOWER= LOWER HALF OF SS, IF NOT SYMMETRIC ! NOTE: Alternate way is to loop over non_zeros in skyline and ! multiply by the corresponding vector DO ROW = 1, N_D_FRE ! LOOP OVER EQUATIONS ! FIND COLUMN NUMBER BOUNDS FOR NON-ZERO COEFFICIENTS IN ROW COL_START = FIRST_COL_AT_SKY_ROW (ROW, I_DIAG) COL_STOP = LAST_COL_AT_SKY_ROW (ROW, I_DIAG) ! CORRECT FOR ANY INITIALLY NULL ROWS (GRAB DIAGONAL) IF ( COL_START > ROW ) COL_START = ROW IF ( COL_STOP < ROW ) COL_STOP = ROW ! MULTIPLY ROW OF NON_ZEROS IN SYSTEM EQUATIONS DOT = 0.d0 DO COL = COL_START, COL_STOP ! column loop IN_SKY = LOCATION_IN_SKY (ROW, COL, I_DIAG) ! locate IF ( IN_SKY > 0 ) THEN ! probably non-zero IF ( ROW <= COL ) THEN ! in upper storage DOT = DOT + S_SKY (IN_SKY) * VECTOR (COL) ELSE ! ROW > COL, ! in lower storage IF ( SYMMETRIC ) THEN ! symmetric upper value DOT = DOT + S_SKY (IN_SKY) * VECTOR (COL) ELSE ! non_sym lower value DOT = DOT + S_SKY_LOWER (IN_SKY) * VECTOR (COL) END IF END IF END IF END DO ! OVER COLUMNS RESULT (ROW) = DOT END DO ! OVER ROWS END SUBROUTINE VECTOR_MULT_SKY_GEN SUBROUTINE VECTOR_MULT_SKY_SYM (S_SKY, VECTOR, I_DIAG, RESULT) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! MULTIPLY A SYMMETRIC SKYLINE ARRAY BY A VECTOR, SKY*VECTOR = RESULT ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! N_D_FRE, N_COEFF, NOT_SYM, MAX_NP, SYMMETRIC Use Interface_Header IMPLICIT NONE REAL(DP), INTENT(IN) :: S_SKY (N_COEFF) ! always used REAL(DP), INTENT(IN) :: VECTOR (N_D_FRE) ! always used INTEGER, INTENT(IN) :: I_DIAG (N_D_FRE) ! always used REAL(DP), INTENT(OUT) :: RESULT (N_D_FRE) ! always used ! Automatic Arrays INTEGER :: ROW ! equation number of dof at node INTEGER :: IN_SKY ! location in skyline, 0 if absent INTEGER :: COL, COL_START, COL_STOP ! column search at row REAL (DP) :: DOT ! row dot product ! I_DIAG = DIAGONAL LOCATION IN SKYLINE VECTOR ! N_COEFF = NUMBER OF TERMS IN UPPER STIFFNESS ! N_G_DOF = NUMBER OF DOF PER NODE ! N_D_FRE = TOTAL NUMBER OF EQUATIONS ! NOT_SYM = NUMBER OF TERMS IN LOWER STIFFNESS, 0 OR N_COEFF ! S_SKY = UPPER HALF OF SS MATRIX IN SKYLINE MODE ! NOTE: Alternate way is to loop over non_zeros in skyline and ! multiply by the corresponding vector DO ROW = 1, N_D_FRE ! LOOP OVER EQUATIONS ! FIND COLUMN NUMBER BOUNDS FOR NON-ZERO COEFFICIENTS IN ROW COL_START = FIRST_COL_AT_SKY_ROW (ROW, I_DIAG) COL_STOP = LAST_COL_AT_SKY_ROW (ROW, I_DIAG) ! CORRECT FOR ANY INITIALLY NULL ROWS (GRAB DIAGONAL) IF ( COL_START > ROW ) COL_START = ROW IF ( COL_STOP < ROW ) COL_STOP = ROW ! MULTIPLY ROW OF NON_ZEROS IN SYSTEM EQUATIONS DOT = 0.d0 DO COL = COL_START, COL_STOP ! column loop IN_SKY = LOCATION_IN_SKY (ROW, COL, I_DIAG) ! locate IF ( IN_SKY > 0 ) DOT = DOT + S_SKY (IN_SKY) * VECTOR (COL) END DO ! OVER COLUMNS RESULT (ROW) = DOT END DO ! OVER ROWS END SUBROUTINE VECTOR_MULT_SKY_SYM SUBROUTINE TYPE_1_COL_ZERO (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ, CC) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ZERO COLUMN AT ROWS WITH EBC, FOR TRANSIENT OR DYNAMIC STEPS ! OF EQS, S*D = CC ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: COEF_C_EQ (MAX_ACT, N_CEQ) INTEGER, INTENT(IN) :: N_RES_EQ (MAX_ACT), & NDEX_C_EQ (MAX_ACT, N_CEQ) REAL(DP), INTENT(INOUT) :: CC (N_D_FRE) INTEGER :: IEQ, NEQ, NTEST ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! CC = SYS. EQ. COL. MATRIX ! N_G_DOF = NO. PARAMETERS PER NODE ! N_RES_EQ = NO OF CONSTRAINT EQS OF EACH TYPE ! COEF_C_EQ (I,J) = CONSTRAINT EQS COEFF I FOR EQ J ! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J ! N_CEQ = NUMBER OF CONSTRAINT EQUATIONS IF ( I_BUG == 1 ) THEN PRINT *, 'Entering TYPE_1_COL_ZERO' END IF ! WORK BACKWARD FROM LARGEST CONSTRAINT TYPE (DO TYPE 1 LAST) IEQ = SUM ( N_RES_EQ ) + 1 NTEST = N_RES_EQ (1) IF ( NTEST > 0 ) THEN ! Apply constraint of type IC !--> TYPE 1 D(L1) = C1 DO NEQ = 1, NTEST IEQ = IEQ - 1 CC (NDEX_C_EQ (1, IEQ)) = 0.d0 END DO END IF END SUBROUTINE TYPE_1_COL_ZERO