! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE JACOBI_MODEL (X, I_BC, NODES, N_RES_EQ, COEF_C_EQ, & NDEX_C_EQ, MASSES) !b NDEX_C_EQ, CC_IN) ! ****************************************************** ! * MODULAR PROGRAMS FOR FINITE ELEMENT ANALYSES * ! * VERSION 1, GENERAL EIGEN ANALYSIS SOLUTION * ! * COPYRIGHT J. E. AKIN, 2004 * ! ****************************************************** ! Fortran 90 version, uses double precision and "automatic arrays" ! ! A SET OF BUILDING BLOCK PROGRAMS ! BY ! DR. J. E. AKIN, P.E. ! DEPT. OF MECHANICAL ENGINEERING & MATERIALS SCIENCE ! RICE UNIVERSITY ! HOUSTON, TEXAS 77251-1892 ! ! email: akin@rice.edu ! Use System_Constants Use Elem_Type_Data Use Sys_Properties_Data Use GEOMETRIC_PROPERTIES Use Interface_Header !b IMPLICIT NONE LOGICAL :: FACT, BACK ! Given Arrays REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP), INTENT(IN) :: COEF_C_EQ (MAX_ACT, N_CEQ) !b REAL(DP), INTENT(IN) :: MASSES (N_D_FRE) !b ! point masses INTEGER, INTENT(IN) :: I_BC (MAX_NP) INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL) INTEGER, INTENT(IN) :: N_RES_EQ (MAX_TYP) INTEGER, INTENT(IN) :: NDEX_C_EQ (MAX_ACT, N_CEQ) !b ! Allocatable Arrays INTEGER, ALLOCATABLE :: L_TO_L_NEIGH (:, :) ! neighbors INTEGER, ALLOCATABLE :: L_TO_N_NEIGH (:, :) ! neighbors ! Allocatable Arrays, JACOBI only REAL(DP), ALLOCATABLE :: M_FULL (:, :) ! full mass REAL(DP), ALLOCATABLE :: S_FULL (:, :) ! full stiffness ! make automatic ?? ! Automatic Arrays REAL(DP) :: CC (N_D_FRE) REAL(DP) :: DD (N_D_FRE), DD_OLD (N_D_FRE) REAL(DP) :: SCP_AVERAGES (MAX_NP, SCP_FIT) REAL(DP) :: EL_ENORM (N_ELEMS), ERR_MAX REAL(DP) :: EL_REFINE (N_ELEMS) !b REAL(DP) :: C_MD (0) ! not in static analysis INTEGER :: L_FIRST (MAX_NP), L_LAST (MAX_NP) INTEGER :: L_TO_L_SUM (N_ELEMS) INTEGER :: L_TO_N_SUM (MAX_NP) !b INTEGER :: INDEX (N_G_DOF) !b passed by Interface_Header INTEGER :: LIST (N_G_DOF) !b passed by Interface_Header INTEGER, SAVE :: BAR_DOF = 1 ! limit bar chart loop INTEGER :: ITER, IBAR, IPARM, IE, OLD_DEL = -1 INTEGER :: IO_1, I, I_WARN, J, K REAL(DP) :: R_TEST, RATIO, OLD_NORM, NEW_NORM, DIFF_NORM ! JACOBI METHOD ONLY INTEGER :: NTYPE, EQ_N ! constraints INTEGER :: NS_MAX, NS_USED ! Jacobi sweeps INTEGER :: EIGEN_ORDER (N_D_FRE) ! eigen sort order !b INTEGER :: EIGEN_LOC (N_D_FRE) ! max abs location REAL(DP) :: VALUE ! fake eigenvalue at BC REAL(DP) :: EIGEN_SCALE (N_D_FRE) ! normalization REAL(DP) :: EIGEN_VAL (N_D_FRE) ! eigenvalues REAL(DP) :: EIGEN_VEC (N_D_FRE, N_D_FRE) ! eigenvectors by column ! SEE NOTATION.F FOR VARIABLE AND ARRAY NOTATIONS NS_MAX = SWEEPS ! keyword jacobi_sweeps (default 15) CALL ALLOCATE_GEOM_PROP ! GENERATE & LIST MESH NEIGHBORS (PATCH) INFORMATION IF ( SCP_NEIGH_FACE ) CALL SET_L_SHAPE_FACE_NEEDS ! for NEEDS !s I_WARN = 1 CALL FIRST_LAST_L_AT_PT (NODES, L_FIRST, L_LAST, I_WARN) !b NEEDS = 1 CALL COUNT_ELEMS_AT_ELEM (N_ELEMS, NOD_PER_EL, MAX_NP, L_FIRST, & L_LAST, NODES, NEEDS, L_TO_L_SUM, N_WARN) NEIGH_L = MAXVAL (L_TO_L_SUM) NEIGH_P = NEIGH_L ! IF ELEMENT BASED PATCHES !b ALLOCATE ( L_TO_L_NEIGH (NEIGH_L, N_ELEMS) ) CALL FORM_ELEMS_AT_EL (N_ELEMS, NOD_PER_EL, MAX_NP, & L_FIRST, L_LAST, NODES, N_SPACE, & L_TO_L_SUM, L_TO_L_NEIGH, & NEIGH_L, NEEDS, ON_BOUNDARY, U_ON_B) !b !b L_TO_L_SUM, L_TO_L_NEIGH, NEIGH_L) !xxx !b force proper on_boundary codes !xxx if ( NEEDS == 0 ) then !xxx OPEN (67, FILE='on_boundary_input', & !b !xxx ACTION='READ', STATUS='OLD', IOSTAT = IO_1) !b !xxx IF ( IO_1 > 0 ) THEN !xxx PRINT *,'OPENED FAILED FOR on_boundary_input' !xxx ELSE !xxx ON_BOUNDARY = .FALSE. !xxx DO ! FOREVER !xxx READ (67, *, IOSTAT = IO_1) IE !xxx IF ( IO_1 /= 0 ) EXIT ! THIS LOOP !xxx ON_BOUNDARY (IE) = .TRUE. !xxx END DO ! FOR FACE FLAGS !xxx END IF ! boundary file provided !xxx end if ! NEEDS !xxx !b MAX_FACES = 2 !b to avoid warning from compiler !!b test facing neighbor option ------------------------------- !! SET DATA FOR FACING PATCHES (MAKE INTO ROUTINE) !SELECT CASE (L_SHAPE) ! SHAPE OF ELEMENT ! CASE (1) ; NEEDS = 1 ; MAX_FACES = 2 ! LINE ELEMENT ! CASE (2) ; NEEDS = 2 ; MAX_FACES = 3 ! TRIANGULAR ELEMENT ! CASE (3) ; NEEDS = 2 ; MAX_FACES = 4 ! QUADRILATERAL ELEMENT ! CASE (4) ; NEEDS = 4 ; MAX_FACES = 6 ! HEXAHEDRA ELEMENT ! CASE (5) ; NEEDS = 3 ; MAX_FACES = 4 ! TETRAHEDRA ELEMENT ! CASE (6) ; NEEDS = 3 ; MAX_FACES = 5 ! WEDGE ELEMENT ! CASE DEFAULT ; NEEDS = 1 ; MAX_FACES = 6 !END SELECT !NEIGH_L = MAX_FACES !b override to facing neighbors !deallocate (L_TO_L_NEIGH) !allocate ( L_TO_L_NEIGH (MAX_FACES, N_ELEMS)) !L_TO_L_NEIGH = 0 !call FORM_FACING_ELEMENTS (N_ELEMS, NOD_PER_EL, MAX_NP, & ! L_FIRST, L_LAST, NODES, NEEDS,& ! MAX_FACES, L_TO_L_NEIGH) !!b end test facing neighbor option ------------------------------- IF ( LIST_EL_TO_EL ) CALL LIST_ELEMENTS_AT_EL (L_TO_L_NEIGH) CALL COUNT_L_AT_NODE (N_ELEMS, NOD_PER_EL, MAX_NP, NODES, & L_TO_N_SUM) IF ( I_BUG > 0 .OR. PT_EL_SUMS .OR. PT_EL_LIST ) THEN !b IF ( NEIGH_L > 0 ) PRINT *, 'L_TO_L Sparseness Value ', & !b float( count(L_TO_L_NEIGH == 0)) / (NEIGH_L*N_ELEMS) CALL LIST_L_AT_NODE_INFO (L_TO_N_SUM, L_FIRST, L_LAST) IF ( PT_EL_LIST .OR. PT_EL_LIST ) THEN NEIGH_N = MAXVAL (L_TO_N_SUM) NEIGH_P = NEIGH_N ! SINCE NODAL BASED PATCHES !b ALLOCATE ( L_TO_N_NEIGH (NEIGH_N, MAX_NP ) ) CALL FORM_L_ADJACENT_NODES (N_ELEMS, NOD_PER_EL, MAX_NP, & NODES, NEIGH_N, L_TO_N_NEIGH) IF ( PT_EL_LIST ) CALL LIST_ELEMENTS_AT_PT (L_TO_N_NEIGH) !b print *, 'L_TO_N Sparseness Value ', & !b float( count(L_TO_N_NEIGH == 0)) / (NEIGH_N*MAX_NP) DEALLOCATE ( L_TO_N_NEIGH ) !s ? END IF END IF ! BUILD NEIGHBORS, TO DEBUG !--> ** INPUT INITIAL DELETED ELEMENT LIST IF ( GET_DELETED_EL ) CALL INPUT_DELETED_ELEMENTS ! INITIALIZE SYSTEM DOF FOR ITERATIVE SOLUTION !b IF ( N_ITER > 1) CALL START_ITERATION (X, DD_OLD) IF ( N_ITER > 1) THEN IF ( RESTART_STEP > 0 ) THEN THIS_ITER = RESTART_STEP CALL GET_STEP_RESULTS (DD_OLD) ELSE THIS_ITER = 0 CALL START_ITERATION (X, DD_OLD) END IF ELSE DD_OLD = 0.d0 END IF IF ( I_BUG > 0 ) THEN print *,'No iteration start debug active' END IF ! *** BEGIN ITERATION LOOP *** DO ITER = 1, N_ITER ! ======> ======> ======> ======> ======> THIS_ITER = RESTART_STEP + ITER ! Set global value IF ( N_ITER > 1 ) THEN PRINT *, '--------------------------------------------------' PRINT *, ' BEGINNING ITERATION NUMBER ', ITER PRINT *, '--------------------------------------------------' END IF IF ( ITER > 1 ) THEN IF ( N_REACT > 0 ) REWIND (N_REACT) IF ( L_REACT > 0 ) REWIND (L_REACT) IF ( U_FLUX > 0 ) REWIND (U_FLUX ) IF ( N_FILE1 > 0 ) REWIND (N_FILE1) IF ( N_FILE2 > 0 ) REWIND (N_FILE2) IF ( N_FILE3 > 0 ) REWIND (N_FILE3) IF ( N_FILE4 > 0 ) REWIND (N_FILE4) END IF N_L_DEL = COUNT (DELETED) IF ( N_L_DEL > 0 ) THEN print *,' ' print *,'NUMBER OF ELEMENTS OMITTED = ', N_L_DEL print *,'** OMITTED ELEMENT LIST ** ' DO IE = 1, L_S_TOT IF ( DELETED (IE) ) print *, IE END DO END IF IF ( N_L_DEL > OLD_DEL ) THEN ! more elements omitted this iter OLD_DEL = N_L_DEL ! BUILD MESH GEOMETRIC PROPERTIES (INERTIA MATRIX) CALL ASSEMBLE_MESH_GEOM (NODES, X) WRITE (6, 5) N_D_FRE 5 FORMAT (/, 'TOTAL NUMBER OF SYSTEM EQUATIONS ....',I10 ) END IF ! number of active elements reduced ! ** STOP IF DATA CHECK ONLY RUN ** IF ( INPUT_ONLY ) THEN IF ( N_WARN == 0 ) THEN PRINT *, 'END OF MODEL_F90 DATA CHECK (NO WARNINGS)' ELSE PRINT *, 'END MODEL_F90 DATA CHECK, WITH ', N_WARN,' WARNINGS' END IF STOP 'REQUESTED input_only CHECK COMPLETED. EXITING' END IF ! Data check only ! ALLOCATE AND ZERO THE SYSTEM MATRICES ALLOCATE ( M_FULL (N_D_FRE, N_D_FRE)) !b ALLOCATE ( S_FULL (N_D_FRE, N_D_FRE)) !b M_FULL = 0.d0 ; S_FULL = 0.d0 IF ( I_BUG > 0 ) PRINT *, ' ALLOCATED FULL REAL ARRAYS' ! INITIALIZE THE SYSTEM SOURCE MATRIX DD = 0.0_DP CC = 0.0_DP ! zero load !b IF ( IN_RHS > 0 ) CC = CC_IN ! initialize load !b IF ( I_BUG > 0 ) CALL SUM_RHS (CC) !--> *** CALCULATE AND ASSEMBLE ELEMENT MATRICES *** !--> *** GENERATE POST SOLUTION MATRICES & STORE *** CALL ASSEMBLE_FULL_EQS (NODES, S_FULL, CC, X, & !b DD_OLD, ITER, M_FULL) IF ( PT_MASSES ) THEN ! Add diagonal mass terms DO J = 1, N_D_FRE M_FULL (J, J) = M_FULL (J, J) + MASSES (J) END DO END IF ! point masses were input !--> ** APPLY BOUNDARY CONSTRAINTS TO NODAL PARAMETERS ** IF ( MAX_ACT > 1 ) THEN PRINT *,'WARNING: MPC NOT IMPLEMENTED FOR EIGENANALYSIS' N_WARN = N_WARN + 1 END IF ! mpc requested !--> TYPE 1 D(L1) = C1 NTYPE = N_RES_EQ (1) ! number of type_1 BC IF ( NTYPE > 0 ) THEN DO J = 1, NTYPE ! zero rows and columns EQ_N = NDEX_C_EQ (1, J) VALUE = N_D_FRE + J ! eigenvalue flag S_FULL (:, EQ_N) = 0.d0 S_FULL (EQ_N, :) = 0.d0 S_FULL (EQ_N, EQ_N) = VALUE M_FULL (:, EQ_N) = 0.d0 M_FULL (EQ_N, :) = 0.d0 M_FULL (EQ_N, EQ_N) = 1.d0 END DO ! J IF ( I_BUG == 1 .AND. N_D_FRE < M_SHOW ) THEN PRINT *,'DEBUG PRINT AFTER BC APPLICATION' CALL RPRINT (S_FULL, N_D_FRE, N_D_FRE, 1) CALL RPRINT (M_FULL, N_D_FRE, N_D_FRE, 1) END IF ! debug IF ( EIGEN_SHOW > (N_D_FRE - NTYPE) ) THEN EIGEN_SHOW = N_D_FRE - NTYPE PRINT *,'NOTE: EIGEN_SHOW REDUCED BY EBC TO ', EIGEN_SHOW END IF END IF ! Type 1 > 0 !--> *** SOLVE FOR EIGEN ANALYSIS RESULTS *** ! CALL JACOBI_GENERAL (S_FULL, M_FULL, EIGEN_VEC, EIGEN_VAL, & N_D_FRE, NS_MAX, NS_USED) ! FIND THEIR ORDER FROM LOWEST TO HIGHEST CALL Real_Vector_Order (N_D_FRE, EIGEN_VAL, EIGEN_ORDER) ! RE_ORDER FOR ESSENTIAL BC IF ( NTYPE > 0 ) CALL SHIFT_ORDER_FOR_CEQ (N_D_FRE, & EIGEN_ORDER, NTYPE, NDEX_C_EQ (1, 1:NTYPE)) ! Re-order once and for all for output, post-processing EIGEN_VAL = EIGEN_VAL (EIGEN_ORDER) EIGEN_VEC (:, :) = EIGEN_VEC (:, EIGEN_ORDER) ! Scale maximum value absolute value to unity !! EIGEN_LOC = MAXLOC ( ABS (EIGEN_VEC), DIM = 1 ) ! xxx need to use nodal dof 1, not slope, etc ! DO J = 1, N_D_FRE ! EIGEN_VEC (:, J) = EIGEN_VEC (:, J) / EIGEN_VEC (EIGEN_LOC(J), J) ! END DO ! J !b do this scaling in plots, may need physical scale ! !--> *** SOLUTION COMPLETE, GET SYSTEM REACTIONS *** ! IF ( I_BUG > 0 ) THEN print *,'No reactions in eigen solution' END IF ! *** PRINT RESULTS *** PRINT *, ' ' PRINT *, '*** EIGEN ANALYSIS RESULTS ***' PRINT *, 'GENERALIZED JACOBI SWEEPS USED = ', NS_USED PRINT *, 'NOTE: EIGEN SAVES, VALUES IN FILE 000, VECTORS IN 001 ...' IF ( EIGEN_ROOT ) THEN ! Use square root PRINT * , EIGEN_SHOW, & ' EIGENVALUE SQUARE ROOTS, Radians/Sec (LOW TO HIGH):' DO J = 1, EIGEN_SHOW IF ( EIGEN_VAL (J) < 0.d0 ) THEN ! possible user error PRINT *,'WARNING: NEGATIVE EIGENVALUE SET TO ZERO' N_WARN = N_WARN + 1 ; EIGEN_VAL (J) = 0.d0 ELSE EIGEN_VAL (J) = SQRT ( EIGEN_VAL (J) ) END IF END DO ELSE ! Use original values PRINT * , EIGEN_SHOW, & ' EIGENVALUES, Radians/Sec (LOW TO HIGH, OMIT EBC):' END IF ! need square root of eigenvalue WRITE (N_PRT,6) (J, EIGEN_VAL (J), J = 1, EIGEN_SHOW) 6 FORMAT ( 4(1X, I3, 1X, 1PE12.5), & (/, 4(1X, I3, 1X, 1PE12.5))) PRINT * , EIGEN_SHOW, ' VALUES IN Cycles/Sec:' WRITE (N_PRT,6) (J, (EIGEN_VAL (J)/6.283185d0), J=1,EIGEN_SHOW) THIS_STEP = 0 ! Global flag for temp file number CALL SAVE_STEP_RESULTS (EIGEN_VAL) PRINT *, ' ' PRINT *, EIGEN_SHOW, ' EIGENVECTORS (FROM LOWEST TO HIGHEST):' DO J = 1, EIGEN_SHOW !b WRITE (N_PRT, 7) J, (EIGEN_VEC (I, J), I = 1, N_D_FRE) !b 7 FORMAT (I6, 5(1X, 1PE12.5), & !b ( (/, 6X, 5(1X, 1PE12.5)))) WRITE (N_PRT, 7) J, N_G_DOF 7 FORMAT ('MODE SHAPE NUMBER ', I3, /, 'NODE, ', I2, ' DOF') DO K = 1, MAX_NP !b INDEX (1:N_G_DOF) = GET_INDEX_AT_PT (K) ! pt size !b WRITE (N_PRT, 8) EIGEN_VEC (INDEX (1:N_G_DOF), J) LIST = GET_INDEX_AT_PT (K) ! pt size WRITE (N_PRT, 8) K, EIGEN_VEC (LIST, J) 8 FORMAT ( I4, 6(ES12.4) ) END DO ! SAVE EIGEN VECTORS TO node_result_###.tmp THIS_STEP = J ! Global flag for temp file number CALL SAVE_STEP_RESULTS (EIGEN_VEC (:, J)) END DO ! J ! ** BAR CHART OF SOLUTION ** IF ( BAR_CHART ) THEN PRINT *, '' ; PRINT *, 'NODAL EIGEN-SOLUTION GRAPHS' IBAR = 1 ! DEFAULT CONTROL BAR_DOF = N_G_DOF ! Plot all components DO J = 1, EIGEN_SHOW PRINT *, ' ' PRINT *, 'EIGENVECTOR NUMBER ', J, ': ', EIGEN_VAL (J) ! EXTRACT THIS EIGEN VECTOR DD = EIGEN_VEC (:, J) !b DO IPARM = 1, N_G_DOF DO IPARM = 1, BAR_DOF CALL NODAL_BAR_CHART (IBAR, IPARM, X, DD) END DO END DO END IF ! bar charts !b CALL MAX_AND_MIN_RANGE (DD) !b IF ( NPT_WRT == 1) THEN !b IF ( USE_EXACT ) THEN !b CALL WRITE_BY_NODES_WITH_EXACT (X, DD) !b ELSE !b CALL WRITE_BY_NODES (X, DD) !b END IF ! EXACT SOLUTION PROVIDED (IN EXACT_SOLUTION) !b END IF ! LIST BY NODES !b IF ( LEM_WRT == 1) CALL WRITE_BY_ELEMENTS (DD, NODES) !b IF ( U_PLT3 > 0) CALL SAVE_RESULTS (DD) IF ( U_PLT8 > 0 .AND. USE_EXACT ) CALL & SAVE_EXACT_SOLUTION_AT_NODES (X) !b IF ( U_PLT9 > 0 .AND. USE_EXACT ) CALL & !b SAVE_EXACT_FLUXES_AT_NODES (X) !b ! RESET THE REACTION DATA AND SQ MATRIX FOR NEXT ITERATION !b IF ( N_REACT > 0) REWIND (N_REACT) !a IF ( .NOT. SYMMETRIC ) DEALLOCATE ( S_SKY_LOWER ) !a DEALLOCATE ( S_SKY ) !a IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED SKYLINE REAL ARRAYS' !b ! ** OPTIONAL ELEMENT REACTION RECOVERY *** !b IF ( L_REACT > 0 ) CALL LIST_REACTIONS_AT_ELEMS (NODES, DD, CC_IN) !--> *** POST SOLUTION CALCULATIONS *** IF ( U_FLUX > 0 .OR. N_FILE1 > 0 .OR. N_FILE2 > 0 ) THEN !sc IF ( U_SCPR > 0 .AND. N_QP > 0 ) THEN ! POSSIBLE SCP RECOVERY, PREPARE RECORDS SCP_PAD = EST_SCP_PAD_SIZE () SCP_RECL = SET_SCP_RECORD_SIZE () CALL INITIALIZE_SCP_RECORDS ! FOR FILLING ELSE ! REMIND USER PRINT *, ' ' PRINT *, 'NOTE: SUPER_CONVERGENT PATCH RECOVERY NOT ACTIVE' PRINT *, ' (UNIT U_SCPR = 0 OR N_QP = 0)' PRINT *, ' ' END IF ! U_SCPR, POSSIBLE SCP RECOVERY ! POST-PROCESS APPLICATION (AND USE SCP AVERAGING) IF ( N_FILE1 > 0 .OR. N_FILE2 > 0 ) CALL & POST_PROCESS_EIGENVECTORS (NODES, EIGEN_VAL, EIGEN_VEC, & EIGEN_SCALE, DD, DD_OLD, ITER, X) !b! SELECT ONE MODE FOR FLUX AVERAGES AND/OR ERROR ESTIMATE !b PRINT *,' ' !b PRINT *,'MODE SELECTED FOR AVERAGES AND/OR ERROR ESTIMATE=',& !b EIGEN_SCP IF ( EIGEN_SCP < 1 ) THEN PRINT *,'NOTE: NO MODE FLUX AVERAGES, ADD KEYWORD eigen_scp' ELSE PRINT *,'MODES FOR AVERAGES AND/OR ERROR ESTIMATE = ', & EIGEN_SCP DO K = 1, EIGEN_SCP ! ----> ----> ----> ----> ----> ----> PRINT *,' PROCESSING MODE NUMBER ', K DD (:) = EIGEN_VEC (:, K) ! and DD_OLD = 0 IF ( DEBUG_SCP ) THEN PRINT *, 'USING MODE SHAPE ', DD END IF ! POST-PROCESS GRADIENTS & GENERATE SCP DATA FOR AVERAGING IF ( U_FLUX > 0 ) CALL POST_PROCESS_GRADS (NODES, DD, ITER) IF ( SCP_RECORDS_SAVED ) THEN ! *** GET SCP NODAL FLUX AVERAGES *** IF ( SCP_NEIGH_PT ) THEN ! NODE BASED PATCH NEIGHBORS PRINT *, 'NOTE: SCP BASED ON NODE PATCHES' N_PATCH = MAX_NP CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_N_NEIGH, & SCP_AVERAGES) ELSE IF ( SCP_NEIGH_FACE ) THEN ! FACE BASED PATCHES PRINT *, 'NOTE: SCP BASED ON ELEMENT FACING PATCHES' N_PATCH = N_ELEMS ! Only face neighbors in L_TO_L_NEIGH CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH,& SCP_AVERAGES) ELSE ! ELEMENT BASED NEIGHBORS (DEFAULT) PRINT *, 'NOTE: SCP BASED ON ELEMENT PATCHES' N_PATCH = N_ELEMS ! All element neighbors in L_TO_L_NEIGH CALL CALC_SCP_AVE_NODE_FLUXES (NODES, X, L_TO_L_NEIGH,& SCP_AVERAGES) END IF ! SCP_NEIGH_PT, BASIS FOR PATCH ELEMENTS IF ( USE_EXACT_FLUX ) THEN IF ( GRAD_BASE_ERROR ) PRINT *, & 'NOTE: FLUX IS ACTUALLY GRADIENT HERE' CALL LIST_SCP_AVE_AND_EXACT_FLUXES (X, SCP_AVERAGES) CALL SAVE_SCP_NODE_AVERAGES (SCP_AVERAGES) IF ( SCP_2ND_DERIV ) THEN CALL LIST_SCP_AVE_EX_FLUX_GRAD (X, SCP_AVERAGES) CALL SAVE_SCP_AVE_EX_FLUX_GRAD (X, SCP_AVERAGES) END IF ! 2ND DERIVATIVES ELSE IF ( GRAD_BASE_ERROR ) THEN CALL LIST_SCP_AVE_GRADIENTS (X, SCP_AVERAGES) ELSE CALL LIST_SCP_AVE_FLUXES (X, SCP_AVERAGES) ! IF ( SCP_2ND_DERIV ) THEN !b ! CALL LIST_SCP_AVE_FLUX_GRAD (X, SCP_AVERAGES) !b ! CALL SAVE_SCP_AVE_FLUX_GRAD (SCP_AVERAGES) !b ! END IF ! 2ND DERIVATIVES !b END IF ! gradient based only CALL SAVE_SCP_NODE_AVERAGES (SCP_AVERAGES) END IF ! USE_EXACT_FLUX, EXACT FLUXES GIVEN ! *** GET ELEMENT ERROR ESTIMATES VIA SCP GRADIENTS *** IF ( .NOT. NO_ERROR_EST ) THEN ! Use the SCP averages CALL SCP_ERROR_ESTIMATES (NODES, X, SCP_AVERAGES, DD, & EL_ENORM, EL_REFINE, ERR_MAX) CALL SAVE_ELEM_ERROR_ESTIMATES (EL_ENORM) CALL SAVE_PT_AVE_ERROR_EST (NODES, EL_ENORM) IF ( BAR_CHART ) THEN IBAR = 1 ; IPARM = 1 PRINT *, '' ; PRINT *,'ELEMENT ERROR ESTIMATE GRAPH' CALL ELEM_BAR_CHART (IBAR, IPARM, IPARM, EL_ENORM) END IF ! BAR_CHART, Print bar charts too END IF ! NO_ERROR_EST, no error est requested END IF ! SCP_RECORDS_SAVED, SCP REQUESTED END DO ! K FOR SCP ON MODES <---- <---- <---- <---- <---- END IF ! EIGEN_SCP, POST PROCESS THE EIGENVALUES END IF ! U_FLUX, N_FILEn, were data saved !b ! *** UPDATE VALUES FOR NEXT ITERATION (IF ANY) *** !b IF ( ITER == 1) THEN ; R_TEST = 1.d0 ; DD_OLD = DD ; END IF !b IF ( ITER > 1) THEN !b WRITE (6, * ) ' ' !b WRITE (6, * ) 'RESULTS FOR ITERATION NUMBER = ', ITER !b CALL CHANGE_IN_SOLUTION (DD, DD_OLD, OLD_NORM, NEW_NORM, & !b DIFF_NORM, RATIO) !b IF ( (RATIO / R_TEST) < CUT_OFF) THEN !b PRINT *, 'SOLUTION HAS CONVERGED, EXIT ITERATIONS' !b EXIT ! ITERATION LOOP !b END IF !b CALL CORRECT_ITERATION (DD, DD_OLD) !b END IF END DO ! N_ITER <====== <====== <====== <====== <====== <====== ! *** PROBLEM COMPLETED, DELETE SCRATCH FILES *** ! IF ( L_REACT > 0) CLOSE (L_REACT, STATUS = "DELETE" ) IF ( N_REACT > 0) CLOSE (N_REACT, STATUS = "DELETE" ) IF ( U_FLUX > 0) CLOSE (U_FLUX , STATUS = "DELETE" ) IF ( U_INTGR > 0) CLOSE (U_INTGR, STATUS = "DELETE" ) IF ( N_FILE1 > 0) CLOSE (N_FILE1, STATUS = "DELETE" ) IF ( N_FILE2 > 0) CLOSE (N_FILE2, STATUS = "DELETE" ) !b IF ( N_FILE3 > 0) CLOSE (N_FILE3) !bc , STATUS = "DELETE" ) !b IF ( N_FILE4 > 0) CLOSE (N_FILE4) !bc , STATUS = "DELETE" ) !b IF ( N_FILE5 > 0) CLOSE (N_FILE5) !bc , STATUS = "DELETE" ) ! **** DEALLOCATE ARRAYS **** CALL DEALLOCATE_GEOM_PROP !b CALL DEALLOCATE_PROPERTIES IF ( NO_PROPERTIES ) THEN !b DEALLOCATE ( DELETED ) !b ELSE !b CALL DEALLOCATE_PROPERTIES !b END IF !b !a IF ( ALLOCATED (I_DIAG) ) DEALLOCATE ( I_DIAG ) !a IF ( ALLOCATED (I_DOF_HI) ) DEALLOCATE ( I_DOF_HI ) !a IF ( I_BUG > 0 ) PRINT *, 'DEALLOCATED SKYLINE INTEGER ARRAYS' IF ( N_WARN > 0 ) PRINT *,'JACOBI_MODEL GAVE ', N_WARN, ' WARNINGS' END SUBROUTINE JACOBI_MODEL SUBROUTINE SHIFT_ORDER_FOR_CEQ (N, ORDER, N_CE, NDEX_C) IMPLICIT NONE INTEGER, INTENT (IN) :: N, N_CE INTEGER, INTENT (IN) :: NDEX_C (N_CE) ! type 1 dof INTEGER, INTENT (INOUT) :: ORDER (N) INTEGER :: J, K, KOUNT, TEMP, SWAPS_MADE KOUNT = N_CE IF ( KOUNT == 0 ) RETURN DO J = 1, N_CE DO K = 1, N-1 IF ( ORDER (K) == NDEX_C(J) ) THEN TEMP = ORDER (K) ORDER (K:N-1) = ORDER (K+1:N) ORDER (N) = TEMP KOUNT = KOUNT - 1 IF ( KOUNT <= 0 ) EXIT END IF END DO ! K IF ( KOUNT <= 0 ) EXIT END DO END SUBROUTINE SHIFT_ORDER_FOR_CEQ SUBROUTINE POST_PROCESS_EIGENVECTORS (NODES, EIGEN_VAL, EIGEN_VEC, & EIGEN_SCALE, DD, DD_OLD, ITER, X) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for PT (LT_PARM, LT_QP), WT (LT_QP), ! G (LT_GEOM, LT_QP), H (LT_N ), ! DLG (LT_PARM, LT_GEOM, LT_QP), ! DLH (LT_PARM, LT_N , LT_QP), ! S (LT_FREE, LT_FREE), C (LT_FREE), ! ELEM_NODES (LT_N) Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Select_Source Use Interface_Header IMPLICIT NONE INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL), ITER REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) REAL(DP), INTENT(IN) :: EIGEN_VAL (N_D_FRE) REAL(DP), INTENT(INOUT) :: EIGEN_VEC (N_D_FRE, N_D_FRE) REAL(DP), INTENT(INOUT) :: DD (N_D_FRE), DD_OLD (N_D_FRE) REAL(DP), INTENT(OUT) :: EIGEN_SCALE (N_D_FRE) INTEGER :: L_PT_PROP (N_NP_FIX) !b ?? INTEGER :: IE, LT, K !b INTEGER :: EIGEN_POST = 3 ! xxx make global INTEGER, SAVE :: NOTE = 0 !xxx make global !b REAL(DP) :: THIS_SCALE = 1.d0 ! Global <---> POST_PROCESS_ELEM ! Automatic Arrays REAL(DP) :: PRT_MAT (MISC_FL) ! Allocatable Arrays REAL(DP), ALLOCATABLE :: PRT_L_PT (:, :) ! VARIABLES: ! COORD = COORDINATES OF ALL NODES ON ELEMENT ! D = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT ! DD = ARRAY OF SYSTEM DEGREES OF FREEDOM ! DGH = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS ! DLG = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION ! DLH = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS ! E = CONSTITUTIVE MATRIX ! EB = PRODUCT OF E*B ! FLT_MISC = SYSTEM STORAGE OF FLOATING PT MISC PROP ! FLUX_LT = FLUX AT ELEMENT NODES FROM A SCP: (SCP_FIT, LT_N) ! G = GEOMETRIC INTERPOLATION FUNCTIONS ! H = SOLUTION INTERPOLATION FUNCTIONS ! H_INTG = INTEGRAL OF INTERPOLATION FUNCTIONS ! INDEX = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT ! ITER = CURRENT ITERATION NUMBER ! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT ! L_PT_PROP = INTEGER PROPERTIES AT EACH ELEMENT NODE ! MAX_NP = TOTAL NUMBER OF NODES ! MISC_FIX = MISCELLANEOUS INTEGER SYSTEM PROPERTIES ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM ! N_ELEMS = NUMBER OF ELEMENTS IN SYSTEM ! N_EL_FRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT ! N_G_DOF = NUMBERS PARAMETERS PER NODE ! N_GEOM = NUMBER OF GEOMETRY NODES ! N_ITER = MAX NUMBER OF ITERATIONS ! N_MAT = NUMBER OF MATERIAL TYPES ! NODES = ELEMENT INCIDENCES OF ALL ELEMENTS ! N_PARM = DIMENSION OF PARAMETRIC SPACE ! N_QP = NUMBER OF QUADRATURE POINTS ! 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 = UNIT FOR POST SOLUTION MATRICES STORAGE ! N_FILE3-5 = OPTIONAL UNITS FOR USER (USED WHEN > 0) ! PRT_L_PT = REAL PROPERTIES AT ELEMENT NODES ! PRT_MAT = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER ! PT = QUADRATURE COORDINATES ! S = ELEMENT SQUARE MATRIX ! SCP_AVERAGES = AVERAGED FLUX ITEMS AT ALL NODES ! SCP_FIT = NUMBER IF TERMS BEING FIT, = N_R_B USUALLY ! STRAIN = STRAIN OR GRADIENT VECTOR ! STRAIN_0 = INITIAL STRAIN OR GRADIENT VECTOR ! STRESS = STRESS VECTOR ! WT = QUADRATURE WEIGHTS ! XYZ = SPACE COORDINATES AT A POINT IF ( EIGEN_POST < 1 ) RETURN ! no post processing requested ! N_FILE1 OR 2 MUST BE > 0 IF ( N_FILE1 == 0 .AND. N_FILE2 == 0 ) THEN PRINT *, 'MISSING KEYWORD post_el OR post_1 OR post_2' STOP 'NO N_FILE1 OR NFILE2 IN POST_PROCESS_EIGENVECTORS' !e ELSE !e IF ( N_FILE1 > 0 ) REWIND (N_FILE1) !e IF ( N_FILE2 > 0 ) REWIND (N_FILE2) END IF WRITE (N_PRT, "(/,'** BEGIN ELEMENT EIGENVECTORS POST PROCESSING **')") WRITE (N_PRT, "('FOR ', I3, ' EIGENVALUES')") EIGEN_POST LT = 1 ; LAST_LT = 0 ! INITIALIZE ELEMENT TYPES EIGEN_SCALE = 1.d0 ! INITIALIZE SCALING DO K = 1, EIGEN_POST ! loop over selected eigenvalues IF ( N_FILE1 > 0 ) REWIND (N_FILE1) IF ( N_FILE2 > 0 ) REWIND (N_FILE2) IF ( U_FLUX > 0 ) THEN REWIND (U_FLUX ) ELSE PRINT *, 'NOTE: UNIT U_FLUX NOT AVAILABLE FOR EIGENVECTORS' END IF PRINT *, ' ' PRINT *, 'EIGEN_VALUE NUMBER ', K, EIGEN_VAL (K) ! COPY CURRENT EIGEN_VECTOR INTO GLOBAL DOF VECTOR DD DD (:) = EIGEN_VEC (:, K) !--> LOOP OVER ELEMENTS DO IE = 1, L_S_TOT ! for elements and any boundary segments CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG IF ( IE == (N_ELEMS + 1) ) WRITE (N_PRT, & "('** BEGIN EIGENVECTORS BOUNDARY SEGMENTS POST PROCESSING **')") !--> 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 ! GET CONTROLS FOR THIS TYPE CALL SET_ELEM_TYPE_INFO (LT) !--> Allocate Arrays for this type IF ( N_NP_FLO > 0 ) THEN !e IF ( ALLOCATED (PRT_L_PT) ) DEALLOCATE ( PRT_L_PT ) !ba ALLOCATE ( PRT_L_PT ( LT_N, N_NP_FLO ) ) !b need the same for integer prop !b IF ( I_BUG > 0 ) PRINT *, ALLOCATED ( PRT_L_PT ) END IF ! n_np_flo !e !--> GET QUADRATURE RULE FOR ELEMENT TYPE AND SHAPE IF ( LT_QP > 0 ) CALL GET_ELEM_QUADRATURES END IF ! a new element type ! INITIALIZE AND CHECK IF ACTIVE IF ( DELETED (IE) ) CYCLE ! to skip deleted element !--> EXTRACT ELEMENT NODE NUMBERS ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES) !9 !--> EXTRACT NODAL COORDINATES CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES) !--> CALCULATE DEGREE OF FREEDOM NUMBERS INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES) !--> EXTRACT NODAL PARAMETERS OF THE ELEMENT D = GET_ELEM_DOF (DD) !--> EXTRACT NODAL POINT PROPERTIES (IF ANY) IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, & PRT_L_PT, L_PT_PROP, ELEM_NODES) !--> EXTRACT MATERIAL PROPERTIES (IF ANY) IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT) ! !--> PERFORM PROBLEM DEPENDENT CALCULATIONS AND OUTPUT ! IF ( IE <= N_ELEMS ) THEN !b IF ( DEBUG_POST_EL ) PRINT *,'Debug IS_STD ', IS_STD !b IF ( POST_EL .AND. IS_ELEMENT ) CALL & IF ( (EIGEN_POST > 0) .AND. IS_ELEMENT ) CALL & SELECT_POST_PROCESS_ELEM (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE, & DD, DD_OLD) !b !b L_PT_PROP, ITER, IE) !b ELSE IF ( POST_MIXED .AND. IS_MIXED_BC ) THEN !b IF ( DEBUG_MIX_SQ ) PRINT *,'Debug IS_MIX ', IS_MIX CALL SELECT_POST_PROCESS_MIXED (PRT_L_PT, PRT_MAT, & !b L_PT_PROP, ITER, IS_MIX) !b ELSE !b IF ( NOTE == 0 ) THEN ; NOTE = 1; PRINT *,'WARNING, SELECT_POST_PROCESS UNKNOWN OPTION' !b PRINT *, THIS_EL, 'NOT A STANDARD OR MIXED SEGMENT ELEMENT' PRINT *,'IS_STD, IS_MIX, IS_SEG', IS_STD, IS_MIX, IS_SEG PRINT *,'POST_EL, IS_ELEMENT', POST_EL, IS_ELEMENT PRINT *,'POST_MIXED, IS_MIXED_BC', POST_MIXED, IS_MIXED_BC N_WARN = N_WARN + 1 !b END IF ! Warn only once END IF ! Element class !b THIS_SCALE = 1.d0 ! xxx have this set in SELECT_POST_PROCESS_* IF ( IE == N_ELEMS ) WRITE (N_PRT, & "('ELEMENT POST PROCESSING COMPLETE.')") IF ( IE == L_S_TOT .AND. (N_MIXED > 0 .OR. N_F_SEG > 0) ) & WRITE (N_PRT, "('PROCESSING COMPLETE FOR EIGENVALUE ', K, /)") END DO ! over all elements and boundary segments ! RE-SCALE THE EIGEN VECTORS EIGEN_SCALE (K) = THIS_SCALE !b EIGEN_VEC (:, K) = THIS_SCALE * EIGEN_VEC (:, K) EIGEN_VEC (:, K) = THIS_SCALE * DD (:) END DO ! for K eigenvalues !--> POST-PROCESSING COMPLETE IF ( ALLOCATED (PRT_L_PT) ) DEALLOCATE ( PRT_L_PT ) END SUBROUTINE POST_PROCESS_EIGENVECTORS