PROGRAM RP_ADAPTIVE_FE_SOLVER ! Copyright 1999 M. Singh, J. E. Akin. All rights reserved. ! Rice University, MEMS Department, Houston, TX USA ! Shareware. Not for commerical use. akin@rice.edu USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE SPR_TYPES_MODULE USE SYSTEM_EQNS_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE) :: NODE_LST TYPE (EDGE_LST_PTR_TYPE) :: EDGE_LST TYPE (ELEMENT_LST_PTR_TYPE) :: ELEMENT_LST TYPE (SPR_DATA_LST_PTR_TYPE) :: SPR_DATA_LST TYPE (NODE_LST_PTR_TYPE) :: INITIALIZE_NODE_LIST TYPE (EDGE_LST_PTR_TYPE) :: INITIALIZE_EDGE_LIST TYPE (ELEMENT_LST_PTR_TYPE) :: INITIALIZE_ELEMENT_LIST TYPE (SPR_DATA_LST_PTR_TYPE) :: INITIALIZE_SPR_DATA_LIST INTEGER, PARAMETER :: NO_SPTD_TYPES = 1 INTEGER, DIMENSION (NO_SPTD_TYPES) :: CNSTRNT_CNT REAL (DP), DIMENSION (:), ALLOCATABLE :: TYPE_1_VALS INTEGER, DIMENSION (:), ALLOCATABLE :: TYPE_1_DOFS INTEGER :: NO_SYS_DOFS INTEGER :: RESET_SYS_DOF_INFO REAL (DP) :: ACOND REAL (DP) :: DETCF INTEGER :: IDETEX INTEGER :: NEGEIG INTEGER :: IFAIL INTEGER :: CRRNT_ITER LOGICAL :: DONE REAL (DP) :: CUTOFF REAL (DP), DIMENSION (:), ALLOCATABLE :: EL_ENORM ! Initialize the necessary data structures NODE_LST = INITIALIZE_NODE_LIST () EDGE_LST = INITIALIZE_EDGE_LIST () ELEMENT_LST = INITIALIZE_ELEMENT_LIST () SPR_DATA_LST = INITIALIZE_SPR_DATA_LIST () ! Read and echo the problem control data CALL READ_CONTROL_DATA () ALLOCATE (EL_ENORM(NO_ELEMENTS)) ! Input the geometry nodes and element topology CALL INPUT_NODES_AND_ELEMENTS (NODE_LST, EDGE_LST, ELEMENT_LST) ! Decrypt each nodal constraint indicator and get number of constraints CALL COUNT_CONSTRAINTS (NODE_LST, NO_SPTD_TYPES, CNSTRNT_CNT) ! Read nodal parameter constraint equations (only type 1 for now) ALLOCATE (TYPE_1_VALS(CNSTRNT_CNT(1))) ALLOCATE (TYPE_1_DOFS(CNSTRNT_CNT(1))) CALL INPUT_TYPE_1_CNSTRNT_DATA (CNSTRNT_CNT, TYPE_1_VALS, TYPE_1_DOFS) ! Read and echo the property data (if furnished) CALL INPUT_PROPERTY_DATA () CRRNT_ITER = 0 DONE = .FALSE. DO IF (DONE) EXIT CRRNT_ITER = CRRNT_ITER + 1 WRITE (6,"(A,I2)") & " I T E R A T I O N # ", CRRNT_ITER WRITE (6,"(A)") & " -----------------------" ! Write out an SDRC IDEAS universal file (dataset #15) ! CALL WRITE_IDEAS_NODE_UNV (NODE_LST,CRRNT_ITER) ! Open any necessary files IF (POST_PROC_IO_UNIT1 > 0) OPEN (POST_PROC_IO_UNIT1,FORM="UNFORMATTED") IF (POST_PROC_IO_UNIT2 > 0) OPEN (POST_PROC_IO_UNIT2,FORM="UNFORMATTED") ! Reset the system degree-of-freedom information NO_SYS_DOFS = RESET_SYS_DOF_INFO (NODE_LST) ! Allocate the system equation data structures CALL ALLOC_SYSTEM_EQNS_DATA (NO_SYS_DOFS, ELEMENT_LST) ! Input the initial forcing or right-hand side vector IF (CRRNT_ITER == 1 .AND. RHS_FLAG > 0) CALL INPUT_RHS (NO_SYS_DOFS, SYS_SRC_VCTR) ! Dump the data structures (for debugging purposes) IF (DEBUG_MODE > 0) THEN CALL PRINT_NODE_LIST (NODE_LST) CALL PRINT_EDGE_LIST (EDGE_LST) CALL PRINT_ELEMENT_LIST (ELEMENT_LST) END IF ! Assemble the system matrices CALL ASSEMBLE_SYSTEM_EQNS (NO_SYS_DOFS,ELEMENT_LST,SPR_DATA_LST) ! Apply the boundary conditions IF (CNSTRNT_CNT(1) > 0) THEN CALL APPLY_TYPE_1_BCS (NO_SYS_DOFS, DIAG_LCTNS, CNSTRNT_CNT(1), & TYPE_1_VALS, TYPE_1_DOFS, SYS_SQR_MTRX, SYS_SRC_VCTR) END IF ! Perform the system solution CALL SKYFAC (NO_SYS_DOFS, DIAG_LCTNS, .TRUE., .TRUE., 2, SYS_SQR_MTRX, & ACOND, DETCF, IDETEX, NEGEIG, IFAIL) IF (DEBUG_MODE > 0) THEN PRINT*,"*** DEBUG INFORMATION" PRINT*,"acond = ",acond PRINT*,"detcf = ",detcf PRINT*,"idetex = ",idetex PRINT*,"negeig = ",idetex PRINT*,"ifail = ",ifail PRINT*,"*********************" END IF CALL SKYSOL (NO_SYS_DOFS, DIAG_LCTNS, SYS_SQR_MTRX, SYS_SRC_VCTR, & SYS_SOL_VCTR) ! Print results IF (NODAL_PRINT_FLAG == 0) CALL OUTPUT_BY_NODES (NODE_LST, NO_SYS_DOFS, & SYS_SOL_VCTR) IF (ELEMENTAL_PRINT_FLAG == 0) CALL OUTPUT_BY_ELEMENTS (ELEMENT_LST, & NO_SYS_DOFS, SYS_SOL_VCTR) ! Perform any post-solution calculations CALL PERFORM_POST_PRCSSNG (ELEMENT_LST,NO_SYS_DOFS,SYS_SOL_VCTR,SPR_DATA_LST) CALL ESTIMATE_ERROR (NO_SYS_DOFS,ELEMENT_LST,NODE_LST,SPR_DATA_LST, & EL_ENORM,CUTOFF) CALL ADAPT_MESH (ELEMENT_LST,EDGE_LST,NODE_LST,EL_ENORM,CUTOFF) CALL DEALLOC_SYSTEM_EQNS_DATA () IF (POST_PROC_IO_UNIT1 > 0) CLOSE (POST_PROC_IO_UNIT1,STATUS="DELETE") IF (POST_PROC_IO_UNIT2 > 0) CLOSE (POST_PROC_IO_UNIT2,STATUS="DELETE") IF (CRRNT_ITER == MAX_ITERS) DONE = .TRUE. END DO ! Clean up DEALLOCATE (EL_ENORM) DEALLOCATE (TYPE_1_VALS) DEALLOCATE (TYPE_1_DOFS) WRITE (6,"(//,A)") & " E N D O F R U N" WRITE (6,"(A)") & " -----------------" END PROGRAM RP_ADAPTIVE_FE_SOLVER SUBROUTINE ADAPT_MESH (EL_LST,EDG_LST,NODE_LST,EL_ENORM,CUTOFF) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDG_LST TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST REAL (DP), DIMENSION (NO_ELEMENTS), INTENT (IN) :: EL_ENORM REAL (DP), INTENT (IN) :: CUTOFF IF (ADAPT_FLAG > 0) THEN CALL P_ADAPT (EL_LST,EDG_LST,NODE_LST,EL_ENORM,CUTOFF) END IF END SUBROUTINE ADAPT_MESH SUBROUTINE ALLOC_ELMNT_ASSMBLY_DATA (ELEMENT_PTR,SPR_DATA_LST) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE ELMNT_ASSMBLY_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER :: NO_NODES_E1 INTEGER :: NO_NODES_E2 INTEGER :: NO_NODES_E3 INTEGER :: NO_NODES_E4 INTEGER :: GET_NO_NODES_ON_EDGE INTEGER :: MAX_NDS_R INTEGER :: MAX_NDS_S ALLOCATE (DOF_MAP(NO_EL_DOFS)) ALLOCATE (NODE_CRDS(NO_AN_NODES,SOL_SPACE_DIM)) IF (NO_REAL_PRPS_PER_NODE > 0) & ALLOCATE (NODE_PRPS_R(NO_AN_NODES,NO_REAL_PRPS_PER_NODE)) IF (NO_INT_PRPS_PER_NODE > 0) & ALLOCATE (NODE_PRPS_I(NO_AN_NODES,NO_INT_PRPS_PER_NODE)) ALLOCATE (AJ(SOL_SPACE_DIM,SOL_SPACE_DIM)) ALLOCATE (AJINV(SOL_SPACE_DIM,SOL_SPACE_DIM)) ALLOCATE (B_MTRX(NO_ROWS_B,NO_EL_DOFS)) ALLOCATE (BODY(SOL_SPACE_DIM)) ALLOCATE (DGH(SOL_SPACE_DIM,NO_AN_NODES)) ALLOCATE (DLG(2,NO_SYS_GEOM_NODES)) ALLOCATE (DLH(SOL_SPACE_DIM,NO_AN_NODES)) ALLOCATE (E_MTRX(NO_ROWS_B,NO_ROWS_B)) ALLOCATE (EB_MTRX(NO_ROWS_B,NO_EL_DOFS)) ALLOCATE (EL_SQR_MTRX(NO_EL_DOFS,NO_EL_DOFS)) IF (EL_COL_FLAG == 0) ALLOCATE (EL_SRC_VCTR(NO_EL_DOFS)) ALLOCATE (G(NO_SYS_GEOM_NODES)) ALLOCATE (H(NO_AN_NODES)) ALLOCATE (HINTG(NO_AN_NODES)) ALLOCATE (STRAIN(NO_ROWS_B+2)) ALLOCATE (STRAN0(NO_ROWS_B)) ALLOCATE (STRESS(NO_ROWS_B+2)) NO_NODES_E1 = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_1) NO_NODES_E3 = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_3) MAX_NDS_R = MAX (NO_NODES_E1,NO_NODES_E3) ! Worst case scenario ! NO_QP_R = 2*MAX_NDS_R - 2 NO_QP_R = MAX_NDS_R ALLOCATE (PTR(NO_QP_R)) ALLOCATE (WTR(NO_QP_R)) NO_NODES_E2 = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_2) NO_NODES_E4 = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_4) MAX_NDS_S = MAX (NO_NODES_E2,NO_NODES_E4) ! Worst case scenario ! NO_QP_S = 2*MAX_NDS_S - 2 NO_QP_S = MAX_NDS_S ALLOCATE (PTS(NO_QP_S)) ALLOCATE (WTS(NO_QP_S)) IF (ADAPT_FLAG > 0) THEN CALL INSERT_SPR_EL_DATA (SPR_DATA_LST,ELEMENT_PTR % ID,MAX_NDS_R, & MAX_NDS_S) END IF END SUBROUTINE ALLOC_ELMNT_ASSMBLY_DATA SUBROUTINE ALLOC_ERROR_EST_DATA (EL_PTR) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE ERROR_EST_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: EL_PTR INTEGER :: NO_NODES_E1 INTEGER :: NO_NODES_E2 INTEGER :: NO_NODES_E3 INTEGER :: NO_NODES_E4 INTEGER :: MAX_NODES_ON_SIDE INTEGER :: GET_NO_NODES_ON_EDGE ALLOCATE (DOF_MAP(NO_EL_DOFS)) ALLOCATE (NODE_CRDS(NO_EL_AN_NODES,SOL_SPACE_DIM)) NO_NODES_E1 = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_1) NO_NODES_E3 = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_3) NO_NODES_E2 = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_2) NO_NODES_E4 = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_4) MAX_NODES_ON_SIDE = MAX (NO_NODES_E1,NO_NODES_E2,NO_NODES_E3,NO_NODES_E4) NO_QP_R = CEILING (((MAX_NODES_ON_SIDE-1.0)*3.0+1.0)/2.0) NO_QP_S = NO_QP_R ALLOCATE (PTR(NO_QP_R)) ALLOCATE (WTR(NO_QP_R)) ALLOCATE (PTS(NO_QP_S)) ALLOCATE (WTS(NO_QP_S)) ALLOCATE (AJ(SOL_SPACE_DIM,SOL_SPACE_DIM)) ALLOCATE (AJINV(SOL_SPACE_DIM,SOL_SPACE_DIM)) ALLOCATE (DGH(SOL_SPACE_DIM,NO_EL_AN_NODES)) ALLOCATE (DLH(SOL_SPACE_DIM,NO_EL_AN_NODES)) ALLOCATE (EL_ERROR(NO_ROWS_B)) ALLOCATE (EL_SOL_VCTR(NO_EL_DOFS)) ALLOCATE (EL_SPR_NDL_DRVS(NO_EL_AN_NODES,NO_ROWS_B)) ALLOCATE (H(NO_EL_AN_NODES)) ALLOCATE (SGMA_STR(NO_ROWS_B)) ALLOCATE (SGMA_HAT(NO_ROWS_B)) END SUBROUTINE ALLOC_ERROR_EST_DATA SUBROUTINE ALLOC_SYSTEM_EQNS_DATA (NO_SYS_DOFS, ELEMENT_LST) USE DATA_STRUCTURES_MODULE USE SYSTEM_EQNS_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST ALLOCATE (SYS_COL_HGHTS(NO_SYS_DOFS)) ALLOCATE (DIAG_LCTNS(NO_SYS_DOFS+1)) CALL PERFORM_SKYLINE_ANALYSIS (NO_SYS_DOFS, ELEMENT_LST, & SYS_COL_HGHTS, DIAG_LCTNS) ALLOCATE (SYS_SQR_MTRX(DIAG_LCTNS(NO_SYS_DOFS+1))) ALLOCATE (SYS_SRC_VCTR(NO_SYS_DOFS)) ALLOCATE (SYS_SOL_VCTR(NO_SYS_DOFS)) CALL ZERO_OUT_VECTOR (NO_SYS_DOFS, SYS_SRC_VCTR) END SUBROUTINE ALLOC_SYSTEM_EQNS_DATA SUBROUTINE APPLY_TYPE_1_BCS (NO_SYS_DOFS, DIAG_LCTNS, NO_CNSTRNTS, & TYPE_1_VALS, TYPE_1_DOFS, SYS_SQR_MTRX, SYS_SRC_VCTR) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS INTEGER, DIMENSION (NO_SYS_DOFS+1), INTENT (IN) :: DIAG_LCTNS INTEGER, INTENT (IN) :: NO_CNSTRNTS REAL (DP), DIMENSION (NO_CNSTRNTS), INTENT (IN) :: TYPE_1_VALS INTEGER, DIMENSION (NO_CNSTRNTS), INTENT (IN) :: TYPE_1_DOFS REAL (DP), DIMENSION (DIAG_LCTNS(NO_SYS_DOFS+1)), INTENT (INOUT) :: SYS_SQR_MTRX REAL (DP), DIMENSION (NO_SYS_DOFS), INTENT (INOUT) :: SYS_SRC_VCTR INTEGER :: DOF REAL (DP) :: VAL INTEGER :: CRRNT_CNSTRNT INTEGER :: CRRNT_DOF INTEGER :: LRGST INTEGER :: INV INTEGER :: ITOP DO CRRNT_CNSTRNT = 1,NO_CNSTRNTS DOF = TYPE_1_DOFS(CRRNT_CNSTRNT) VAL = TYPE_1_VALS(CRRNT_CNSTRNT) DO CRRNT_DOF = 1,NO_SYS_DOFS LRGST = MAX (CRRNT_DOF,DOF) INV = DIAG_LCTNS(LRGST+1) - ABS(CRRNT_DOF-DOF) ITOP = 1 IF (LRGST > 1) ITOP = DIAG_LCTNS(LRGST) + 1 ! is it outside the skyline IF (INV >= ITOP) THEN SYS_SRC_VCTR(CRRNT_DOF) = SYS_SRC_VCTR(CRRNT_DOF) - & VAL * SYS_SQR_MTRX(INV) SYS_SQR_MTRX(INV) = 0.d0 END IF END DO SYS_SQR_MTRX(DIAG_LCTNS(DOF+1)) = 1.d0 SYS_SRC_VCTR(DOF) = VAL END DO END SUBROUTINE APPLY_TYPE_1_BCS SUBROUTINE ASSEMBLE_SYSTEM_EQNS (NO_SYS_DOFS,ELEMENT_LST,SPR_DATA_LST) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE ELMNT_ASSMBLY_MODULE USE SPR_TYPES_MODULE USE SYSTEM_EQNS_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ELEMENT_PTR INTEGER :: GET_NO_ANALYSIS_NODES IF (NO_REAL_PRPS_PER_EL > 0) ALLOCATE (EL_PRPS_R(NO_REAL_PRPS_PER_EL)) IF (NO_INT_PRPS_PER_EL > 0) ALLOCATE (EL_PRPS_I(NO_INT_PRPS_PER_EL)) CALL ZERO_OUT_VECTOR (DIAG_LCTNS(NO_SYS_DOFS+1),SYS_SQR_MTRX) ELEMENT_PTR => ELEMENT_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (ELEMENT_PTR)) EXIT NO_AN_NODES = GET_NO_ANALYSIS_NODES (ELEMENT_PTR) NO_EL_DOFS = NO_AN_NODES * NO_PARAMS_PER_NODE CALL ALLOC_ELMNT_ASSMBLY_DATA (ELEMENT_PTR,SPR_DATA_LST) CALL GET_SYS_DOFS_FOR_ELEMENT (ELEMENT_PTR,NO_EL_DOFS,DOF_MAP) CALL GET_NODAL_COORDINATES (ELEMENT_PTR,NO_AN_NODES,NODE_CRDS) IF (NO_REAL_PRPS_PER_NODE > 0) & CALL EXTRACT_REAL_NODAL_PT_PRPS (ELEMENT_PTR, NO_AN_NODES, NODE_PRPS_R) IF (NO_INT_PRPS_PER_NODE > 0) & CALL EXTRACT_INT_NODAL_PT_PRPS (ELEMENT_PTR, NO_AN_NODES, NODE_PRPS_I) IF (NO_REAL_PRPS_PER_EL > 0) & CALL EXTRACT_REAL_ELEMENT_PRPS (ELEMENT_PTR, EL_PRPS_R) IF (NO_INT_PRPS_PER_EL > 0) & CALL EXTRACT_INT_ELEMENT_PRPS (ELEMENT_PTR, EL_PRPS_I) CALL GENERATE_EL_SQR_MTRX (ELEMENT_PTR) IF (EL_COL_FLAG == 0) CALL GENERATE_EL_SRC_VCTR () CALL STORE_FOR_POST_PRCSSNG () CALL UPDATE_SYS_SQR_MTRX (NO_EL_DOFS,NO_SYS_DOFS,DOF_MAP,DIAG_LCTNS, & EL_SQR_MTRX, SYS_SQR_MTRX) IF (EL_COL_FLAG == 0) CALL UPDATE_SYS_SRC_VCTR (NO_EL_DOFS, & NO_SYS_DOFS, DOF_MAP, EL_SRC_VCTR, SYS_SRC_VCTR) CALL DEALLOC_ELMNT_ASSMBLY_DATA () ELEMENT_PTR => ELEMENT_PTR % NEXT END DO IF (ALLOCATED (EL_PRPS_R)) DEALLOCATE (EL_PRPS_R) IF (ALLOCATED (EL_PRPS_I)) DEALLOCATE (EL_PRPS_I) END SUBROUTINE ASSEMBLE_SYSTEM_EQNS SUBROUTINE COMPUTE_ERROR (EL_LST,SPR_DATA_LST) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR INTEGER, PARAMETER :: MAX_PTCH_ELS = 9 INTEGER, DIMENSION (MAX_PTCH_ELS) :: PTCH_IDS INTEGER :: M INTEGER :: N REAL (DP), DIMENSION (:,:), ALLOCATABLE :: A REAL (DP), DIMENSION (:), ALLOCATABLE :: Y EL_PTR => EL_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_PTR)) EXIT CALL GET_ELS_IN_PATCH (EL_PTR,PTCH_DFN_FLAG,MAX_PTCH_ELS,PTCH_IDS) CALL DETERMINE_M_AND_N (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS,M,N) print*,"m = ",m print*,"n = ",n ALLOCATE (A(M,N)) ALLOCATE (Y(M)) DEALLOCATE (A) DEALLOCATE (Y) EL_PTR => EL_PTR % NEXT END DO END SUBROUTINE COMPUTE_ERROR SUBROUTINE COMPUTE_GLOBAL_DRVS (NO_AN_NODES, AJINV, DELTA, GLOBAL) USE DBL_PRECISION_MODULE USE CONTROL_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_AN_NODES REAL (DP), DIMENSION (SOL_SPACE_DIM,SOL_SPACE_DIM), INTENT (IN) :: AJINV REAL (DP), DIMENSION (SOL_SPACE_DIM,NO_AN_NODES), INTENT (IN) :: DELTA REAL (DP), DIMENSION (SOL_SPACE_DIM,NO_AN_NODES), INTENT (OUT) :: GLOBAL INTEGER :: I INTEGER :: J INTEGER :: K REAL (DP) :: SUM DO I = 1,SOL_SPACE_DIM DO J = 1,NO_AN_NODES SUM = 0.d0 DO K = 1,SOL_SPACE_DIM SUM = SUM + AJINV(I,K)*DELTA(K,J) END DO GLOBAL(I,J) = SUM END DO END DO END SUBROUTINE COMPUTE_GLOBAL_DRVS SUBROUTINE COMPUTE_JACOBIAN (NO_AN_NODES, DELTA, NODE_CRDS, AJ) USE DBL_PRECISION_MODULE USE CONTROL_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_AN_NODES REAL (DP), DIMENSION (SOL_SPACE_DIM,NO_AN_NODES), INTENT (IN) :: DELTA REAL (DP), DIMENSION (NO_AN_NODES,SOL_SPACE_DIM), INTENT (IN) :: NODE_CRDS REAL (DP), DIMENSION (SOL_SPACE_DIM,SOL_SPACE_DIM), INTENT (OUT) :: AJ INTEGER :: ROW INTEGER :: COL INTEGER :: NDE REAL (DP) :: SUM DO ROW = 1,SOL_SPACE_DIM DO COL = 1,SOL_SPACE_DIM SUM = 0.d0 DO NDE = 1,NO_AN_NODES SUM = SUM + DELTA(ROW,NDE)*NODE_CRDS(NDE,COL) END DO AJ(ROW,COL) = SUM END DO END DO END SUBROUTINE COMPUTE_JACOBIAN SUBROUTINE COMPUTE_JCBN_INVRS_DTRMNT (AJ, AJINV, DTRMNT) USE DBL_PRECISION_MODULE USE CONTROL_MODULE IMPLICIT NONE REAL (DP), DIMENSION (SOL_SPACE_DIM,SOL_SPACE_DIM), INTENT (IN) :: AJ REAL (DP), DIMENSION (SOL_SPACE_DIM,SOL_SPACE_DIM), INTENT (OUT) :: AJINV REAL (DP), INTENT (OUT) :: DTRMNT INTEGER :: ROW INTEGER :: COL IF (SOL_SPACE_DIM == 1) THEN DTRMNT = AJ(1,1) AJINV(1,1) = 1.d0/DTRMNT ELSE IF (SOL_SPACE_DIM == 2) THEN DTRMNT = AJ(1,1)*AJ(2,2) - AJ(1,2)*AJ(2,1) AJINV(1,1) = AJ(2,2)/DTRMNT AJINV(1,2) = -AJ(1,2)/DTRMNT AJINV(2,1) = -AJ(2,1)/DTRMNT AJINV(2,2) = AJ(1,1)/DTRMNT ELSE IF (SOL_SPACE_DIM == 3) THEN AJINV(1,1) = AJ(2,2)*AJ(3,3) - AJ(3,2)*AJ(2,3) AJINV(2,1) = -AJ(2,1)*AJ(3,3) + AJ(3,1)*AJ(2,3) AJINV(3,1) = AJ(2,1)*AJ(3,2) - AJ(3,1)*AJ(2,2) AJINV(1,2) = -AJ(1,2)*AJ(3,3) + AJ(3,2)*AJ(1,3) AJINV(2,2) = AJ(1,1)*AJ(3,3) - AJ(3,1)*AJ(1,3) AJINV(3,2) = -AJ(1,1)*AJ(3,2) + AJ(3,1)*AJ(1,2) AJINV(1,3) = AJ(1,2)*AJ(2,3) - AJ(2,2)*AJ(1,3) AJINV(2,3) = -AJ(1,1)*AJ(2,3) + AJ(2,1)*AJ(1,3) AJINV(3,3) = AJ(1,1)*AJ(2,2) - AJ(2,1)*AJ(1,2) DTRMNT = AJ(1,1)*AJINV(1,1) + AJ(1,2)*AJINV(2,1) + AJ(1,3)*AJINV(3,1) DO ROW = 1,3 DO COL = 1,3 AJINV(ROW,COL) = AJINV(ROW,COL)/DTRMNT END DO END DO ELSE WRITE (6,"(/A/)") "*** Error in COMPUTE_JCBN_INVRS_DTRMNT ***" DTRMNT = 0.d0 END IF END SUBROUTINE COMPUTE_JCBN_INVRS_DTRMNT SUBROUTINE COMPUTE_SRNDPTY_FNS_DRVS (ELEMENT_PTR, NO_AN_NODES, RPT, SPT, H, DLH) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, INTENT (IN) :: NO_AN_NODES REAL (DP), INTENT (IN) :: RPT REAL (DP), INTENT (IN) :: SPT REAL (DP), DIMENSION (NO_AN_NODES), INTENT (OUT) :: H REAL (DP), DIMENSION (SOL_SPACE_DIM,NO_AN_NODES), INTENT (OUT) :: DLH INTEGER :: NODE INTEGER :: DOF_INFO TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR LOGICAL :: RIGHT INTEGER :: NURFN INTEGER :: NSRFN REAL (DP) :: SFVAL REAL (DP), DIMENSION (3) :: DEVAL = 0.d0 INTEGER, DIMENSION (12) :: NNOSI = 0 REAL (DP), DIMENSION (3) :: CDRPT = 0.d0 INTEGER :: GET_NO_NODES_ON_EDGE CDRPT(1) = RPT CDRPT(2) = SPT NNOSI(1) = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_1) NNOSI(2) = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_2) NNOSI(3) = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_3) NNOSI(4) = GET_NO_NODES_ON_EDGE (ELEMENT_PTR % EDG_4) NODE = 1 HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF NURFN = 0 NSRFN = 1 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 DO IF (RIGHT) THEN IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT ELSE IF (ASSOCIATED (TEMP_PTR % LPTR, HEAD_PTR)) EXIT END IF DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NURFN = NURFN + 1 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO NURFN = 0 NSRFN = 2 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF NSRFN = 2 DO IF (RIGHT) THEN IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT ELSE IF (ASSOCIATED (TEMP_PTR % LPTR, HEAD_PTR)) EXIT END IF DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NURFN = NURFN + 1 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO NURFN = 0 NSRFN = 3 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF NSRFN = 3 DO IF (RIGHT) THEN IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT ELSE IF (ASSOCIATED (TEMP_PTR % LPTR, HEAD_PTR)) EXIT END IF DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NURFN = NURFN + 1 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO NURFN = 0 NSRFN = 4 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF NSRFN = 4 DO IF (RIGHT) THEN IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT ELSE IF (ASSOCIATED (TEMP_PTR % LPTR, HEAD_PTR)) EXIT END IF DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NURFN = NURFN + 1 CALL DERSHAFN (NNOSI, NURFN, NSRFN, 2, CDRPT, SFVAL, DEVAL) H(NODE) = SFVAL DLH(1,NODE) = DEVAL(1) DLH(2,NODE) = DEVAL(2) NODE = NODE + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO END SUBROUTINE COMPUTE_SRNDPTY_FNS_DRVS SUBROUTINE CONSTRUCT_COMPLETE_AT_NODE (N,DGR,X,Y,P) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: DGR REAL (DP), INTENT (IN) :: X REAL (DP), INTENT (IN) :: Y REAL (DP), DIMENSION (N), INTENT (OUT) :: P INTEGER :: I INTEGER :: J INTEGER :: INDX REAL (DP) :: XPRT REAL (DP) :: YPRT INDX = 0 DO I = 0,DGR DO J = 0,I INDX = INDX + 1 IF ((I-J) == 0) THEN XPRT = 1.d0 ELSE IF (X == 0.d0) THEN XPRT = 0.d0 ELSE XPRT = X**(I-J) END IF IF ((J) == 0) THEN YPRT = 1.d0 ELSE IF (Y == 0.d0) THEN YPRT = 0.d0 ELSE YPRT = Y**J END IF P(INDX) = XPRT * YPRT END DO END DO END SUBROUTINE CONSTRUCT_COMPLETE_AT_NODE SUBROUTINE CONSTRUCT_SHPFNCTN_AT_NODE (N,MAX_NODES_R,MAX_NODES_S,X,Y,P) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: MAX_NODES_R INTEGER, INTENT (IN) :: MAX_NODES_S REAL (DP), INTENT (IN) :: X REAL (DP), INTENT (IN) :: Y REAL (DP), DIMENSION (N), INTENT (OUT) :: P INTEGER :: I INTEGER :: INDX INDX = 0 DO I = 0,MAX_NODES_R-1 INDX = INDX + 1 IF (I == 0) THEN P(INDX) = 1.d0 ELSE IF (X == 0.d0) THEN P(INDX) = 0.d0 ELSE P(INDX) = X**I END IF END DO DO I = 0,MAX_NODES_R-1 INDX = INDX + 1 IF (I == 0) THEN P(INDX) = Y ELSE IF (X == 0.0 .OR. Y == 0) THEN P(INDX) = 0.d0 ELSE P(INDX) = X**I * Y END IF END DO DO I = 3,MAX_NODES_S INDX = INDX + 1 IF (Y == 0.d0) THEN P(INDX) = 0.d0 ELSE P(INDX) = Y**(I-1) END IF END DO DO I = 3,MAX_NODES_S INDX = INDX + 1 IF (Y == 0.0 .OR. X == 0) THEN P(INDX) = 0.d0 ELSE P(INDX) = X * Y**(I-1) END IF END DO END SUBROUTINE CONSTRUCT_SHPFNCTN_AT_NODE MODULE CONTROL_MODULE USE DBL_PRECISION_MODULE CHARACTER (LEN = 80) :: TITLE ! General problem description INTEGER :: NO_SYS_GEOM_NODES ! Number of geometry nodes in the system INTEGER :: NO_ELEMENTS ! Number of elements in the system INTEGER :: NO_PARAMS_PER_NODE ! Number of parameters per analysis node INTEGER :: UNUSED_1 INTEGER :: SOL_SPACE_DIM ! Dimension of the solution space INTEGER :: UNUSED_2 INTEGER :: UNUSED_3 INTEGER :: UNUSED_4 INTEGER :: UNUSED_5 INTEGER :: RHS_FLAG ! Initial system column vector flag INTEGER :: NO_USER_REMARKS ! Number of user remark lines INTEGER :: NO_ROWS_B ! Number of rows in the element B matrix INTEGER :: NO_QUAD_PTS ! Number of quadrature points INTEGER :: UNUSED_6 INTEGER :: UNUSED_7 INTEGER :: DEBUG_MODE INTEGER :: NO_INT_PRPS_PER_NODE ! Number of integer properties per node INTEGER :: NO_REAL_PRPS_PER_NODE ! Number of real properties per node INTEGER :: NO_INT_PRPS_PER_EL ! Number of integer properties per el INTEGER :: NO_REAL_PRPS_PER_EL ! Number of real properties per element INTEGER :: NO_INT_PRPS_MISC ! Number of miscellaneous integer props INTEGER :: NO_REAL_PRPS_MISC ! Number of miscellaneous real props INTEGER :: NODAL_PRPS_FLAG ! Nodal properties flag INTEGER :: ELEMENTAL_PRPS_FLAG ! Elemental properties flag INTEGER :: NODAL_PRINT_FLAG ! Nodal parameter print flag INTEGER :: ELEMENTAL_PRINT_FLAG ! Elemental parameter print flag INTEGER :: POST_PROC_IO_UNIT1 ! I/O unit number used in post-processing INTEGER :: POST_PROC_IO_UNIT2 ! I/O unit number used in post-processing INTEGER :: NO_FLUX_COMPS ! Number of generalized flux components INTEGER :: EL_COL_FLAG ! Element column matrix flag INTEGER :: NO_INT_PRPS_PER_SEG ! Number of segment integer properties INTEGER :: NO_REAL_PRPS_PER_SEG ! Number of segment real properties INTEGER :: ADAPT_FLAG ! Type of adaptivity to be used ! 0 => no adaptivity ! 1 => p -adaptivity ! 2 => r -adaptivity ! 3 => rp-adaptivity INTEGER :: PTCH_DFN_FLAG ! SPR patch definition indicator ! 0 => adjacent elements only ! 1 => all surrounding elements INTEGER :: PLYNML_FLAG ! Polynomial for least squares fit ! 0 => complete for highest degree ! 1 => shape function ! 2 => complete for highest degree + 1 INTEGER :: MAX_ITERS ! Upper bound on adaptive iterations REAL (DP) :: PRCNT_CUTOFF ! Error percentage cutoff INTEGER :: NEXT_UNUSED_NODE END MODULE CONTROL_MODULE SUBROUTINE COUNT_CONSTRAINTS (NODE_LST, NO_SPTD_TYPES, CNSTRNT_CNT) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER, INTENT (IN) :: NO_SPTD_TYPES INTEGER, DIMENSION (NO_SPTD_TYPES), INTENT (OUT) :: CNSTRNT_CNT TYPE (NODE_LST_NDE_TYPE), POINTER :: TEMP_PTR INTEGER :: PARAM_CNT INTEGER :: SPTD_TYPE INTEGER :: CNSTRNT_TYPE INTEGER, DIMENSION (NO_SPTD_TYPES) :: COUNTER INTEGER, DIMENSION (NO_PARAMS_PER_NODE) :: NODAL_CNSTRNTS TEMP_PTR => NODE_LST % HDR % NEXT_PTR DO SPTD_TYPE = 1,NO_SPTD_TYPES COUNTER(SPTD_TYPE) = 0 END DO ! SPTD_TYPE DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (TEMP_PTR % CNSTRNT_INDICATOR > 0) THEN CALL EXTRACT_CONSTRAINTS (TEMP_PTR, NODAL_CNSTRNTS) DO PARAM_CNT = 1,NO_PARAMS_PER_NODE CNSTRNT_TYPE = NODAL_CNSTRNTS(PARAM_CNT) IF (CNSTRNT_TYPE > 0) COUNTER(CNSTRNT_TYPE) = COUNTER(CNSTRNT_TYPE) + 1 END DO ! PARAM_CNT ENDIF TEMP_PTR => TEMP_PTR % NEXT_PTR END DO WRITE (6,"(//,A,//)") & " C O N S T R A I N T E Q U A T I O N D A T A" DO SPTD_TYPE = 1,NO_SPTD_TYPES CNSTRNT_CNT(SPTD_TYPE) = COUNTER(SPTD_TYPE) / SPTD_TYPE END DO ! SPTD_TYPE END SUBROUTINE COUNT_CONSTRAINTS MODULE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE ! Node list declarations TYPE NODE_LST_NDE_TYPE INTEGER :: ID REAL (DP) :: X REAL (DP) :: Y INTEGER :: DOF_INFO INTEGER :: CNSTRNT_INDICATOR TYPE (NODE_LST_NDE_TYPE), POINTER :: PREV_PTR TYPE (NODE_LST_NDE_TYPE), POINTER :: NEXT_PTR END TYPE NODE_LST_NDE_TYPE TYPE NODE_LST_PTR_TYPE TYPE (NODE_LST_NDE_TYPE), POINTER :: HDR END TYPE NODE_LST_PTR_TYPE ! Edge list declarations TYPE NODE_SUBLST_NDE_TYPE TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: LPTR TYPE (NODE_LST_NDE_TYPE), POINTER :: NPTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: RPTR END TYPE NODE_SUBLST_NDE_TYPE TYPE EDGE_LST_NDE_TYPE INTEGER :: ID TYPE (NODE_SUBLST_NDE_TYPE) :: NS_HDR TYPE (EDGE_LST_NDE_TYPE), POINTER :: NEXT END TYPE EDGE_LST_NDE_TYPE TYPE EDGE_LST_PTR_TYPE TYPE (EDGE_LST_NDE_TYPE), POINTER :: HDR END TYPE EDGE_LST_PTR_TYPE ! Element list declarations TYPE ELEMENT_LST_NDE_TYPE INTEGER :: ID INTEGER, DIMENSION (1:4) :: CRNR_SEQ TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDG_1 TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ADJ_1 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDG_2 TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ADJ_2 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDG_3 TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ADJ_3 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDG_4 TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ADJ_4 TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: NEXT END TYPE ELEMENT_LST_NDE_TYPE TYPE ELEMENT_LST_PTR_TYPE TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: HDR END TYPE ELEMENT_LST_PTR_TYPE END MODULE DATA_STRUCTURES_MODULE MODULE DBL_PRECISION_MODULE INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND (15,307) END MODULE DBL_PRECISION_MODULE SUBROUTINE DEALLOC_ELMNT_ASSMBLY_DATA () USE ELMNT_ASSMBLY_MODULE IMPLICIT NONE DEALLOCATE (DOF_MAP) DEALLOCATE (NODE_CRDS) IF (ALLOCATED (NODE_PRPS_R)) DEALLOCATE (NODE_PRPS_R) IF (ALLOCATED (NODE_PRPS_I)) DEALLOCATE (NODE_PRPS_I) DEALLOCATE (AJ) DEALLOCATE (AJINV) DEALLOCATE (B_MTRX) DEALLOCATE (BODY) DEALLOCATE (DGH) DEALLOCATE (DLG) DEALLOCATE (DLH) DEALLOCATE (E_MTRX) DEALLOCATE (EB_MTRX) DEALLOCATE (EL_SQR_MTRX) IF (ALLOCATED (EL_SRC_VCTR)) DEALLOCATE (EL_SRC_VCTR) DEALLOCATE (G) DEALLOCATE (H) DEALLOCATE (HINTG) DEALLOCATE (STRAIN) DEALLOCATE (STRAN0) DEALLOCATE (STRESS) DEALLOCATE (PTR) DEALLOCATE (WTR) DEALLOCATE (PTS) DEALLOCATE (WTS) END SUBROUTINE DEALLOC_ELMNT_ASSMBLY_DATA SUBROUTINE DEALLOC_ERROR_EST_DATA () USE ERROR_EST_MODULE IMPLICIT NONE DEALLOCATE (DOF_MAP) DEALLOCATE (NODE_CRDS) DEALLOCATE (PTR) DEALLOCATE (WTR) DEALLOCATE (PTS) DEALLOCATE (WTS) DEALLOCATE (AJ) DEALLOCATE (AJINV) DEALLOCATE (DGH) DEALLOCATE (DLH) DEALLOCATE (EL_ERROR) DEALLOCATE (EL_SOL_VCTR) DEALLOCATE (EL_SPR_NDL_DRVS) DEALLOCATE (H) DEALLOCATE (SGMA_STR) DEALLOCATE (SGMA_HAT) END SUBROUTINE DEALLOC_ERROR_EST_DATA SUBROUTINE DEALLOC_SYSTEM_EQNS_DATA () USE SYSTEM_EQNS_MODULE IMPLICIT NONE DEALLOCATE (SYS_COL_HGHTS) DEALLOCATE (DIAG_LCTNS) DEALLOCATE (SYS_SQR_MTRX) DEALLOCATE (SYS_SRC_VCTR) DEALLOCATE (SYS_SOL_VCTR) END SUBROUTINE DEALLOC_SYSTEM_EQNS_DATA SUBROUTINE DERSHAFN (NNOSI, NURFN, NSRFN, MCORD, CDRPT, SFVAL, DEVAL) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, DIMENSION (12), INTENT (IN) :: NNOSI INTEGER, INTENT (IN) :: NURFN INTEGER, INTENT (IN) :: NSRFN INTEGER, INTENT (IN) :: MCORD REAL (DP), DIMENSION (3), INTENT (IN) :: CDRPT REAL (DP), INTENT (OUT) :: SFVAL REAL (DP), DIMENSION (3), INTENT (OUT) :: DEVAL REAL (DP) :: POLI1 REAL (DP), DIMENSION (3) :: POLI2 REAL (DP), DIMENSION (3,8) :: VERCD INTEGER, DIMENSION (3,8) :: NSVEN REAL (DP), DIMENSION (3) :: CDRFN INTEGER, DIMENSION (3,12) :: NSOPV REAL (DP), DIMENSION (3) :: FRSID INTEGER, DIMENSION (3) :: MNODE INTEGER :: ICORD REAL (DP) :: CPNUL INTEGER :: NSIDE INTEGER :: INODE REAL (DP), DIMENSION (3,300) :: CDNIS INTEGER :: NOPV1 INTEGER :: NOPV2 INTEGER :: ISRFN REAL (DP) :: PLAN2 REAL (DP) :: POLI3 INTEGER :: ICOR1 REAL (DP) :: DPOL1 INTEGER :: ICOR2 REAL (DP) :: DPOL2 INTEGER :: INOD1 REAL (DP) :: DETP2 INTEGER :: INOD2 REAL (DP) :: DPLA2 REAL (DP) :: DPOL3 REAL (DP) :: DETP3 VERCD(1,1:8) = (/-1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0/) VERCD(2,1:8) = (/-1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0/) VERCD(3,1:8) = (/-1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0/) NSVEN(1,1:8) = (/ 1, 1, 3, 3, 5, 5, 7, 7/) NSVEN(2,1:8) = (/ 4, 2, 2, 4, 8, 6, 6, 8/) NSVEN(3,1:8) = (/ 9, 10, 11, 12, 9, 10, 11, 12/) NSOPV(1,1:12) = (/ 7, 8, 5, 6, 3, 4, 1, 2, 3, 4, 1, 2/) NSOPV(2,1:12) = (/ 8, 5, 6, 7, 4, 1, 2, 3, 7, 8, 5, 6/) NSOPV(3,1:12) = (/-1, -2, 1, 2, -1, -2, 1, 2, -3, -3, -3, -3/) POLI1 = 1.d0 IF (NURFN == 0) THEN DO ICORD = 1,MCORD POLI1 = POLI1*(CDRPT(ICORD)+VERCD(ICORD,NSRFN)) / & (2.0 * VERCD(ICORD,NSRFN)) END DO CPNUL = 1.d0 POLI2(1) = 0.d0 POLI2(2) = 0.d0 POLI2(3) = 0.d0 DO ICORD = 1,MCORD NSIDE = NSVEN(ICORD,NSRFN) MNODE(ICORD) = NNOSI(NSIDE) - 2 IF (MNODE(ICORD) > 0) THEN CPNUL = CPNUL - 1.d0 POLI2(ICORD) = 1.d0 FRSID(ICORD) = 2.d0 / (NNOSI(NSIDE)-1) DO INODE = 1,MNODE(ICORD) CDNIS(ICORD,INODE) = -1.d0 + FRSID(ICORD) * INODE POLI2(ICORD) = POLI2(ICORD) * (CDRPT(ICORD) - & CDNIS(ICORD,INODE))/(VERCD(ICORD,NSRFN) - & CDNIS(ICORD,INODE)) END DO END IF END DO SFVAL = POLI1 * (POLI2(1) + POLI2(2) + POLI2(3) + CPNUL) ELSE NOPV1 = NSOPV(1,NSRFN) NOPV2 = NSOPV(2,NSRFN) ISRFN = ABS (NSOPV(3,NSRFN)) FRSID(1) = 2.0 / (NNOSI(NSRFN) - 1) CDRFN(1) = -VERCD(1,NOPV1) CDRFN(2) = -VERCD(2,NOPV1) CDRFN(3) = -VERCD(3,NOPV1) CDRFN(ISRFN) = (1.d0-FRSID(1)*NURFN)*NSOPV(3,NSRFN)/ISRFN DO ICORD = 1,MCORD POLI1 = POLI1*(CDRPT(ICORD)-VERCD(ICORD,NOPV1)) / & (CDRFN(ICORD)-VERCD(ICORD,NOPV1)) END DO PLAN2 = (CDRPT(ISRFN)-VERCD(ISRFN,NOPV2)) / & (CDRFN(ISRFN)-VERCD(ISRFN,NOPV2)) POLI3 = 1.d0 MNODE(1) = NNOSI(NSRFN)-2 IF (MNODE(1) > 0) THEN DO INODE = 1,MNODE(1) CDNIS(1,INODE) = -1.d0+FRSID(1)*INODE IF (ABS (CDNIS(1,INODE)-CDRFN(ISRFN)) > 0.0001) THEN POLI3 = POLI3*(CDRPT(ISRFN)-CDNIS(1,INODE)) / & (CDRFN(ISRFN)-CDNIS(1,INODE)) END IF END DO END IF SFVAL = POLI1*PLAN2*POLI3 ENDIF DO ICOR1 = 1,MCORD IF (NURFN == 0) THEN DPOL1 = POLI2(1)+POLI2(2)+POLI2(3)+CPNUL DO ICOR2 = 1,MCORD IF (ICOR2 /= ICOR1) THEN DPOL1 = DPOL1*(CDRPT(ICOR2)+VERCD(ICOR2,NSRFN)) / & (2.0*VERCD(ICOR2,NSRFN)) ELSE DPOL1 = DPOL1/(2.0*VERCD(ICOR2,NSRFN)) END IF END DO DPOL2 = 0.d0 DO INOD1 = 1,MNODE(ICOR1) DETP2 = 1.d0 DO INOD2 = 1,MNODE(ICOR1) IF (INOD2 /= INOD1) THEN DETP2 = DETP2*(CDRPT(ICOR1)-CDNIS(ICOR1,INOD2)) / & (VERCD(ICOR1,NSRFN)-CDNIS(ICOR1,INOD2)) ELSE DETP2 = DETP2/(VERCD(ICOR1,NSRFN)-CDNIS(ICOR1,INOD2)) END IF END DO DPOL2 = DPOL2+DETP2 END DO DPOL2 = DPOL2*POLI1 DEVAL(ICOR1) = DPOL1+DPOL2 ELSE DPOL1 = POLI3*PLAN2 DO ICOR2 = 1,MCORD IF (ICOR2 /= ICOR1) THEN DPOL1 = DPOL1*(CDRPT(ICOR2)-VERCD(ICOR2,NOPV1)) / & (CDRFN(ICOR2)-VERCD(ICOR2,NOPV1)) ELSE DPOL1 = DPOL1/(CDRFN(ICOR2)-VERCD(ICOR2,NOPV1)) END IF END DO DPLA2 = 0.d0 DPOL3 = 0.d0 IF (ICOR1 == ISRFN) THEN DPLA2 = POLI1*POLI3/(CDRFN(ISRFN)-VERCD(ISRFN,NOPV2)) DO INOD1 = 1,MNODE(1) IF (ABS(CDNIS(1,INOD1)-CDRFN(ISRFN)) > 0.0001) THEN DETP3 = 1.d0 DO INOD2 = 1,MNODE(1) IF (ABS(CDNIS(1,INOD2)-CDRFN(ISRFN)) > 0.0001) THEN IF (INOD2 /= INOD1) THEN DETP3 = DETP3*(CDRPT(ISRFN)-CDNIS(1,INOD2)) / & (CDRFN(ISRFN)-CDNIS(1,INOD2)) ELSE DETP3 = DETP3/(CDRFN(ISRFN)-CDNIS(1,INOD2)) END IF END IF END DO DPOL3 = DPOL3+DETP3 END IF END DO DPOL3 = DPOL3*POLI1*PLAN2 END IF DEVAL(ICOR1) = DPOL1+DPLA2+DPOL3 END IF END DO END SUBROUTINE DERSHAFN SUBROUTINE DETERMINE_M_AND_N (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS,M,N) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS INTEGER, INTENT (OUT) :: M INTEGER, INTENT (OUT) :: N INTEGER :: INDEX TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: TEMP_PTR M = 0 N = 0 DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN TEMP_PTR => SPR_DATA_LST % HDR %NEXT DO IF (TEMP_PTR % EL_ID == PTCH_IDS(INDEX)) EXIT TEMP_PTR => TEMP_PTR % NEXT END DO M = M + TEMP_PTR % NO_SP N = MAX (N,TEMP_PTR % HGHST_ORDR) END IF END DO END SUBROUTINE DETERMINE_M_AND_N SUBROUTINE DETERMINE_PTCH_BNDS (MAX_PTCH_ELS,PTCH_IDS,EL_LST,XMIN,XMAX, & YMIN,YMAX) USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST REAL (DP), INTENT (OUT) :: XMIN REAL (DP), INTENT (OUT) :: XMAX REAL (DP), INTENT (OUT) :: YMIN REAL (DP), INTENT (OUT) :: YMAX INTEGER :: INDEX TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: NS_HDR_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: EDG_NDE_PTR XMIN = HUGE (XMIN) XMAX = -HUGE (XMAX) YMIN = HUGE (YMIN) YMAX = -HUGE (YMAX) DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN EL_PTR => EL_LST % HDR % NEXT DO IF (EL_PTR % ID == PTCH_IDS(INDEX)) EXIT EL_PTR => EL_PTR % NEXT END DO NS_HDR_PTR => EL_PTR % EDG_1 % NS_HDR EDG_NDE_PTR => NS_HDR_PTR % RPTR DO IF (ASSOCIATED (EDG_NDE_PTR,NS_HDR_PTR)) EXIT IF (EDG_NDE_PTR % NPTR % X < XMIN) THEN XMIN = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % X > XMAX) THEN XMAX = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % Y < YMIN) THEN YMIN = EDG_NDE_PTR % NPTR % Y END IF IF (EDG_NDE_PTR % NPTR % Y > YMAX) THEN YMAX = EDG_NDE_PTR % NPTR % Y END IF EDG_NDE_PTR => EDG_NDE_PTR % RPTR END DO NS_HDR_PTR => EL_PTR % EDG_2 % NS_HDR EDG_NDE_PTR => NS_HDR_PTR % RPTR DO IF (ASSOCIATED (EDG_NDE_PTR,NS_HDR_PTR)) EXIT IF (EDG_NDE_PTR % NPTR % X < XMIN) THEN XMIN = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % X > XMAX) THEN XMAX = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % Y < YMIN) THEN YMIN = EDG_NDE_PTR % NPTR % Y END IF IF (EDG_NDE_PTR % NPTR % Y > YMAX) THEN YMAX = EDG_NDE_PTR % NPTR % Y END IF EDG_NDE_PTR => EDG_NDE_PTR % RPTR END DO NS_HDR_PTR => EL_PTR % EDG_3 % NS_HDR EDG_NDE_PTR => NS_HDR_PTR % RPTR DO IF (ASSOCIATED (EDG_NDE_PTR,NS_HDR_PTR)) EXIT IF (EDG_NDE_PTR % NPTR % X < XMIN) THEN XMIN = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % X > XMAX) THEN XMAX = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % Y < YMIN) THEN YMIN = EDG_NDE_PTR % NPTR % Y END IF IF (EDG_NDE_PTR % NPTR % Y > YMAX) THEN YMAX = EDG_NDE_PTR % NPTR % Y END IF EDG_NDE_PTR => EDG_NDE_PTR % RPTR END DO NS_HDR_PTR => EL_PTR % EDG_4 % NS_HDR EDG_NDE_PTR => NS_HDR_PTR % RPTR DO IF (ASSOCIATED (EDG_NDE_PTR,NS_HDR_PTR)) EXIT IF (EDG_NDE_PTR % NPTR % X < XMIN) THEN XMIN = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % X > XMAX) THEN XMAX = EDG_NDE_PTR % NPTR % X END IF IF (EDG_NDE_PTR % NPTR % Y < YMIN) THEN YMIN = EDG_NDE_PTR % NPTR % Y END IF IF (EDG_NDE_PTR % NPTR % Y > YMAX) THEN YMAX = EDG_NDE_PTR % NPTR % Y END IF EDG_NDE_PTR => EDG_NDE_PTR % RPTR END DO END IF END DO END SUBROUTINE DETERMINE_PTCH_BNDS FUNCTION DOTPRD (VCTR_1, VCTR_2, LNGTH) RESULT (VALUE) USE DBL_PRECISION_MODULE IMPLICIT NONE REAL (DP), DIMENSION (:), INTENT (IN) :: VCTR_1 REAL (DP), DIMENSION (:), INTENT (IN) :: VCTR_2 INTEGER, INTENT (IN) :: LNGTH REAL (DP) :: VALUE INTEGER :: INDEX VALUE = 0.d0 DO INDEX = 1,LNGTH VALUE = VALUE + VCTR_1(INDEX) * VCTR_2(INDEX) END DO END FUNCTION DOTPRD MODULE ELMNT_ASSMBLY_MODULE USE DBL_PRECISION_MODULE INTEGER :: NO_AN_NODES INTEGER :: NO_EL_DOFS REAL (DP) :: DTRMNT INTEGER :: NO_QP_R INTEGER :: NO_QP_S INTEGER, DIMENSION (:), ALLOCATABLE :: DOF_MAP REAL (DP), DIMENSION (:,:), ALLOCATABLE :: NODE_CRDS REAL (DP), DIMENSION (:,:), ALLOCATABLE :: NODE_PRPS_R REAL (DP), DIMENSION (:), ALLOCATABLE :: EL_PRPS_R INTEGER, DIMENSION (:,:), ALLOCATABLE :: NODE_PRPS_I INTEGER, DIMENSION (:), ALLOCATABLE :: EL_PRPS_I REAL (DP), DIMENSION (:,:), ALLOCATABLE :: AJ REAL (DP), DIMENSION (:,:), ALLOCATABLE :: AJINV REAL (DP), DIMENSION (:,:), ALLOCATABLE :: B_MTRX REAL (DP), DIMENSION (:), ALLOCATABLE :: BODY REAL (DP), DIMENSION (:,:), ALLOCATABLE :: DGH REAL (DP), DIMENSION (:,:), ALLOCATABLE :: DLG REAL (DP), DIMENSION (:,:), ALLOCATABLE :: DLH REAL (DP), DIMENSION (:,:), ALLOCATABLE :: E_MTRX REAL (DP), DIMENSION (:,:), ALLOCATABLE :: EB_MTRX REAL (DP), DIMENSION (:,:), ALLOCATABLE :: EL_SQR_MTRX REAL (DP), DIMENSION (:), ALLOCATABLE :: EL_SRC_VCTR REAL (DP), DIMENSION (:), ALLOCATABLE :: G REAL (DP), DIMENSION (:), ALLOCATABLE :: H REAL (DP), DIMENSION (:), ALLOCATABLE :: HINTG REAL (DP), DIMENSION (:), ALLOCATABLE :: STRAIN REAL (DP), DIMENSION (:), ALLOCATABLE :: STRAN0 REAL (DP), DIMENSION (:), ALLOCATABLE :: STRESS REAL (DP), DIMENSION (:), ALLOCATABLE :: PTR REAL (DP), DIMENSION (:), ALLOCATABLE :: WTR REAL (DP), DIMENSION (:), ALLOCATABLE :: PTS REAL (DP), DIMENSION (:), ALLOCATABLE :: WTS END MODULE ELMNT_ASSMBLY_MODULE MODULE ERROR_EST_MODULE USE DBL_PRECISION_MODULE INTEGER :: NO_EL_AN_NODES INTEGER :: NO_EL_DOFS REAL (DP) :: DTRMNT INTEGER :: NO_QP_R INTEGER :: NO_QP_S INTEGER, DIMENSION (:), ALLOCATABLE :: DOF_MAP REAL (DP), DIMENSION (:,:), ALLOCATABLE :: NODE_CRDS REAL (DP), DIMENSION (:), ALLOCATABLE :: PTR REAL (DP), DIMENSION (:), ALLOCATABLE :: WTR REAL (DP), DIMENSION (:), ALLOCATABLE :: PTS REAL (DP), DIMENSION (:), ALLOCATABLE :: WTS REAL (DP), DIMENSION (:,:), ALLOCATABLE :: AJ REAL (DP), DIMENSION (:,:), ALLOCATABLE :: AJINV REAL (DP), DIMENSION (:,:), ALLOCATABLE :: DGH REAL (DP), DIMENSION (:,:), ALLOCATABLE :: DLH REAL (DP), DIMENSION (:), ALLOCATABLE :: EL_ERROR REAL (DP), DIMENSION (:), ALLOCATABLE :: EL_SOL_VCTR REAL (DP), DIMENSION (:,:), ALLOCATABLE :: EL_SPR_NDL_DRVS REAL (DP), DIMENSION (:), ALLOCATABLE :: H REAL (DP), DIMENSION (:), ALLOCATABLE :: SGMA_STR REAL (DP), DIMENSION (:), ALLOCATABLE :: SGMA_HAT END MODULE ERROR_EST_MODULE SUBROUTINE ESTIMATE_ERROR (NO_SYS_DOFS,EL_LST,NODE_LST,SPR_DATA_LST,EL_ENORM, & CUTOFF) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE USE ERROR_EST_MODULE USE SPR_TYPES_MODULE USE SYSTEM_EQNS_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST REAL (DP), DIMENSION (NO_ELEMENTS), INTENT (OUT) :: EL_ENORM REAL (DP), INTENT (OUT) :: CUTOFF INTEGER :: GET_NO_ANALYSIS_NODES INTEGER :: NO_SYS_AN_NODES REAL (DP), DIMENSION (:,:), ALLOCATABLE :: SPRCNVRGNT_NDL_DRVS REAL (DP) :: GL_GNORM REAL (DP) :: GL_ENORM REAL (DP) :: EL_GNORM REAL (DP) :: GL_ERROR REAL (DP) :: RMS_NORM TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR INTEGER :: QPR INTEGER :: QPS REAL (DP) :: DETWT INTEGER :: CMPNT INTEGER :: I REAL (DP) :: SUM1 REAL (DP) :: SUM2 REAL (DP), DIMENSION (NO_ROWS_B) :: DIFF NO_SYS_AN_NODES = NO_SYS_DOFS/NO_PARAMS_PER_NODE ALLOCATE (SPRCNVRGNT_NDL_DRVS(NO_SYS_AN_NODES,NO_ROWS_B)) CALL RECOVER_SPRCNVRGNT_DERIVS (NO_SYS_AN_NODES,EL_LST,SPR_DATA_LST, & SPRCNVRGNT_NDL_DRVS) CALL OUTPUT_NDL_GRADS (NO_SYS_AN_NODES,NO_ROWS_B,NODE_LST, & SPRCNVRGNT_NDL_DRVS) GL_GNORM = 0.d0 GL_ENORM = 0.d0 EL_PTR => EL_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_PTR)) EXIT NO_EL_AN_NODES = GET_NO_ANALYSIS_NODES (EL_PTR) NO_EL_DOFS = NO_EL_AN_NODES * NO_PARAMS_PER_NODE CALL ALLOC_ERROR_EST_DATA (EL_PTR) CALL GET_SYS_DOFS_FOR_ELEMENT (EL_PTR,NO_EL_DOFS,DOF_MAP) CALL GET_NODAL_COORDINATES (EL_PTR,NO_EL_AN_NODES,NODE_CRDS) DO I = 1,NO_EL_DOFS EL_SOL_VCTR(I) = SYS_SOL_VCTR(DOF_MAP(I)) END DO CALL GET_EL_SPRCNVRGNT_NDL_DRVS (EL_PTR,NO_SYS_AN_NODES,NO_EL_AN_NODES, & SPRCNVRGNT_NDL_DRVS,EL_SPR_NDL_DRVS) CALL GET_1D_GAUSS_DATA (NO_QP_R,PTR,WTR) CALL GET_1D_GAUSS_DATA (NO_QP_S,PTS,WTS) EL_ENORM (EL_PTR % ID) = 0.d0 EL_GNORM = 0.d0 DO QPR = 1,NO_QP_R DO QPS = 1,NO_QP_S CALL COMPUTE_SRNDPTY_FNS_DRVS (EL_PTR,NO_EL_AN_NODES, & PTR(QPR),PTS(QPS),H,DLH) CALL COMPUTE_JACOBIAN (NO_EL_AN_NODES,DLH,NODE_CRDS,AJ) CALL COMPUTE_JCBN_INVRS_DTRMNT (AJ,AJINV,DTRMNT) CALL COMPUTE_GLOBAL_DRVS (NO_EL_AN_NODES,AJINV,DLH,DGH) DETWT = DTRMNT*WTR(QPR)*WTS(QPS) DO CMPNT = 1,NO_ROWS_B SGMA_STR (CMPNT) = 0.d0 SGMA_HAT (CMPNT) = 0.d0 DO I = 1,NO_EL_AN_NODES SGMA_STR (CMPNT) = SGMA_STR (CMPNT) + & EL_SPR_NDL_DRVS (I,CMPNT) * H(I) SGMA_HAT (CMPNT) = SGMA_HAT (CMPNT) + & EL_SOL_VCTR(I) * DGH (CMPNT,I) END DO DIFF (CMPNT) = SGMA_STR (CMPNT) - SGMA_HAT (CMPNT) END DO SUM1 = 0.d0 SUM2 = 0.d0 DO CMPNT = 1,NO_ROWS_B SUM1 = SUM1 + DIFF (CMPNT) * DIFF (CMPNT) SUM2 = SUM2 + SGMA_STR (CMPNT) * SGMA_STR (CMPNT) END DO EL_ENORM (EL_PTR % ID) = EL_ENORM (EL_PTR % ID) + SUM1 * DETWT EL_GNORM = EL_GNORM + SUM2 * DETWT END DO END DO GL_ENORM = GL_ENORM + EL_ENORM (EL_PTR % ID) GL_GNORM = GL_GNORM + EL_GNORM EL_ENORM (EL_PTR % ID) = SQRT (EL_ENORM (EL_PTR % ID)) CALL DEALLOC_ERROR_EST_DATA () EL_PTR => EL_PTR % NEXT END DO GL_ERROR = SQRT (GL_ENORM/(GL_GNORM+GL_ENORM)) * 100.d0 RMS_NORM = SQRT ((GL_GNORM+GL_ENORM)/NO_ELEMENTS) CUTOFF = RMS_NORM * (PRCNT_CUTOFF/100.d0) GL_ENORM = SQRT (GL_ENORM) GL_GNORM = SQRT (GL_GNORM) WRITE (6,"(//,A,//)") & " E R R O R E S T I M A T E D A T A" WRITE (6,"(2A,E15.8)") " GLOBAL PERCENT ERROR................", & "...................",GL_ERROR WRITE (6,"(2A,E15.8)") " AVERAGE (RMS) ELEMENT NORM..........", & "...................",RMS_NORM WRITE (6,"(2A,E15.8)") " CUTOFF GRADIENT ERROR NORM..........", & "...................",CUTOFF WRITE (6,"(2A,E15.8)") " GLOBAL GRADIENT ERROR NORM..........", & "...................",GL_ENORM WRITE (6,"(2A,E15.8)") " GLOBAL GRADIENT NORM................", & "...................",GL_GNORM WRITE (6,"(/)") WRITE (6,"(A/)") & " ELEMENT ID ERROR NORM RELATIVE ERROR (%) REFINE?" EL_PTR => EL_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_PTR)) EXIT IF (EL_ENORM(EL_PTR % ID) > CUTOFF) THEN WRITE (6,"(I18,6X,E15.8,4X,E15.8,6X,A)") & EL_PTR % ID, EL_ENORM(EL_PTR % ID), & (EL_ENORM(EL_PTR % ID)/RMS_NORM)*100.0, "**" ELSE WRITE (6,"(I18,6X,E15.8,4X,E15.8)") & EL_PTR % ID, EL_ENORM(EL_PTR % ID), & (EL_ENORM(EL_PTR % ID)/RMS_NORM)*100.0 END IF EL_PTR => EL_PTR % NEXT END DO CALL PURGE_SPR_DATA_LIST (SPR_DATA_LST) DEALLOCATE (SPRCNVRGNT_NDL_DRVS) END SUBROUTINE ESTIMATE_ERROR SUBROUTINE EVALUATE_COMPLETE_NDL_DRVS (N,DGR,MAX_PTCH_ELS,NO_SYS_AN_NODES, & EL_LST,PTCH_IDS,A,CNTS,SUMS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: DGR INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, INTENT (IN) :: NO_SYS_AN_NODES TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS REAL (DP), DIMENSION (N), INTENT (IN) :: A INTEGER, DIMENSION (NO_SYS_AN_NODES), INTENT (OUT) :: CNTS REAL (DP), DIMENSION (NO_SYS_AN_NODES,NO_ROWS_B), INTENT (OUT) :: SUMS INTEGER :: INDEX TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR LOGICAL :: RIGHT TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR REAL (DP), DIMENSION (N) :: P REAL (DP) :: X REAL (DP) :: Y REAL (DP) :: XMIN REAL (DP) :: XMAX REAL (DP) :: YMIN REAL (DP) :: YMAX INTEGER :: INDX LOGICAL, DIMENSION (NO_SYS_AN_NODES) :: PRVSLY_PRCSSD CALL DETERMINE_PTCH_BNDS (MAX_PTCH_ELS,PTCH_IDS,EL_LST,XMIN,XMAX,YMIN,YMAX) PRVSLY_PRCSSD = .FALSE. DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN EL_PTR => EL_LST % HDR % NEXT DO IF (EL_PTR % ID == PTCH_IDS(INDEX)) EXIT EL_PTR => EL_PTR % NEXT END DO HEAD_PTR => EL_PTR % EDG_1 % NS_HDR IF (EL_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0+2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0+2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_COMPLETE_AT_NODE (N,DGR,X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_2 % NS_HDR IF (EL_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0+2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0+2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_COMPLETE_AT_NODE (N,DGR,X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_3 % NS_HDR IF (EL_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0+2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0+2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_COMPLETE_AT_NODE (N,DGR,X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_4 % NS_HDR IF (EL_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == EL_PTR % CRNR_SEQ(1)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0+2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0+2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_COMPLETE_AT_NODE (N,DGR,X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO END IF END DO END SUBROUTINE EVALUATE_COMPLETE_NDL_DRVS SUBROUTINE EVALUATE_SHPFNCTN_NDL_DRVS (N,MAX_NODES_R,MAX_NODES_S,MAX_PTCH_ELS, & NO_SYS_AN_NODES,EL_LST,PTCH_IDS,A,CNTS,SUMS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: MAX_NODES_R INTEGER, INTENT (IN) :: MAX_NODES_S INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, INTENT (IN) :: NO_SYS_AN_NODES TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS REAL (DP), DIMENSION (N), INTENT (IN) :: A INTEGER, DIMENSION (NO_SYS_AN_NODES), INTENT (OUT) :: CNTS REAL (DP), DIMENSION (NO_SYS_AN_NODES,NO_ROWS_B), INTENT (OUT) :: SUMS INTEGER :: INDEX TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR LOGICAL :: RIGHT TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR REAL (DP), DIMENSION (N) :: P REAL (DP) :: X REAL (DP) :: Y REAL (DP) :: XMIN REAL (DP) :: XMAX REAL (DP) :: YMIN REAL (DP) :: YMAX INTEGER :: INDX LOGICAL, DIMENSION (NO_SYS_AN_NODES) :: PRVSLY_PRCSSD CALL DETERMINE_PTCH_BNDS (MAX_PTCH_ELS,PTCH_IDS,EL_LST,XMIN,XMAX,YMIN,YMAX) PRVSLY_PRCSSD = .FALSE. DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN EL_PTR => EL_LST % HDR % NEXT DO IF (EL_PTR % ID == PTCH_IDS(INDEX)) EXIT EL_PTR => EL_PTR % NEXT END DO HEAD_PTR => EL_PTR % EDG_1 % NS_HDR IF (EL_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0 + 2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0 + 2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_SHPFNCTN_AT_NODE (N,MAX_NODES_R,MAX_NODES_S, & X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_2 % NS_HDR IF (EL_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0 + 2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0 + 2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_SHPFNCTN_AT_NODE (N,MAX_NODES_R,MAX_NODES_S, & X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_3 % NS_HDR IF (EL_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0 + 2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0 + 2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_SHPFNCTN_AT_NODE (N,MAX_NODES_R,MAX_NODES_S, & X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_4 % NS_HDR IF (EL_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == EL_PTR % CRNR_SEQ(1)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN INDX = ABS (TEMP_PTR % NPTR % DOF_INFO) IF (.NOT. PRVSLY_PRCSSD(INDX)) THEN X = -1.d0 + 2.d0*((TEMP_PTR % NPTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0 + 2.d0*((TEMP_PTR % NPTR % Y - YMIN)/(YMAX-YMIN)) CALL CONSTRUCT_SHPFNCTN_AT_NODE (N,MAX_NODES_R,MAX_NODES_S, & X,Y,P) CALL UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX, & CNTS,SUMS) PRVSLY_PRCSSD(INDX) = .TRUE. END IF END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO END IF END DO END SUBROUTINE EVALUATE_SHPFNCTN_NDL_DRVS SUBROUTINE EXTRACT_CONSTRAINTS (NODE_PTR, CNSTRNTS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_NDE_TYPE), INTENT (IN) :: NODE_PTR INTEGER, DIMENSION (1:NO_PARAMS_PER_NODE), INTENT (OUT) :: CNSTRNTS INTEGER :: PARAM_CNT INTEGER :: OLD INTEGER :: NEW INTEGER :: SUM INTEGER :: CNSTRNT_TYPE SUM = 0 OLD = NODE_PTR % CNSTRNT_INDICATOR DO PARAM_CNT = 1,NO_PARAMS_PER_NODE NEW = OLD / 10 CNSTRNT_TYPE = OLD - (NEW * 10) CNSTRNTS((NO_PARAMS_PER_NODE+1)-PARAM_CNT) = CNSTRNT_TYPE SUM = SUM + CNSTRNT_TYPE * 10**(PARAM_CNT-1) OLD = NEW END DO ! PARAM_CNT IF (NODE_PTR % CNSTRNT_INDICATOR > SUM) WRITE (6,"(A,I10)") & "WARNING: NODAL CONSTRAINT INDICATOR NOT RIGHT JUSTIFIED FOR NODE", & NODE_PTR % ID END SUBROUTINE EXTRACT_CONSTRAINTS SUBROUTINE EXTRACT_INT_ELEMENT_PRPS (ELEMENT_PTR, EL_PRPS_I) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE PROPERTIES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, DIMENSION (NO_INT_PRPS_PER_EL), INTENT (OUT) :: EL_PRPS_I INTEGER :: INT_PRP_CNT IF (ELEMENTAL_PRPS_FLAG > 0) THEN DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_EL EL_PRPS_I(INT_PRP_CNT) = EL_PRPS_INT(1,INT_PRP_CNT) END DO ELSE DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_EL EL_PRPS_I(INT_PRP_CNT) = EL_PRPS_INT(ELEMENT_PTR % ID,INT_PRP_CNT) END DO ENDIF END SUBROUTINE EXTRACT_INT_ELEMENT_PRPS SUBROUTINE EXTRACT_INT_NODAL_PT_PRPS (ELEMENT_PTR, NO_AN_NODES, EL_NODE_PRPS_I) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE PROPERTIES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, INTENT (IN) :: NO_AN_NODES INTEGER, DIMENSION (NO_AN_NODES,NO_INT_PRPS_PER_NODE), INTENT (OUT) :: EL_NODE_PRPS_I INTEGER :: DOF_INFO TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR LOGICAL :: RIGHT INTEGER :: NODE_CNT INTEGER :: INT_PRP_CNT INTEGER :: COUNT IF (NODAL_PRPS_FLAG > 0) THEN DO NODE_CNT = 1,NO_AN_NODES DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(NODE_CNT,INT_PRP_CNT) = & NODE_PRPS_INT(1,INT_PRP_CNT) END DO END DO ELSE COUNT = 1 HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF ((TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) .OR. & (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(2))) THEN DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & NODE_PRPS_INT(TEMP_PTR % NPTR % ID,INT_PRP_CNT) END DO ELSE DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & HUGE (NODE_PRPS_INT(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF ((TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(2)) .OR. & (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(3))) THEN DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & NODE_PRPS_INT(TEMP_PTR % NPTR % ID,INT_PRP_CNT) END DO ELSE DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & HUGE (NODE_PRPS_INT(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF ((TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(3)) .OR. & (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(4))) THEN DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & NODE_PRPS_INT(TEMP_PTR % NPTR % ID,INT_PRP_CNT) END DO ELSE DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & HUGE (NODE_PRPS_INT(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) EXIT IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(3)) THEN DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & NODE_PRPS_INT(TEMP_PTR % NPTR % ID,INT_PRP_CNT) END DO ELSE DO INT_PRP_CNT = 1,NO_INT_PRPS_PER_NODE EL_NODE_PRPS_I(COUNT,INT_PRP_CNT) = & HUGE (NODE_PRPS_INT(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO ! nodal properties interpolation should probably go here ENDIF END SUBROUTINE EXTRACT_INT_NODAL_PT_PRPS SUBROUTINE EXTRACT_REAL_ELEMENT_PRPS (ELEMENT_PTR, EL_PRPS_R) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE PROPERTIES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR REAL (DP), DIMENSION (NO_REAL_PRPS_PER_EL), INTENT (OUT) :: EL_PRPS_R INTEGER :: REAL_PRP_CNT IF (ELEMENTAL_PRPS_FLAG > 0) THEN DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_EL EL_PRPS_R(REAL_PRP_CNT) = EL_PRPS_REAL(1,REAL_PRP_CNT) END DO ELSE DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_EL EL_PRPS_R(REAL_PRP_CNT) = EL_PRPS_REAL(ELEMENT_PTR % ID,REAL_PRP_CNT) END DO ENDIF END SUBROUTINE EXTRACT_REAL_ELEMENT_PRPS SUBROUTINE EXTRACT_REAL_NODAL_PT_PRPS (ELEMENT_PTR, NO_AN_NODES, EL_NODE_PRPS_R) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE PROPERTIES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, INTENT (IN) :: NO_AN_NODES REAL (DP), DIMENSION (NO_AN_NODES,NO_REAL_PRPS_PER_NODE), & INTENT (OUT) :: EL_NODE_PRPS_R INTEGER :: DOF_INFO TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR LOGICAL :: RIGHT INTEGER :: NODE_CNT INTEGER :: REAL_PRP_CNT INTEGER :: COUNT IF (NODAL_PRPS_FLAG > 0) THEN DO NODE_CNT = 1,NO_AN_NODES DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(NODE_CNT,REAL_PRP_CNT) = & NODE_PRPS_REAL(1,REAL_PRP_CNT) END DO END DO ELSE COUNT = 1 HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF ((TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) .OR. & (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(2))) THEN DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & NODE_PRPS_REAL(TEMP_PTR % NPTR % ID,REAL_PRP_CNT) END DO ELSE DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & HUGE (NODE_PRPS_REAL(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF ((TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(2)) .OR. & (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(3))) THEN DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & NODE_PRPS_REAL(TEMP_PTR % NPTR % ID,REAL_PRP_CNT) END DO ELSE DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & HUGE (NODE_PRPS_REAL(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF ((TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(3)) .OR. & (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(4))) THEN DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & NODE_PRPS_REAL(TEMP_PTR % NPTR % ID,REAL_PRP_CNT) END DO ELSE DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & HUGE (NODE_PRPS_REAL(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) EXIT IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(3)) THEN DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & NODE_PRPS_REAL(TEMP_PTR % NPTR % ID,REAL_PRP_CNT) END DO ELSE DO REAL_PRP_CNT = 1,NO_REAL_PRPS_PER_NODE EL_NODE_PRPS_R(COUNT,REAL_PRP_CNT) = & HUGE (NODE_PRPS_REAL(1,1)) END DO ENDIF COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO ! nodal properties interpolation should probably go here ENDIF END SUBROUTINE EXTRACT_REAL_NODAL_PT_PRPS SUBROUTINE GENERATE_EL_SQR_MTRX (ELEMENT_PTR) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE ELMNT_ASSMBLY_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER :: QPR INTEGER :: QPS INTEGER :: I INTEGER :: J REAL (DP) :: DETWT REAL (DP) :: XK REAL (DP) :: YK REAL (DP), DIMENSION (SOL_SPACE_DIM) :: QP_XY XK = EL_PRPS_R (1) YK = EL_PRPS_R (2) CALL ZERO_OUT_MATRIX (NO_EL_DOFS, NO_EL_DOFS, EL_SQR_MTRX) ! Determine the quadrature data CALL GET_1D_GAUSS_DATA (NO_QP_R,PTR,WTR) CALL GET_1D_GAUSS_DATA (NO_QP_S,PTS,WTS) IF (POST_PROC_IO_UNIT1 > 0) WRITE (POST_PROC_IO_UNIT1) NO_QP_R*NO_QP_S DO QPR = 1,NO_QP_R DO QPS = 1,NO_QP_S CALL COMPUTE_SRNDPTY_FNS_DRVS (ELEMENT_PTR, NO_AN_NODES, & PTR(QPR), PTS(QPS), H, DLH) CALL MULTIPLY_MATRICES (1,NO_AN_NODES,SOL_SPACE_DIM,H,NODE_CRDS,QP_XY) CALL COMPUTE_JACOBIAN (NO_AN_NODES, DLH, NODE_CRDS, AJ) CALL COMPUTE_JCBN_INVRS_DTRMNT (AJ, AJINV, DTRMNT) CALL COMPUTE_GLOBAL_DRVS (NO_AN_NODES, AJINV, DLH, DGH) IF (POST_PROC_IO_UNIT1 > 0) WRITE (POST_PROC_IO_UNIT1) & QP_XY, DGH DETWT = DTRMNT*WTR(QPR)*WTS(QPS) DO J = 1,NO_AN_NODES DO I = 1,NO_AN_NODES EL_SQR_MTRX(I,J) = EL_SQR_MTRX(I,J)+DETWT* & (XK*DGH(1,I)*DGH(1,J) + YK*DGH(2,I)*DGH(2,J)) END DO END DO END DO END DO END SUBROUTINE GENERATE_EL_SQR_MTRX SUBROUTINE GENERATE_EL_SRC_VCTR () USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE ELMNT_ASSMBLY_MODULE IMPLICIT NONE END SUBROUTINE GENERATE_EL_SRC_VCTR SUBROUTINE GET_1D_GAUSS_DATA (NO_PTS, QUAD_PTS, QUAD_WTS) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_PTS REAL (DP), DIMENSION (NO_PTS), INTENT (OUT) :: QUAD_PTS REAL (DP), DIMENSION (NO_PTS), INTENT (OUT) :: QUAD_WTS IF (NO_PTS == 2) THEN QUAD_PTS(1) = 0.577350269189625764509149d0 QUAD_PTS(2) = -QUAD_PTS(1) QUAD_WTS(1) = 1.00000000000000000000000d0 QUAD_WTS(2) = 1.00000000000000000000000d0 ELSE IF (NO_PTS == 3) THEN QUAD_PTS(1) = 0.774596669241483377035853d0 QUAD_PTS(2) = 0.000000000000000000000000d0 QUAD_PTS(3) = -QUAD_PTS(1) QUAD_WTS(1) = 0.55555555555555555555556d0 QUAD_WTS(2) = 0.88888888888888888888889d0 QUAD_WTS(3) = QUAD_WTS(1) ELSE IF (NO_PTS == 4) THEN QUAD_PTS(1) = 0.861136311594052575223946d0 QUAD_PTS(2) = 0.339981043584856264802666d0 QUAD_PTS(3) = -QUAD_PTS(2) QUAD_PTS(4) = -QUAD_PTS(1) QUAD_WTS(1) = 0.34785484513745385737306d0 QUAD_WTS(2) = 0.65214515486254614262694d0 QUAD_WTS(3) = QUAD_WTS(2) QUAD_WTS(4) = QUAD_WTS(1) ELSE IF (NO_PTS == 5) THEN QUAD_PTS(1) = 0.906179845938663992797627d0 QUAD_PTS(2) = 0.538469310105683091036314d0 QUAD_PTS(3) = 0.000000000000000000000000d0 QUAD_PTS(4) = -QUAD_PTS(2) QUAD_PTS(5) = -QUAD_PTS(1) QUAD_WTS(1) = 0.23692688505618908751426d0 QUAD_WTS(2) = 0.47862867049936646804129d0 QUAD_WTS(3) = 0.56888888888888888888889d0 QUAD_WTS(4) = QUAD_WTS(2) QUAD_WTS(5) = QUAD_WTS(1) ELSE IF (NO_PTS == 6) THEN QUAD_PTS(1) = 0.932469514203152027812302d0 QUAD_PTS(2) = 0.661209386466264513661400d0 QUAD_PTS(3) = 0.238619186083196908630502d0 QUAD_PTS(4) = -QUAD_PTS(3) QUAD_PTS(5) = -QUAD_PTS(2) QUAD_PTS(6) = -QUAD_PTS(1) QUAD_WTS(1) = 0.17132449237917034504030d0 QUAD_WTS(2) = 0.36076157304813860756983d0 QUAD_WTS(3) = 0.46791393457269104738987d0 QUAD_WTS(4) = QUAD_WTS(3) QUAD_WTS(5) = QUAD_WTS(2) QUAD_WTS(6) = QUAD_WTS(1) ELSE IF (NO_PTS == 7) THEN QUAD_PTS(1) = 0.949107912342758524526190d0 QUAD_PTS(2) = 0.741531185599394439863865d0 QUAD_PTS(3) = 0.405845151377397166906607d0 QUAD_PTS(4) = 0.000000000000000000000000d0 QUAD_PTS(5) = -QUAD_PTS(3) QUAD_PTS(6) = -QUAD_PTS(2) QUAD_PTS(7) = -QUAD_PTS(1) QUAD_WTS(1) = 0.12948496616886969327061d0 QUAD_WTS(2) = 0.27970539148927666790147d0 QUAD_WTS(3) = 0.38183005050511894495037d0 QUAD_WTS(4) = 0.41795918367346938775510d0 QUAD_WTS(5) = QUAD_WTS(3) QUAD_WTS(6) = QUAD_WTS(2) QUAD_WTS(7) = QUAD_WTS(1) ELSE IF (NO_PTS == 8) THEN QUAD_PTS(1) = 0.960289856497536231683561d0 QUAD_PTS(2) = 0.796666477413626739591554d0 QUAD_PTS(3) = 0.525532409916328985817739d0 QUAD_PTS(4) = 0.183434642495649804939476d0 QUAD_PTS(5) = -QUAD_PTS(4) QUAD_PTS(6) = -QUAD_PTS(3) QUAD_PTS(7) = -QUAD_PTS(2) QUAD_PTS(8) = -QUAD_PTS(1) QUAD_WTS(1) = 0.10122853629037625915253d0 QUAD_WTS(2) = 0.22238103445337447054436d0 QUAD_WTS(3) = 0.31370664587788728733796d0 QUAD_WTS(4) = 0.36268378337836198296515d0 QUAD_WTS(5) = QUAD_WTS(4) QUAD_WTS(6) = QUAD_WTS(3) QUAD_WTS(7) = QUAD_WTS(2) QUAD_WTS(8) = QUAD_WTS(1) ELSE IF (NO_PTS == 9) THEN QUAD_PTS(1) = 0.968160239507626089835576d0 QUAD_PTS(2) = 0.836031107326635794299430d0 QUAD_PTS(3) = 0.613371432700590397308702d0 QUAD_PTS(4) = 0.324253423403808929038538d0 QUAD_PTS(5) = 0.000000000000000000000000d0 QUAD_PTS(6) = -QUAD_PTS(4) QUAD_PTS(7) = -QUAD_PTS(3) QUAD_PTS(8) = -QUAD_PTS(2) QUAD_PTS(9) = -QUAD_PTS(1) QUAD_WTS(1) = 0.08127438836157441197189d0 QUAD_WTS(2) = 0.18064816069485740405847d0 QUAD_WTS(3) = 0.26061069640293546231874d0 QUAD_WTS(4) = 0.31234707704000284006863d0 QUAD_WTS(5) = 0.33023935500125976316453d0 QUAD_WTS(6) = QUAD_WTS(4) QUAD_WTS(7) = QUAD_WTS(3) QUAD_WTS(8) = QUAD_WTS(2) QUAD_WTS(9) = QUAD_WTS(1) ELSE IF (NO_PTS == 10) THEN QUAD_PTS(1) = 0.973906528517171720077964d0 QUAD_PTS(2) = 0.865063366688984510732097d0 QUAD_PTS(3) = 0.679409568299024406234327d0 QUAD_PTS(4) = 0.433395394129247190799266d0 QUAD_PTS(5) = 0.148874338981631210884826d0 QUAD_PTS(6) = -QUAD_PTS(5) QUAD_PTS(7) = -QUAD_PTS(4) QUAD_PTS(8) = -QUAD_PTS(3) QUAD_PTS(9) = -QUAD_PTS(2) QUAD_PTS(10) = -QUAD_PTS(1) QUAD_WTS(1) = 0.06667134430868813759357d0 QUAD_WTS(2) = 0.14945134915058059314578d0 QUAD_WTS(3) = 0.21908636251598204399554d0 QUAD_WTS(4) = 0.26926671930999635509123d0 QUAD_WTS(5) = 0.29552422471475287017389d0 QUAD_WTS(6) = QUAD_WTS(5) QUAD_WTS(7) = QUAD_WTS(4) QUAD_WTS(8) = QUAD_WTS(3) QUAD_WTS(9) = QUAD_WTS(2) QUAD_WTS(10) = QUAD_WTS(1) ELSE IF (NO_PTS == 11) THEN QUAD_PTS(1) = 0.987228658146056992803938d0 QUAD_PTS(2) = 0.887062599768095299075158d0 QUAD_PTS(3) = 0.730152005574049324093416d0 QUAD_PTS(4) = 0.519096129206811815925726d0 QUAD_PTS(5) = 0.269543155952344972331532d0 QUAD_PTS(6) = 0.000000000000000000000000d0 QUAD_PTS(7) = -QUAD_PTS(5) QUAD_PTS(8) = -QUAD_PTS(4) QUAD_PTS(9) = -QUAD_PTS(3) QUAD_PTS(10) = -QUAD_PTS(2) QUAD_PTS(11) = -QUAD_PTS(1) QUAD_WTS(1) = 0.05566856711617366648275d0 QUAD_WTS(2) = 0.12558036946490462463469d0 QUAD_WTS(3) = 0.18629021092773425142610d0 QUAD_WTS(4) = 0.23319376459199047991852d0 QUAD_WTS(5) = 0.26280454451024666218069d0 QUAD_WTS(6) = 0.27292508677790063071448d0 QUAD_WTS(7) = QUAD_WTS(5) QUAD_WTS(8) = QUAD_WTS(4) QUAD_WTS(9) = QUAD_WTS(3) QUAD_WTS(10) = QUAD_WTS(2) QUAD_WTS(11) = QUAD_WTS(1) ELSE IF (NO_PTS == 12) THEN QUAD_PTS(1) = 0.981560634246719250690549d0 QUAD_PTS(2) = 0.904117256370474856678466d0 QUAD_PTS(3) = 0.769902674194304687036894d0 QUAD_PTS(4) = 0.587317954286617447296702d0 QUAD_PTS(5) = 0.367831498998180193752692d0 QUAD_PTS(6) = 0.125233408511468915472441d0 QUAD_PTS(7) = -QUAD_PTS(6) QUAD_PTS(8) = -QUAD_PTS(5) QUAD_PTS(9) = -QUAD_PTS(4) QUAD_PTS(10) = -QUAD_PTS(3) QUAD_PTS(11) = -QUAD_PTS(2) QUAD_PTS(12) = -QUAD_PTS(1) QUAD_WTS(1) = 0.04717533638651182719462d0 QUAD_WTS(2) = 0.10693932599531843096025d0 QUAD_WTS(3) = 0.16007832854334622633465d0 QUAD_WTS(4) = 0.20316742672306592174906d0 QUAD_WTS(5) = 0.23349253653835480876085d0 QUAD_WTS(6) = 0.24914704581340278500056d0 QUAD_WTS(7) = QUAD_WTS(6) QUAD_WTS(8) = QUAD_WTS(5) QUAD_WTS(9) = QUAD_WTS(4) QUAD_WTS(10) = QUAD_WTS(3) QUAD_WTS(11) = QUAD_WTS(2) QUAD_WTS(12) = QUAD_WTS(1) ELSE IF (NO_PTS == 13) THEN QUAD_PTS(1) = 0.984183054718588149472829d0 QUAD_PTS(2) = 0.917598399222977965206548d0 QUAD_PTS(3) = 0.801578090733309912794207d0 QUAD_PTS(4) = 0.642349339440340220643985d0 QUAD_PTS(5) = 0.448492751036446852877913d0 QUAD_PTS(6) = 0.230458315955134794065528d0 QUAD_PTS(7) = 0.000000000000000000000000d0 QUAD_PTS(8) = -QUAD_PTS(6) QUAD_PTS(9) = -QUAD_PTS(5) QUAD_PTS(10) = -QUAD_PTS(4) QUAD_PTS(11) = -QUAD_PTS(3) QUAD_PTS(12) = -QUAD_PTS(2) QUAD_PTS(13) = -QUAD_PTS(1) QUAD_WTS(1) = 0.04048400476531587952002d0 QUAD_WTS(2) = 0.09212149983772844791442d0 QUAD_WTS(3) = 0.13887351021978723846360d0 QUAD_WTS(4) = 0.17814598076194573828005d0 QUAD_WTS(5) = 0.20781604753688850231252d0 QUAD_WTS(6) = 0.22628318026289723841209d0 QUAD_WTS(7) = 0.23255155323087391019459d0 QUAD_WTS(8) = QUAD_WTS(6) QUAD_WTS(9) = QUAD_WTS(5) QUAD_WTS(10) = QUAD_WTS(4) QUAD_WTS(11) = QUAD_WTS(3) QUAD_WTS(12) = QUAD_WTS(2) QUAD_WTS(13) = QUAD_WTS(1) ELSE IF (NO_PTS == 14) THEN QUAD_PTS(1) = 0.986283808696812338841597d0 QUAD_PTS(2) = 0.928434883663573517336391d0 QUAD_PTS(3) = 0.827201315069764993189795d0 QUAD_PTS(4) = 0.687292904811685470148020d0 QUAD_PTS(5) = 0.515248636358154091965291d0 QUAD_PTS(6) = 0.319112368927889760435672d0 QUAD_PTS(7) = 0.108054948707343662066254d0 QUAD_PTS(8) = -QUAD_PTS(7) QUAD_PTS(9) = -QUAD_PTS(6) QUAD_PTS(10) = -QUAD_PTS(5) QUAD_PTS(11) = -QUAD_PTS(4) QUAD_PTS(12) = -QUAD_PTS(3) QUAD_PTS(13) = -QUAD_PTS(2) QUAD_PTS(14) = -QUAD_PTS(1) QUAD_WTS(1) = 0.03511946033175186303183d0 QUAD_WTS(2) = 0.08015808715976020980563d0 QUAD_WTS(3) = 0.12151857068790318468942d0 QUAD_WTS(4) = 0.15720316715819353456960d0 QUAD_WTS(5) = 0.18553839747793781374172d0 QUAD_WTS(6) = 0.20519846372129560396592d0 QUAD_WTS(7) = 0.21526385346315779019588d0 QUAD_WTS(8) = QUAD_WTS(7) QUAD_WTS(9) = QUAD_WTS(6) QUAD_WTS(10) = QUAD_WTS(5) QUAD_WTS(11) = QUAD_WTS(4) QUAD_WTS(12) = QUAD_WTS(3) QUAD_WTS(13) = QUAD_WTS(2) QUAD_WTS(14) = QUAD_WTS(1) ELSE IF (NO_PTS == 15) THEN QUAD_PTS(1) = 0.987992518020485428489566d0 QUAD_PTS(2) = 0.937273392400705904307759d0 QUAD_PTS(3) = 0.848206583410427216200648d0 QUAD_PTS(4) = 0.724417731360170047416186d0 QUAD_PTS(5) = 0.570972172608538847537227d0 QUAD_PTS(6) = 0.394151347077563369897207d0 QUAD_PTS(7) = 0.201194093997434522300628d0 QUAD_PTS(8) = 0.000000000000000000000000d0 QUAD_PTS(9) = -QUAD_PTS(7) QUAD_PTS(10) = -QUAD_PTS(6) QUAD_PTS(11) = -QUAD_PTS(5) QUAD_PTS(12) = -QUAD_PTS(4) QUAD_PTS(13) = -QUAD_PTS(3) QUAD_PTS(14) = -QUAD_PTS(2) QUAD_PTS(15) = -QUAD_PTS(1) QUAD_WTS(1) = 0.03075324199611726835463d0 QUAD_WTS(2) = 0.07036604748810812470927d0 QUAD_WTS(3) = 0.10715922046717193501187d0 QUAD_WTS(4) = 0.13957067792615431444781d0 QUAD_WTS(5) = 0.16626920581699393355320d0 QUAD_WTS(6) = 0.18616100001556221102680d0 QUAD_WTS(7) = 0.19843148532711157645612d0 QUAD_WTS(8) = 0.20257824192556127288062d0 QUAD_WTS(9) = QUAD_WTS(7) QUAD_WTS(10) = QUAD_WTS(6) QUAD_WTS(11) = QUAD_WTS(5) QUAD_WTS(12) = QUAD_WTS(4) QUAD_WTS(13) = QUAD_WTS(3) QUAD_WTS(14) = QUAD_WTS(2) QUAD_WTS(15) = QUAD_WTS(1) ELSE IF (NO_PTS == 16) THEN QUAD_PTS(1) = 0.989400934991649932596154d0 QUAD_PTS(2) = 0.944575023073232576077988d0 QUAD_PTS(3) = 0.865631202387831743880468d0 QUAD_PTS(4) = 0.755404408355003033895101d0 QUAD_PTS(5) = 0.617876244402643748446672d0 QUAD_PTS(6) = 0.458016777657227386342420d0 QUAD_PTS(7) = 0.281603550779258913230461d0 QUAD_PTS(8) = 0.095012509837637440185319d0 QUAD_PTS(9) = -QUAD_PTS(8) QUAD_PTS(10) = -QUAD_PTS(7) QUAD_PTS(11) = -QUAD_PTS(6) QUAD_PTS(12) = -QUAD_PTS(5) QUAD_PTS(13) = -QUAD_PTS(4) QUAD_PTS(14) = -QUAD_PTS(3) QUAD_PTS(15) = -QUAD_PTS(2) QUAD_PTS(16) = -QUAD_PTS(1) QUAD_WTS(1) = 0.02715245941175409485178d0 QUAD_WTS(2) = 0.06225352393864789286284d0 QUAD_WTS(3) = 0.09515851168249278480993d0 QUAD_WTS(4) = 0.12462897125553387205248d0 QUAD_WTS(5) = 0.14959598881657673208150d0 QUAD_WTS(6) = 0.16915651939500253818931d0 QUAD_WTS(7) = 0.18260341504492358886676d0 QUAD_WTS(8) = 0.18945061045506849628540d0 QUAD_WTS(9) = QUAD_WTS(8) QUAD_WTS(10) = QUAD_WTS(7) QUAD_WTS(11) = QUAD_WTS(6) QUAD_WTS(12) = QUAD_WTS(5) QUAD_WTS(13) = QUAD_WTS(4) QUAD_WTS(14) = QUAD_WTS(3) QUAD_WTS(15) = QUAD_WTS(2) QUAD_WTS(16) = QUAD_WTS(1) ELSE IF (NO_PTS == 17) THEN QUAD_PTS(1) = 0.990575475314417335675434d0 QUAD_PTS(2) = 0.950675521768767761222717d0 QUAD_PTS(3) = 0.880239153726985902122956d0 QUAD_PTS(4) = 0.781514003896801406925230d0 QUAD_PTS(5) = 0.657671159216690765850302d0 QUAD_PTS(6) = 0.512690537086476967886247d0 QUAD_PTS(7) = 0.351231763453876315297186d0 QUAD_PTS(8) = 0.178484181495847855850678d0 QUAD_PTS(9) = 0.000000000000000000000000d0 QUAD_PTS(10) = -QUAD_PTS(8) QUAD_PTS(11) = -QUAD_PTS(7) QUAD_PTS(12) = -QUAD_PTS(6) QUAD_PTS(13) = -QUAD_PTS(5) QUAD_PTS(14) = -QUAD_PTS(4) QUAD_PTS(15) = -QUAD_PTS(3) QUAD_PTS(16) = -QUAD_PTS(2) QUAD_PTS(17) = -QUAD_PTS(1) QUAD_WTS(1) = 0.02414830286854793196011d0 QUAD_WTS(2) = 0.05545952937398720112944d0 QUAD_WTS(3) = 0.08503614831717918088354d0 QUAD_WTS(4) = 0.11188384719340397109479d0 QUAD_WTS(5) = 0.13513636846852547328632d0 QUAD_WTS(6) = 0.15404576107681028808143d0 QUAD_WTS(7) = 0.16800410215645004450997d0 QUAD_WTS(8) = 0.17656270536699264632527d0 QUAD_WTS(9) = 0.17944647035620652545827d0 QUAD_WTS(10) = QUAD_WTS(8) QUAD_WTS(11) = QUAD_WTS(7) QUAD_WTS(12) = QUAD_WTS(6) QUAD_WTS(13) = QUAD_WTS(5) QUAD_WTS(14) = QUAD_WTS(4) QUAD_WTS(15) = QUAD_WTS(3) QUAD_WTS(16) = QUAD_WTS(2) QUAD_WTS(17) = QUAD_WTS(1) ELSE IF (NO_PTS == 18) THEN QUAD_PTS(1) = 0.991565168420930946730016d0 QUAD_PTS(2) = 0.955823949571397755181196d0 QUAD_PTS(3) = 0.892602466497555739206061d0 QUAD_PTS(4) = 0.803704958972523115682418d0 QUAD_PTS(5) = 0.691687043060353207874891d0 QUAD_PTS(6) = 0.559770831073947534607872d0 QUAD_PTS(7) = 0.411751161462842646035932d0 QUAD_PTS(8) = 0.251886225691505509588973d0 QUAD_PTS(9) = 0.084775013041735301242262d0 QUAD_PTS(10) = -QUAD_PTS(9) QUAD_PTS(11) = -QUAD_PTS(8) QUAD_PTS(12) = -QUAD_PTS(7) QUAD_PTS(13) = -QUAD_PTS(6) QUAD_PTS(14) = -QUAD_PTS(5) QUAD_PTS(15) = -QUAD_PTS(4) QUAD_PTS(16) = -QUAD_PTS(3) QUAD_PTS(17) = -QUAD_PTS(2) QUAD_PTS(18) = -QUAD_PTS(1) QUAD_WTS(1) = 0.02161601352648331031334d0 QUAD_WTS(2) = 0.04971454889496979645333d0 QUAD_WTS(3) = 0.07642573025488905652913d0 QUAD_WTS(4) = 0.10094204410628716556281d0 QUAD_WTS(5) = 0.12255520671147846018452d0 QUAD_WTS(6) = 0.14064291467065065120473d0 QUAD_WTS(7) = 0.15468467512626524492542d0 QUAD_WTS(8) = 0.16427648374583272298605d0 QUAD_WTS(9) = 0.16914238296314359184066d0 QUAD_WTS(10) = QUAD_WTS(9) QUAD_WTS(11) = QUAD_WTS(8) QUAD_WTS(12) = QUAD_WTS(7) QUAD_WTS(13) = QUAD_WTS(6) QUAD_WTS(14) = QUAD_WTS(5) QUAD_WTS(15) = QUAD_WTS(4) QUAD_WTS(16) = QUAD_WTS(3) QUAD_WTS(17) = QUAD_WTS(2) QUAD_WTS(18) = QUAD_WTS(1) ELSE IF (NO_PTS == 19) THEN QUAD_PTS(1) = 0.992406843843584403189018d0 QUAD_PTS(2) = 0.960208152134830030852779d0 QUAD_PTS(3) = 0.903155903614817901642661d0 QUAD_PTS(4) = 0.822714656537142824978923d0 QUAD_PTS(5) = 0.720966177335229378617096d0 QUAD_PTS(6) = 0.600545304661681023469638d0 QUAD_PTS(7) = 0.464570741375960945717267d0 QUAD_PTS(8) = 0.316564099963629831990117d0 QUAD_PTS(9) = 0.160358645640225375868096d0 QUAD_PTS(10) = 0.000000000000000000000000d0 QUAD_PTS(11) = -QUAD_PTS(9) QUAD_PTS(12) = -QUAD_PTS(8) QUAD_PTS(13) = -QUAD_PTS(7) QUAD_PTS(14) = -QUAD_PTS(6) QUAD_PTS(15) = -QUAD_PTS(5) QUAD_PTS(16) = -QUAD_PTS(4) QUAD_PTS(17) = -QUAD_PTS(3) QUAD_PTS(18) = -QUAD_PTS(2) QUAD_PTS(19) = -QUAD_PTS(1) QUAD_WTS(1) = 0.01946178822972647703631d0 QUAD_WTS(2) = 0.04481422676569960033284d0 QUAD_WTS(3) = 0.06904454273764122658071d0 QUAD_WTS(4) = 0.09149002162244999946446d0 QUAD_WTS(5) = 0.11156664554733399471602d0 QUAD_WTS(6) = 0.12875396253933622767552d0 QUAD_WTS(7) = 0.14260670217360661177575d0 QUAD_WTS(8) = 0.15276604206585966677886d0 QUAD_WTS(9) = 0.15896884339395434764996d0 QUAD_WTS(10) = 0.16105444984878369597916d0 QUAD_WTS(11) = QUAD_WTS(9) QUAD_WTS(12) = QUAD_WTS(8) QUAD_WTS(13) = QUAD_WTS(7) QUAD_WTS(14) = QUAD_WTS(6) QUAD_WTS(15) = QUAD_WTS(5) QUAD_WTS(16) = QUAD_WTS(4) QUAD_WTS(17) = QUAD_WTS(3) QUAD_WTS(18) = QUAD_WTS(2) QUAD_WTS(19) = QUAD_WTS(1) ELSE IF (NO_PTS == 20) THEN QUAD_PTS(1) = 0.993128599185094924786122d0 QUAD_PTS(2) = 0.963971927277913791267666d0 QUAD_PTS(3) = 0.912234428251325905867753d0 QUAD_PTS(4) = 0.839116971822218823394529d0 QUAD_PTS(5) = 0.746331906460150792614305d0 QUAD_PTS(6) = 0.636053680726515025452837d0 QUAD_PTS(7) = 0.510867001950827098004364d0 QUAD_PTS(8) = 0.373706088715419560672548d0 QUAD_PTS(9) = 0.227785851141645078080496d0 QUAD_PTS(10) = 0.076526521133497333754640d0 QUAD_PTS(11) = -QUAD_PTS(10) QUAD_PTS(12) = -QUAD_PTS(9) QUAD_PTS(13) = -QUAD_PTS(8) QUAD_PTS(14) = -QUAD_PTS(7) QUAD_PTS(15) = -QUAD_PTS(6) QUAD_PTS(16) = -QUAD_PTS(5) QUAD_PTS(17) = -QUAD_PTS(4) QUAD_PTS(18) = -QUAD_PTS(3) QUAD_PTS(19) = -QUAD_PTS(2) QUAD_PTS(20) = -QUAD_PTS(1) QUAD_WTS(1) = 0.01761400713915211831186d0 QUAD_WTS(2) = 0.04060142980038694133104d0 QUAD_WTS(3) = 0.06267204833410906356951d0 QUAD_WTS(4) = 0.08327674157670474872476d0 QUAD_WTS(5) = 0.10193011981724043503675d0 QUAD_WTS(6) = 0.11819453196151841731238d0 QUAD_WTS(7) = 0.13168863844917662689849d0 QUAD_WTS(8) = 0.14209610931838205132930d0 QUAD_WTS(9) = 0.14917298647260374678783d0 QUAD_WTS(10) = 0.15275338713072585069808d0 QUAD_WTS(11) = QUAD_WTS(10) QUAD_WTS(12) = QUAD_WTS(9) QUAD_WTS(13) = QUAD_WTS(8) QUAD_WTS(14) = QUAD_WTS(7) QUAD_WTS(15) = QUAD_WTS(6) QUAD_WTS(16) = QUAD_WTS(5) QUAD_WTS(17) = QUAD_WTS(4) QUAD_WTS(18) = QUAD_WTS(3) QUAD_WTS(19) = QUAD_WTS(2) QUAD_WTS(20) = QUAD_WTS(1) ELSE WRITE (6,"(/A/)") & "QUADRATURE DATA NOT TABULATED!" END IF END SUBROUTINE GET_1D_GAUSS_DATA FUNCTION GET_EDGE_PTR (EDGE_LST, END_1, END_2) RESULT (EDGE_PTR) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST INTEGER, INTENT (IN) :: END_1 INTEGER, INTENT (IN) :: END_2 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDGE_PTR TYPE (EDGE_LST_NDE_TYPE), POINTER :: TEMP_PTR TEMP_PTR => EDGE_LST % HDR % NEXT DO IF (((TEMP_PTR % NS_HDR % LPTR % NPTR % ID .EQ. END_1) .AND. & (TEMP_PTR % NS_HDR % RPTR % NPTR % ID .EQ. END_2)) .OR. & ((TEMP_PTR % NS_HDR % LPTR % NPTR % ID .EQ. END_2) .AND. & (TEMP_PTR % NS_HDR % RPTR % NPTR % ID .EQ. END_1))) EXIT TEMP_PTR => TEMP_PTR % NEXT END DO EDGE_PTR => TEMP_PTR END FUNCTION GET_EDGE_PTR SUBROUTINE GET_EL_SPRCNVRGNT_NDL_DRVS (EL_PTR,NO_SYS_AN_NODES,NO_EL_AN_NODES, & SYS_SPR_NDL_DRVS,EL_SPR_NDL_DRVS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: EL_PTR INTEGER, INTENT (IN) :: NO_SYS_AN_NODES INTEGER, INTENT (IN) :: NO_EL_AN_NODES REAL (DP), DIMENSION (NO_SYS_AN_NODES,NO_ROWS_B), INTENT (IN) :: & SYS_SPR_NDL_DRVS REAL (DP), DIMENSION (NO_EL_AN_NODES,NO_ROWS_B), INTENT (OUT) :: & EL_SPR_NDL_DRVS INTEGER :: CNT INTEGER :: DOF_INFO TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR INTEGER :: CMPNT LOGICAL :: RIGHT CNT = 1 HEAD_PTR => EL_PTR % EDG_1 % NS_HDR IF (EL_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO CMPNT = 1,NO_ROWS_B EL_SPR_NDL_DRVS(CNT,CMPNT) = SYS_SPR_NDL_DRVS(ABS(DOF_INFO),CMPNT) END DO ! CMPNT CNT = CNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_2 % NS_HDR IF (EL_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO CMPNT = 1,NO_ROWS_B EL_SPR_NDL_DRVS(CNT,CMPNT) = SYS_SPR_NDL_DRVS(ABS(DOF_INFO),CMPNT) END DO ! CMPNT CNT = CNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_3 % NS_HDR IF (EL_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO CMPNT = 1,NO_ROWS_B EL_SPR_NDL_DRVS(CNT,CMPNT) = SYS_SPR_NDL_DRVS(ABS(DOF_INFO),CMPNT) END DO ! CMPNT CNT = CNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => EL_PTR % EDG_4 % NS_HDR IF (EL_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == EL_PTR % CRNR_SEQ(1)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO CMPNT = 1,NO_ROWS_B EL_SPR_NDL_DRVS(CNT,CMPNT) = SYS_SPR_NDL_DRVS(ABS(DOF_INFO),CMPNT) END DO ! CMPNT CNT = CNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO END SUBROUTINE GET_EL_SPRCNVRGNT_NDL_DRVS SUBROUTINE GET_ELEMENT_COLUMN_HEIGHTS (NO_EL_DOFS, DOF_MAP, EL_COLUMN_HEIGHTS) IMPLICIT NONE INTEGER, INTENT (IN) :: NO_EL_DOFS INTEGER, DIMENSION (1:NO_EL_DOFS), INTENT (IN) :: DOF_MAP INTEGER, DIMENSION (1:NO_EL_DOFS), INTENT (OUT) :: EL_COLUMN_HEIGHTS INTEGER :: MIN_SYS_DOF INTEGER :: EL_DOF INTEGER :: SYS_DOF EL_COLUMN_HEIGHTS = 0 MIN_SYS_DOF = MINVAL (DOF_MAP) DO EL_DOF = 1,NO_EL_DOFS SYS_DOF = DOF_MAP(EL_DOF) EL_COLUMN_HEIGHTS(EL_DOF) = SYS_DOF - MIN_SYS_DOF + 1 END DO END SUBROUTINE GET_ELEMENT_COLUMN_HEIGHTS SUBROUTINE GET_ELS_IN_PATCH (EL_LST,EL_PTR,DFN_FLAG,MAX_PTCH_ELS,PTCH_IDS) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: EL_PTR INTEGER, INTENT (IN) :: DFN_FLAG INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (OUT) :: PTCH_IDS INTEGER :: INDEX TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: TEMP_PTR INTEGER :: CNR PTCH_IDS = 0 PTCH_IDS(1) = EL_PTR % ID IF (DFN_FLAG == 0) THEN IF (ASSOCIATED (EL_PTR % ADJ_1)) PTCH_IDS(2) = EL_PTR % ADJ_1 % ID IF (ASSOCIATED (EL_PTR % ADJ_2)) PTCH_IDS(3) = EL_PTR % ADJ_2 % ID IF (ASSOCIATED (EL_PTR % ADJ_3)) PTCH_IDS(4) = EL_PTR % ADJ_3 % ID IF (ASSOCIATED (EL_PTR % ADJ_4)) PTCH_IDS(5) = EL_PTR % ADJ_4 % ID ELSE INDEX = 2 TEMP_PTR => EL_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (TEMP_PTR) .OR. INDEX > MAX_PTCH_ELS) EXIT IF (TEMP_PTR % ID /= PTCH_IDS(1)) THEN DO CNR = 1,4 IF ((TEMP_PTR % CRNR_SEQ(CNR) == EL_PTR % CRNR_SEQ(1)) .OR. & (TEMP_PTR % CRNR_SEQ(CNR) == EL_PTR % CRNR_SEQ(2)) .OR. & (TEMP_PTR % CRNR_SEQ(CNR) == EL_PTR % CRNR_SEQ(3)) .OR. & (TEMP_PTR % CRNR_SEQ(CNR) == EL_PTR % CRNR_SEQ(4))) THEN PTCH_IDS(INDEX) = TEMP_PTR % ID INDEX = INDEX + 1 EXIT END IF END DO END IF TEMP_PTR => TEMP_PTR % NEXT END DO END IF END SUBROUTINE GET_ELS_IN_PATCH FUNCTION GET_HGHST_PTCH_DGR (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS) RESULT (DGR) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS INTEGER :: DGR INTEGER :: INDEX INTEGER :: MAX_EDG_NDS INTEGER :: MAX_NDS_R INTEGER :: MAX_NDS_S TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: TEMP_PTR MAX_EDG_NDS = 0 DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN TEMP_PTR => SPR_DATA_LST % HDR %NEXT DO IF (TEMP_PTR % EL_ID == PTCH_IDS(INDEX)) EXIT TEMP_PTR => TEMP_PTR % NEXT END DO MAX_NDS_R = TEMP_PTR % MAX_NDS_R MAX_NDS_S = TEMP_PTR % MAX_NDS_S MAX_EDG_NDS = MAX (MAX_EDG_NDS,MAX_NDS_R,MAX_NDS_S) END IF END DO DGR = MAX_EDG_NDS - 1 END FUNCTION GET_HGHST_PTCH_DGR SUBROUTINE GET_MATRIX_TRANSPOSE (M, N, A, AT) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: M INTEGER, INTENT (IN) :: N REAL (DP), DIMENSION (M,N), INTENT (IN) :: A REAL (DP), DIMENSION (N,M), INTENT (OUT) :: AT INTEGER :: I INTEGER :: J DO I = 1,N DO J = 1,M AT(I,J) = A(J,I) END DO END DO END SUBROUTINE GET_MATRIX_TRANSPOSE SUBROUTINE GET_MAX_NODES_IN_R_S (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS, & MAX_NODES_R,MAX_NODES_S) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS INTEGER, INTENT (OUT) :: MAX_NODES_R INTEGER, INTENT (OUT) :: MAX_NODES_S INTEGER :: INDEX TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: TEMP_PTR MAX_NODES_R = 0 MAX_NODES_S = 0 DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN TEMP_PTR => SPR_DATA_LST % HDR %NEXT DO IF (TEMP_PTR % EL_ID == PTCH_IDS(INDEX)) EXIT TEMP_PTR => TEMP_PTR % NEXT END DO IF (TEMP_PTR % MAX_NDS_R > MAX_NODES_R) & MAX_NODES_R = TEMP_PTR % MAX_NDS_R IF (TEMP_PTR % MAX_NDS_S > MAX_NODES_S) & MAX_NODES_S = TEMP_PTR % MAX_NDS_S END IF END DO END SUBROUTINE GET_MAX_NODES_IN_R_S FUNCTION GET_NEW_EDGE_ID (EDGE_LST, END_1, END_2, NEXT_UNUSED_EDGE_ID) & RESULT (NEW_EDGE_ID) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST INTEGER, INTENT (IN) :: END_1 INTEGER, INTENT (IN) :: END_2 INTEGER, INTENT (INOUT) :: NEXT_UNUSED_EDGE_ID INTEGER :: NEW_EDGE_ID TYPE (EDGE_LST_NDE_TYPE), POINTER :: TEMP_PTR TEMP_PTR => EDGE_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (((TEMP_PTR % NS_HDR % LPTR % NPTR % ID .EQ. END_1) .AND. & (TEMP_PTR % NS_HDR % RPTR % NPTR % ID .EQ. END_2)) .OR. & ((TEMP_PTR % NS_HDR % LPTR % NPTR % ID .EQ. END_2) .AND. & (TEMP_PTR % NS_HDR % RPTR % NPTR % ID .EQ. END_1))) EXIT TEMP_PTR => TEMP_PTR % NEXT END DO IF (.NOT. ASSOCIATED (TEMP_PTR)) THEN NEW_EDGE_ID = NEXT_UNUSED_EDGE_ID NEXT_UNUSED_EDGE_ID = NEXT_UNUSED_EDGE_ID + 1 ELSE NEW_EDGE_ID = 0 ENDIF END FUNCTION GET_NEW_EDGE_ID FUNCTION GET_NO_ANALYSIS_NODES (ELEMENT_PTR) RESULT (NO_AN_NODES) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER :: NO_AN_NODES TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR NO_AN_NODES = 4 HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR TEMP_PTR => HEAD_PTR % RPTR % RPTR DO IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN NO_AN_NODES = NO_AN_NODES + 1 END IF TEMP_PTR => TEMP_PTR % RPTR END DO HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR TEMP_PTR => HEAD_PTR % RPTR % RPTR DO IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN NO_AN_NODES = NO_AN_NODES + 1 END IF TEMP_PTR => TEMP_PTR % RPTR END DO HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR TEMP_PTR => HEAD_PTR % RPTR % RPTR DO IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN NO_AN_NODES = NO_AN_NODES + 1 END IF TEMP_PTR => TEMP_PTR % RPTR END DO HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR TEMP_PTR => HEAD_PTR % RPTR % RPTR DO IF (ASSOCIATED (TEMP_PTR % RPTR, HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN NO_AN_NODES = NO_AN_NODES + 1 END IF TEMP_PTR => TEMP_PTR % RPTR END DO END FUNCTION GET_NO_ANALYSIS_NODES FUNCTION GET_NO_NODES_ON_EDGE (EDGE_PTR) RESULT (NO_NODES) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDGE_PTR INTEGER :: NO_NODES TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: NS_TEMP_PTR NS_TEMP_PTR => EDGE_PTR % NS_HDR % RPTR NO_NODES = 0 DO IF (ASSOCIATED (NS_TEMP_PTR, EDGE_PTR % NS_HDR)) EXIT IF (NS_TEMP_PTR % NPTR % DOF_INFO /= 0) & NO_NODES = NO_NODES + 1 NS_TEMP_PTR => NS_TEMP_PTR % RPTR END DO END FUNCTION GET_NO_NODES_ON_EDGE FUNCTION GET_NO_PTCH_GAUSS_PTS (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS) & RESULT (NO_PTCH_GPS) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS INTEGER :: NO_PTCH_GPS INTEGER :: INDEX TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: TEMP_PTR NO_PTCH_GPS = 0 DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN TEMP_PTR => SPR_DATA_LST % HDR %NEXT DO IF (TEMP_PTR % EL_ID == PTCH_IDS(INDEX)) EXIT TEMP_PTR => TEMP_PTR % NEXT END DO NO_PTCH_GPS = NO_PTCH_GPS + TEMP_PTR % NO_SP END IF END DO END FUNCTION GET_NO_PTCH_GAUSS_PTS SUBROUTINE GET_NODAL_COORDINATES (ELEMENT_PTR, NO_NODES, NODE_CRDS) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, INTENT (IN) :: NO_NODES REAL (DP), DIMENSION (1:NO_NODES,1:SOL_SPACE_DIM), INTENT (OUT) :: NODE_CRDS INTEGER :: COUNT INTEGER :: DOF_INFO TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR LOGICAL :: RIGHT COUNT = 1 HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NODE_CRDS(COUNT,1) = TEMP_PTR % NPTR % X NODE_CRDS(COUNT,2) = TEMP_PTR % NPTR % Y COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NODE_CRDS(COUNT,1) = TEMP_PTR % NPTR % X NODE_CRDS(COUNT,2) = TEMP_PTR % NPTR % Y COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NODE_CRDS(COUNT,1) = TEMP_PTR % NPTR % X NODE_CRDS(COUNT,2) = TEMP_PTR % NPTR % Y COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN NODE_CRDS(COUNT,1) = TEMP_PTR % NPTR % X NODE_CRDS(COUNT,2) = TEMP_PTR % NPTR % Y COUNT = COUNT + 1 END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO END SUBROUTINE GET_NODAL_COORDINATES FUNCTION GET_NODE_PTR (NODE_LST, NODE_ID) RESULT (NODE_PTR) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER, INTENT (IN) :: NODE_ID TYPE (NODE_LST_NDE_TYPE), POINTER :: NODE_PTR TYPE (NODE_LST_NDE_TYPE), POINTER :: TEMP_PTR TEMP_PTR => NODE_LST % HDR % NEXT_PTR DO IF (TEMP_PTR % ID .EQ. NODE_ID) EXIT TEMP_PTR => TEMP_PTR % NEXT_PTR END DO NODE_PTR => TEMP_PTR END FUNCTION GET_NODE_PTR SUBROUTINE GET_SYS_DIAG_LCTNS (NO_SYS_DOFS, SYS_COL_HGHTS, DIAG_LCTNS) IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS INTEGER, DIMENSION (NO_SYS_DOFS), INTENT (IN) :: SYS_COL_HGHTS INTEGER, DIMENSION (NO_SYS_DOFS+1), INTENT (OUT) :: DIAG_LCTNS INTEGER :: DIAG_INDX INTEGER :: COL DIAG_INDX = 0 DIAG_LCTNS(1) = 0 DO COL = 1,NO_SYS_DOFS DIAG_INDX = DIAG_INDX + SYS_COL_HGHTS(COL) DIAG_LCTNS(COL+1) = DIAG_INDX END DO END SUBROUTINE GET_SYS_DIAG_LCTNS SUBROUTINE GET_SYS_DOFS_FOR_ELEMENT (ELEMENT_PTR, NO_EL_DOFS, DOF_MAP) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, INTENT (IN) :: NO_EL_DOFS INTEGER, DIMENSION (1:NO_EL_DOFS), INTENT (OUT) :: DOF_MAP INTEGER :: COUNT INTEGER :: DOF_INFO TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR INTEGER :: PARAM_CNT LOGICAL :: RIGHT COUNT = 1 HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO PARAM_CNT = 1,NO_PARAMS_PER_NODE DOF_MAP(COUNT) = (ABS(DOF_INFO) - 1) * NO_PARAMS_PER_NODE & + PARAM_CNT COUNT = COUNT + 1 END DO ! PARAM_CNT END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO PARAM_CNT = 1,NO_PARAMS_PER_NODE DOF_MAP(COUNT) = (ABS(DOF_INFO) - 1) * NO_PARAMS_PER_NODE & + PARAM_CNT COUNT = COUNT + 1 END DO ! PARAM_CNT END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR, HEAD_PTR)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO PARAM_CNT = 1,NO_PARAMS_PER_NODE DOF_MAP(COUNT) = (ABS(DOF_INFO) - 1) * NO_PARAMS_PER_NODE & + PARAM_CNT COUNT = COUNT + 1 END DO ! PARAM_CNT END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) EXIT DOF_INFO = TEMP_PTR % NPTR % DOF_INFO IF (DOF_INFO /= 0) THEN DO PARAM_CNT = 1,NO_PARAMS_PER_NODE DOF_MAP(COUNT) = (ABS(DOF_INFO) - 1) * NO_PARAMS_PER_NODE & + PARAM_CNT COUNT = COUNT + 1 END DO ! PARAM_CNT END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO END SUBROUTINE GET_SYS_DOFS_FOR_ELEMENT FUNCTION INITIALIZE_EDGE_LIST () RESULT (NEW_EDGE_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_PTR_TYPE) :: NEW_EDGE_LST ALLOCATE (NEW_EDGE_LST % HDR) NULLIFY (NEW_EDGE_LST % HDR % NEXT) END FUNCTION INITIALIZE_EDGE_LIST FUNCTION INITIALIZE_ELEMENT_LIST () RESULT (NEW_ELEMENT_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE) :: NEW_ELEMENT_LST ALLOCATE (NEW_ELEMENT_LST % HDR) NULLIFY (NEW_ELEMENT_LST % HDR % NEXT) END FUNCTION INITIALIZE_ELEMENT_LIST FUNCTION INITIALIZE_NODE_LIST () RESULT (NEW_NODE_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE) :: NEW_NODE_LST ALLOCATE (NEW_NODE_LST % HDR) NULLIFY (NEW_NODE_LST % HDR % PREV_PTR) NULLIFY (NEW_NODE_LST % HDR % NEXT_PTR) END FUNCTION INITIALIZE_NODE_LIST FUNCTION INITIALIZE_SPR_DATA_LIST () RESULT (NEW_SPR_DATA_LST) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE) :: NEW_SPR_DATA_LST ALLOCATE (NEW_SPR_DATA_LST % HDR) NULLIFY (NEW_SPR_DATA_LST % HDR % NEXT) END FUNCTION INITIALIZE_SPR_DATA_LIST SUBROUTINE INPUT_NODES_AND_ELEMENTS (NODE_LST, EDGE_LST, ELEMENT_LST) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST INTEGER :: NODE_CNT ! Loop control INTEGER :: NODE_ID ! Node identifier INTEGER :: CNSTRNT_INDICATOR ! Nodal constraint indicator REAL (DP) :: X ! X coordinate REAL (DP) :: Y ! Y coordinate INTEGER :: DOF_INFO = 0 INTEGER :: ELEMENT_CNT ! Loop control INTEGER :: ELEMENT_ID ! Element identifier INTEGER, DIMENSION (1:4) :: CORNERS ! Corner nodes of element INTEGER, DIMENSION (1:4) :: SIDES ! Midside nodes of element INTEGER :: INDEX ! Index used to read into arrays TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDGE_1 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDGE_2 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDGE_3 TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDGE_4 TYPE (EDGE_LST_NDE_TYPE), POINTER :: GET_EDGE_PTR INTEGER :: NEXT_UNUSED_EDGE_ID = 1 WRITE (6,"(//,A,//)") & " N O D E P O I N T D A T A" WRITE (6,"(A,/)") & " NODE ID CONSTRAINT INDICATOR X-COORDINATE Y-COORDINATE" ! Read in the nodal data and update the node list DO NODE_CNT = 1, NO_SYS_GEOM_NODES READ (5,"(2I10,2F10.4)") NODE_ID, & CNSTRNT_INDICATOR, & X, & Y WRITE (6,"(I14,I17,12X,F12.4,3X,F12.4)") NODE_ID, & CNSTRNT_INDICATOR, & X, & Y ! Insert node into node data structure CALL INSERT_NODE (NODE_LST, NODE_ID, X, Y, DOF_INFO, CNSTRNT_INDICATOR) END DO ! NODE_CNT WRITE (6,"(//,A,//)") & " E L E M E N T C O N N E C T I V I T Y D A T A" WRITE (6,"(A,/)") & " ELEMENT ID NODAL INCIDENCES" ! Read in the elements and update all the lists DO ELEMENT_CNT = 1, NO_ELEMENTS READ (5,"(9I5)") ELEMENT_ID, & (CORNERS(INDEX), INDEX=1,4), & (SIDES(INDEX), INDEX=1,4) WRITE (6,"(I16,6X,8I6)") ELEMENT_ID, & (CORNERS(INDEX), INDEX=1,4), & (SIDES(INDEX), INDEX=1,4) ! Update the edge list CALL UPDATE_EDGE_LIST (NODE_LST, EDGE_LST, CORNERS, SIDES, NEXT_UNUSED_EDGE_ID) ! Tag the analysis/geometry nodes CALL TAG_ANALYSIS_NODES (NODE_LST, CORNERS) ! Update the element list EDGE_1 => GET_EDGE_PTR (EDGE_LST, CORNERS(1), CORNERS(2)) EDGE_2 => GET_EDGE_PTR (EDGE_LST, CORNERS(2), CORNERS(3)) EDGE_3 => GET_EDGE_PTR (EDGE_LST, CORNERS(3), CORNERS(4)) EDGE_4 => GET_EDGE_PTR (EDGE_LST, CORNERS(4), CORNERS(1)) CALL INSERT_ELEMENT (ELEMENT_LST, ELEMENT_ID, CORNERS, EDGE_1, & EDGE_2, EDGE_3, EDGE_4) END DO ! ELEMENT_CNT CALL SET_ADJACENT_EL_PTRS (ELEMENT_LST) END SUBROUTINE INPUT_NODES_AND_ELEMENTS SUBROUTINE INPUT_PROPERTY_DATA () USE CONTROL_MODULE USE PROPERTIES_MODULE IMPLICIT NONE INTEGER :: MAX_PROPS INTEGER :: PROPERTY_CNT INTEGER :: J INTEGER :: K ! Process nodal point properties (if furnished) IF ((NO_INT_PRPS_PER_NODE > 0) .OR. (NO_REAL_PRPS_PER_NODE > 0)) THEN WRITE (6,"(//,A,//)") & " N O D A L P R O P E R T I E S" WRITE (6,"(A31)") " NODE PROPERTY VALUE" IF (NODAL_PRPS_FLAG == 1) THEN MAX_PROPS = 1 ELSE MAX_PROPS = NO_SYS_GEOM_NODES END IF IF (NO_INT_PRPS_PER_NODE > 0) THEN ALLOCATE (NODE_PRPS_INT(1:NO_SYS_GEOM_NODES,0:NO_INT_PRPS_PER_NODE)) DO PROPERTY_CNT = 1, MAX_PROPS READ (5,"(I10,7I10)") J, (NODE_PRPS_INT(J,K), & K=1,NO_INT_PRPS_PER_NODE) END DO ! PROPERTY_CNT WRITE (6,"(2I10,3X,I10)") ((J, K, NODE_PRPS_INT(J,K), & J=1,MAX_PROPS), K=1,NO_INT_PRPS_PER_NODE) WRITE (6,"(A31,/)") "END INTEGER PROPERTIES OF NODES" END IF IF (NO_REAL_PRPS_PER_NODE > 0) THEN ALLOCATE (NODE_PRPS_REAL(1:NO_SYS_GEOM_NODES,0:NO_REAL_PRPS_PER_NODE)) DO PROPERTY_CNT = 1, MAX_PROPS READ (5,"(I10,7F10.4)") J, (NODE_PRPS_REAL(J,K), & K=1,NO_REAL_PRPS_PER_NODE) END DO ! PROPERTY_CNT WRITE (6,"(2I10,3X,E12.5)") ((J, K, NODE_PRPS_REAL(J,K), & J=1,MAX_PROPS), K=1,NO_REAL_PRPS_PER_NODE) WRITE (6,"(A28,/)") "END REAL PROPERTIES OF NODES" END IF END IF ! Process element properties (if furnished) IF ((NO_INT_PRPS_PER_EL > 0) .OR. (NO_REAL_PRPS_PER_EL > 0)) THEN WRITE (6,"(//,A,//)") & " E L E M E N T P R O P E R T I E S" IF (ELEMENTAL_PRPS_FLAG == 1) THEN MAX_PROPS = 1 ELSE MAX_PROPS = NO_ELEMENTS END IF IF (NO_INT_PRPS_PER_EL > 0) THEN WRITE (6,"(A,/)") & " INTEGER PROPERTIES" WRITE (6,"(A,/)") & " ELEMENT ID PROPERTY ID VALUE" ALLOCATE (EL_PRPS_INT(1:NO_ELEMENTS,0:NO_INT_PRPS_PER_EL)) DO PROPERTY_CNT = 1, MAX_PROPS READ (5,"(I10,7I10)") J, (EL_PRPS_INT(J,K), & K=1,NO_INT_PRPS_PER_EL) END DO ! PROPERTY_CNT WRITE (6,"(I16,I25,I23)") ((J, K, EL_PRPS_INT(J,K), J=1,MAX_PROPS), & K=1,NO_INT_PRPS_PER_EL) WRITE (6,"(/)") END IF IF (NO_REAL_PRPS_PER_EL > 0) THEN WRITE (6,"(A,/)") & " REAL PROPERTIES" WRITE (6,"(A,/)") & " ELEMENT ID PROPERTY ID VALUE" ALLOCATE (EL_PRPS_REAL(1:NO_ELEMENTS,0:NO_REAL_PRPS_PER_EL)) DO PROPERTY_CNT = 1, MAX_PROPS READ (5,"(I10,7F10.4)") J, (EL_PRPS_REAL(J,K), & K=1,NO_REAL_PRPS_PER_EL) END DO ! PROPERTY_CNT WRITE (6,"(I16,I25,16X,E12.5)") ((J, K, EL_PRPS_REAL(J,K), & J=1,MAX_PROPS), K=1,NO_REAL_PRPS_PER_EL) WRITE (6,"(/)") END IF END IF ! Process segment properties (if furnished) ! IF ((NO_INT_PRPS_PER_SEG > 0) .OR. (NO_REAL_PRPS_PER_SEG > 0)) THEN ! WRITE (6,"(/,A31,/)") "*** SEGMENT PROPERTIES ***" ! WRITE (6,"(A28)") "SEGMENT PROPERTY VALUE" ! IF (ELEMENTAL_PRPS_FLAG == 1) THEN ! MAX_PROPS = 1 ! ELSE ! MAX_PROPS = NO_SEGMENTS ! END IF ! IF (NO_INT_PRPS_PER_SEG > 0) THEN ! ALLOCATE (SEG_PRPS_INT(0:NO_SEGMENTS,0:NO_INT_PRPS_PER_SEG)) ! DO PROPERTY_CNT = 1, MAX_PROPS ! READ (5,"(I10,7I10)") J, (SEG_PRPS_INT(J,K), & ! K=1,NO_INT_PRPS_PER_SEG) ! END DO ! PROPERTY_CNT ! WRITE (6,"(2I10,3X,I10)") ((J, K, SEG_PRPS_INT(J,K), J=1,MAX_PROPS), & ! K=1,NO_INT_PRPS_PER_SEG) ! WRITE (6,"(A34,/)") "END INTEGER PROPERTIES OF SEGMENTS" ! END IF ! IF (NO_REAL_PRPS_PER_SEG > 0) THEN ! ALLOCATE (SEG_PRPS_REAL(0:NO_SEGMENTS,0:NO_REAL_PRPS_PER_SEG)) ! DO PROPERTY_CNT = 1, MAX_PROPS ! READ (5,"(I10,7F10.4)") J, (SEG_PRPS_REAL(J,K), & ! K=1,NO_REAL_PRPS_PER_SEG) ! END DO ! PROPERTY_CNT ! WRITE (6,"(2I10,3X,E12.5)") ((J, K, SEG_PRPS_REAL(J,K), & ! J=1,MAX_PROPS), K=1,NO_REAL_PRPS_PER_SEG) ! WRITE (6,"(A31,/)") "END REAL PROPERTIES OF SEGMENTS" ! END IF ! END IF ! Process miscellaneous system properties (if furnished) IF ((NO_INT_PRPS_MISC > 0) .OR. (NO_REAL_PRPS_MISC > 0)) THEN WRITE (6,"(//,A,//)") & " M I S C E L L A N E O U S P R O P E R T I E S" WRITE (6,"(A18)") "PROPERTY VALUE" IF (NO_INT_PRPS_MISC > 0) THEN ALLOCATE (MISC_PRPS_INT(0:NO_INT_PRPS_MISC)) READ (5,"(8I10)") (MISC_PRPS_INT(K), K=1,NO_INT_PRPS_MISC) DO PROPERTY_CNT = 1, NO_INT_PRPS_MISC WRITE (6,"(I8,3X,I10)") PROPERTY_CNT, MISC_PRPS_INT(PROPERTY_CNT) END DO ! PROPERTY_CNT WRITE (6,"(A32,/)") "END INTEGER PROPERTIES OF SYSTEM" END IF IF (NO_REAL_PRPS_MISC > 0) THEN ALLOCATE (MISC_PRPS_REAL(0:NO_REAL_PRPS_MISC)) READ (5,"(8F10.4)") (MISC_PRPS_REAL(K), K=1,NO_REAL_PRPS_MISC) DO PROPERTY_CNT = 1, NO_REAL_PRPS_MISC WRITE (6,"(I8,3X,E13.5)") PROPERTY_CNT, MISC_PRPS_REAL(PROPERTY_CNT) END DO ! PROPERTY_CNT WRITE (6,"(A29,/)") "END REAL PROPERTIES OF SYSTEM" END IF END IF END SUBROUTINE INPUT_PROPERTY_DATA SUBROUTINE INPUT_RHS (NO_SYS_DOFS, RHS_VECTOR) USE DBL_PRECISION_MODULE USE CONTROL_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS REAL (DP), DIMENSION (1:NO_SYS_DOFS), INTENT (OUT) :: RHS_VECTOR REAL (DP), DIMENSION (1:NO_PARAMS_PER_NODE) :: RESULTANTS INTEGER :: SYS_DOF INTEGER :: NODE INTEGER :: PARAM REAL (DP) :: VALUE INTEGER :: INDEX INTEGER :: PARAMETER CALL ZERO_OUT_VECTOR (NO_SYS_DOFS, RHS_VECTOR) CALL ZERO_OUT_VECTOR (NO_PARAMS_PER_NODE, RESULTANTS) WRITE (6,"(/,A35)") "*** INITIAL FORCING VECTOR DATA ***" WRITE (6,"(A43)") " NODE PARAMETER VALUE EQUATION" DO SYS_DOF = 1, NO_SYS_DOFS READ (5,"(2I10,F15.4)") NODE, PARAM, VALUE INDEX = NO_PARAMS_PER_NODE * (NODE - 1) + PARAM RHS_VECTOR(INDEX) = VALUE RESULTANTS(PARAM) = RESULTANTS(PARAM) + VALUE WRITE (6,"(2I10,1X,1PE13.5,I8)") NODE, PARAM, VALUE, INDEX END DO WRITE (6,"(A12,/,A13)") "*RESULTANTS*","DOF SUM" DO PARAMETER = 1, NO_PARAMS_PER_NODE WRITE (6,"(I3,2X,1PE12.4)") PARAMETER, RESULTANTS(PARAMETER) END DO END SUBROUTINE INPUT_RHS SUBROUTINE INPUT_TYPE_1_CNSTRNT_DATA (CNSTRNT_CNT, TYPE_1_VALS, TYPE_1_DOFS) USE DBL_PRECISION_MODULE USE CONTROL_MODULE IMPLICIT NONE INTEGER, DIMENSION (1), INTENT (IN) :: CNSTRNT_CNT REAL (DP), DIMENSION (CNSTRNT_CNT(1)), INTENT (OUT) :: TYPE_1_VALS INTEGER, DIMENSION (CNSTRNT_CNT(1)), INTENT (OUT) :: TYPE_1_DOFS INTEGER :: EQN_INDEX INTEGER :: EQN_CNT INTEGER :: NODE1 INTEGER :: IPAR1 REAL (DP) :: A1 EQN_INDEX = 0 IF (CNSTRNT_CNT(1) /= 0) THEN WRITE (6,"(A)") & " CONSTRAINT TYPE # 1" WRITE (6,"(/,A,/)") & " EQUATION ID NODE ID PARAMETER VALUE" DO EQN_CNT = 1,CNSTRNT_CNT(1) EQN_INDEX = EQN_INDEX + 1 READ (5,"(2I10,F10.0)") NODE1, IPAR1, A1 WRITE (6,"(I16,I17,I16,9X,E12.5)") EQN_INDEX, NODE1, IPAR1, A1 TYPE_1_DOFS(EQN_INDEX) = NO_PARAMS_PER_NODE * (NODE1 - 1) + IPAR1 TYPE_1_VALS(EQN_INDEX) = A1 END DO ! EQN_CNT ENDIF END SUBROUTINE INPUT_TYPE_1_CNSTRNT_DATA SUBROUTINE INSERT_EDGE (EDGE_LST, EDGE_ID, END_1, MID, END_2) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST INTEGER, INTENT (IN) :: EDGE_ID TYPE (NODE_LST_NDE_TYPE), TARGET :: END_1 TYPE (NODE_LST_NDE_TYPE), POINTER :: MID TYPE (NODE_LST_NDE_TYPE), TARGET :: END_2 TYPE (EDGE_LST_NDE_TYPE), POINTER :: TRAIL_PTR TYPE (EDGE_LST_NDE_TYPE), POINTER :: TEMP_PTR ! Find location to insert new node TRAIL_PTR => EDGE_LST % HDR TEMP_PTR => TRAIL_PTR % NEXT DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (EDGE_ID < TEMP_PTR % ID) EXIT TRAIL_PTR => TEMP_PTR TEMP_PTR => TEMP_PTR % NEXT END DO ! Insert new node ALLOCATE (TRAIL_PTR % NEXT) TRAIL_PTR % NEXT % ID = EDGE_ID TRAIL_PTR % NEXT % NS_HDR % LPTR => TRAIL_PTR % NEXT % NS_HDR TRAIL_PTR % NEXT % NS_HDR % RPTR => TRAIL_PTR % NEXT % NS_HDR CALL NODE_SUBLST_INSERT (END_1, TRAIL_PTR % NEXT % NS_HDR) CALL NODE_SUBLST_INSERT (END_2, TRAIL_PTR % NEXT % NS_HDR % LPTR) IF (ASSOCIATED (MID)) THEN CALL NODE_SUBLST_INSERT (MID, TRAIL_PTR % NEXT % NS_HDR % RPTR) ENDIF TRAIL_PTR % NEXT % NEXT => TEMP_PTR END SUBROUTINE INSERT_EDGE SUBROUTINE INSERT_ELEMENT (ELEMENT_LST, ELEMENT_ID, CORNERS, EDGE_1, EDGE_2, & EDGE_3, EDGE_4) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST INTEGER, INTENT (IN) :: ELEMENT_ID INTEGER, DIMENSION (1:4), INTENT (IN) :: CORNERS TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDGE_1 TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDGE_2 TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDGE_3 TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDGE_4 TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: TRAIL_PTR TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: TEMP_PTR ! Find location to insert new node TRAIL_PTR => ELEMENT_LST % HDR TEMP_PTR => TRAIL_PTR % NEXT DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (ELEMENT_ID < TEMP_PTR % ID) EXIT TRAIL_PTR => TEMP_PTR TEMP_PTR => TEMP_PTR % NEXT END DO ! Insert new node ALLOCATE (TRAIL_PTR % NEXT) TRAIL_PTR % NEXT % ID = ELEMENT_ID TRAIL_PTR % NEXT % CRNR_SEQ (1) = CORNERS (1) TRAIL_PTR % NEXT % CRNR_SEQ (2) = CORNERS (2) TRAIL_PTR % NEXT % CRNR_SEQ (3) = CORNERS (3) TRAIL_PTR % NEXT % CRNR_SEQ (4) = CORNERS (4) TRAIL_PTR % NEXT % EDG_1 => EDGE_1 TRAIL_PTR % NEXT % EDG_2 => EDGE_2 TRAIL_PTR % NEXT % EDG_3 => EDGE_3 TRAIL_PTR % NEXT % EDG_4 => EDGE_4 NULLIFY (TRAIL_PTR % NEXT % ADJ_1) NULLIFY (TRAIL_PTR % NEXT % ADJ_2) NULLIFY (TRAIL_PTR % NEXT % ADJ_3) NULLIFY (TRAIL_PTR % NEXT % ADJ_4) TRAIL_PTR % NEXT % NEXT => TEMP_PTR END SUBROUTINE INSERT_ELEMENT SUBROUTINE INSERT_NODE (NODE_LST, NODE_ID, X, Y, DOF_INFO, CNSTRNT_INDICATOR) USE DBL_PRECISION_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER, INTENT (IN) :: NODE_ID REAL (DP), INTENT (IN) :: X REAL (DP), INTENT (IN) :: Y INTEGER, INTENT (IN) :: DOF_INFO INTEGER, INTENT (IN) :: CNSTRNT_INDICATOR TYPE (NODE_LST_NDE_TYPE), POINTER :: TEMP_PTR TYPE (NODE_LST_NDE_TYPE), POINTER :: TRAIL_PTR ! Find location to insert new node TRAIL_PTR => NODE_LST % HDR TEMP_PTR => TRAIL_PTR % NEXT_PTR DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (NODE_ID < TEMP_PTR % ID) EXIT TRAIL_PTR => TEMP_PTR TEMP_PTR => TEMP_PTR % NEXT_PTR END DO ! Insert new node ALLOCATE (TRAIL_PTR % NEXT_PTR) TRAIL_PTR % NEXT_PTR % ID = NODE_ID TRAIL_PTR % NEXT_PTR % X = X TRAIL_PTR % NEXT_PTR % Y = Y TRAIL_PTR % NEXT_PTR % DOF_INFO = DOF_INFO TRAIL_PTR % NEXT_PTR % CNSTRNT_INDICATOR = CNSTRNT_INDICATOR IF (.NOT. ASSOCIATED (TEMP_PTR)) THEN NULLIFY (TRAIL_PTR % NEXT_PTR % NEXT_PTR) TRAIL_PTR % NEXT_PTR % PREV_PTR => TRAIL_PTR ELSE TRAIL_PTR % NEXT_PTR % NEXT_PTR => TEMP_PTR TRAIL_PTR % NEXT_PTR % PREV_PTR => TRAIL_PTR TEMP_PTR % PREV_PTR => TRAIL_PTR % NEXT_PTR END IF END SUBROUTINE INSERT_NODE SUBROUTINE INSERT_SPR_EL_DATA (SPR_DATA_LST,EL_ID,MAX_NDS_R,MAX_NDS_S) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: EL_ID INTEGER, INTENT (IN) :: MAX_NDS_R INTEGER, INTENT (IN) :: MAX_NDS_S TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: TRAIL_PTR TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: TEMP_PTR ! Find location to insert new node TRAIL_PTR => SPR_DATA_LST % HDR TEMP_PTR => TRAIL_PTR % NEXT DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (EL_ID < TEMP_PTR % EL_ID) EXIT TRAIL_PTR => TEMP_PTR TEMP_PTR => TEMP_PTR % NEXT END DO ! Insert new node ALLOCATE (TRAIL_PTR % NEXT) TRAIL_PTR % NEXT % EL_ID = EL_ID TRAIL_PTR % NEXT % MAX_NDS_R = MAX_NDS_R TRAIL_PTR % NEXT % MAX_NDS_S = MAX_NDS_S TRAIL_PTR % NEXT % NO_SP = 0 NULLIFY (TRAIL_PTR % NEXT % SD_HDR % NEXT) TRAIL_PTR % NEXT % NEXT => TEMP_PTR END SUBROUTINE INSERT_SPR_EL_DATA SUBROUTINE INSERT_SPR_SMPLNG_DATA (SPR_DATA_LST,EL_ID,X,Y,DERIVS) USE CONTROL_MODULE USE DBL_PRECISION_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: EL_ID REAL (DP), INTENT (IN) :: X REAL (DP), INTENT (IN) :: Y REAL (DP), DIMENSION (NO_ROWS_B), INTENT (IN) :: DERIVS TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: EL_DATA_TEMP_PTR TYPE (SMPLNG_DATA_LST_NDE_TYPE), POINTER :: SMPLNG_DATA_TEMP_PTR ! Find correct element slot EL_DATA_TEMP_PTR => SPR_DATA_LST % HDR DO IF (EL_DATA_TEMP_PTR % EL_ID == EL_ID) EXIT EL_DATA_TEMP_PTR => EL_DATA_TEMP_PTR % NEXT END DO ! Add sampling data to the front of the sampling data list SMPLNG_DATA_TEMP_PTR => EL_DATA_TEMP_PTR % SD_HDR % NEXT ALLOCATE (EL_DATA_TEMP_PTR % SD_HDR % NEXT) EL_DATA_TEMP_PTR % SD_HDR % NEXT % X = X EL_DATA_TEMP_PTR % SD_HDR % NEXT % Y = Y ! ALLOCATE (EL_DATA_TEMP_PTR % SD_HDR % NEXT % DERIVS(NO_ROWS_B)) EL_DATA_TEMP_PTR % SD_HDR % NEXT % DERIVS = DERIVS EL_DATA_TEMP_PTR % SD_HDR % NEXT % NEXT => SMPLNG_DATA_TEMP_PTR EL_DATA_TEMP_PTR % NO_SP = EL_DATA_TEMP_PTR % NO_SP + 1 END SUBROUTINE INSERT_SPR_SMPLNG_DATA SUBROUTINE MODIFY_EDGE (EDG_PTR,PNEW,NODE_LST) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDG_PTR INTEGER, INTENT (IN) :: PNEW TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST REAL (DP) :: X1 REAL (DP) :: Y1 REAL (DP) :: X2 REAL (DP) :: Y2 REAL (DP) :: X3 REAL (DP) :: Y3 REAL (DP) :: DELTA REAL (DP) :: RVAL REAL (DP) :: EPS = 0.001 LOGICAL :: ADD TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR REAL (DP) :: NEWX REAL (DP) :: NEWY REAL (DP) :: CURX REAL (DP) :: CURY TYPE (NODE_LST_NDE_TYPE), POINTER :: GET_NODE_PTR TYPE (NODE_LST_NDE_TYPE), POINTER :: NODE_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: REND TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: LEND INTEGER :: INITIAL_CNT REND => EDG_PTR % NS_HDR % RPTR LEND => EDG_PTR % NS_HDR % LPTR X1 = REND % NPTR % X Y1 = REND % NPTR % Y X3 = LEND % NPTR % X Y3 = LEND % NPTR % Y IF (ASSOCIATED (LEND % LPTR,REND)) THEN INITIAL_CNT = 2 X2 = (X1 + X3)/2.0 Y2 = (Y1 + Y3)/2.0 ELSE INITIAL_CNT = 3 X2 = LEND % LPTR % NPTR % X Y2 = LEND % LPTR % NPTR % Y END IF DELTA = 1.d0/PNEW RVAL = 0.d0 DO IF (RVAL > 1.d0+EPS) EXIT NEWX = (2.0*RVAL**2-3.0*RVAL+1.0)*X1 + (4.0*RVAL-4.0*RVAL**2)*X2 + & (2.0*RVAL**2-RVAL)*X3 NEWY = (2.0*RVAL**2-3.0*RVAL+1.0)*Y1 + (4.0*RVAL-4.0*RVAL**2)*Y2 + & (2.0*RVAL**2-RVAL)*Y3 ADD = .TRUE. TEMP_PTR => EDG_PTR % NS_HDR % RPTR DO IF (ASSOCIATED (TEMP_PTR,EDG_PTR % NS_HDR)) EXIT CURX = TEMP_PTR % NPTR % X CURY = TEMP_PTR % NPTR % Y IF (ABS(NEWX-CURX) < EPS .AND. ABS(NEWY-CURY) < EPS) THEN TEMP_PTR % NPTR % DOF_INFO = -1 ADD = .FALSE. EXIT END IF TEMP_PTR => TEMP_PTR % RPTR END DO IF (ADD) THEN CALL INSERT_NODE (NODE_LST,NEXT_UNUSED_NODE,NEWX,NEWY,1,0) NODE_PTR => GET_NODE_PTR (NODE_LST,NEXT_UNUSED_NODE) IF (INITIAL_CNT == 2) THEN CALL NODE_SUBLST_INSERT (NODE_PTR,EDG_PTR % NS_HDR % LPTR % LPTR) ELSE IF (RVAL < 0.5) THEN CALL NODE_SUBLST_INSERT (NODE_PTR, & EDG_PTR % NS_HDR % LPTR % LPTR % LPTR) ELSE CALL NODE_SUBLST_INSERT (NODE_PTR, & EDG_PTR % NS_HDR % LPTR % LPTR) END IF END IF NEXT_UNUSED_NODE = NEXT_UNUSED_NODE + 1 END IF RVAL = RVAL + DELTA END DO END SUBROUTINE MODIFY_EDGE SUBROUTINE MULTIPLY_MATRICES (L, M, N, A, B, C) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: L INTEGER, INTENT (IN) :: M INTEGER, INTENT (IN) :: N REAL (DP), DIMENSION (L,M), INTENT (IN) :: A REAL (DP), DIMENSION (M,N), INTENT (IN) :: B REAL (DP), DIMENSION (L,N), INTENT (OUT) :: C INTEGER :: I INTEGER :: J INTEGER :: K REAL (DP) :: SUM DO I = 1,L DO J = 1,N SUM = 0.d0 DO K = 1,M SUM = SUM + B(K,J) * A(I,K) END DO C(I,J) = SUM END DO END DO END SUBROUTINE MULTIPLY_MATRICES SUBROUTINE NODE_SUBLST_INSERT (NODE_PTR, LOCATION_PTR) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_NDE_TYPE), TARGET :: NODE_PTR TYPE (NODE_SUBLST_NDE_TYPE), TARGET :: LOCATION_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR ! Insert new node ALLOCATE (TEMP_PTR) TEMP_PTR % LPTR => LOCATION_PTR TEMP_PTR % NPTR => NODE_PTR TEMP_PTR % RPTR => LOCATION_PTR % RPTR LOCATION_PTR % RPTR % LPTR => TEMP_PTR LOCATION_PTR % RPTR => TEMP_PTR END SUBROUTINE NODE_SUBLST_INSERT SUBROUTINE OUTPUT_BY_ELEMENTS (ELEMENT_LST, NO_SYS_DOFS, SYS_SOL_VCTR) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST INTEGER, INTENT (IN) :: NO_SYS_DOFS REAL (DP), DIMENSION (NO_SYS_DOFS), INTENT (IN) :: SYS_SOL_VCTR TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ELEMENT_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: HEAD_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR LOGICAL :: RIGHT INTEGER :: BASE INTEGER :: PARAM WRITE (6,"(//,A,//)") & " E L E M E N T A L R E S U L T S" WRITE (6,"(A,I2,A)") "ELEMENT, NODE, XY COORDINATES, ",NO_PARAMS_PER_NODE, & " PARAMETER(S)" ELEMENT_PTR => ELEMENT_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (ELEMENT_PTR)) EXIT HEAD_PTR => ELEMENT_PTR % EDG_1 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(1) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN BASE = NO_PARAMS_PER_NODE*(ABS(TEMP_PTR % NPTR % DOF_INFO)-1) WRITE (6,"(2I5,(9(1X,1PE12.5)))") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & (SYS_SOL_VCTR(BASE+PARAM), PARAM=1,NO_PARAMS_PER_NODE) ELSE WRITE (6,"(2I5,(2(1X,1PE12.5)),A)") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & " ------------------" END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_2 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(2) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN BASE = NO_PARAMS_PER_NODE*(ABS(TEMP_PTR % NPTR % DOF_INFO)-1) WRITE (6,"(2I5,(9(1X,1PE12.5)))") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & (SYS_SOL_VCTR(BASE+PARAM), PARAM=1,NO_PARAMS_PER_NODE) ELSE WRITE (6,"(2I5,(2(1X,1PE12.5)),A)") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & " ------------------" END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_3 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(3) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (ASSOCIATED (TEMP_PTR,HEAD_PTR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN BASE = NO_PARAMS_PER_NODE*(ABS(TEMP_PTR % NPTR % DOF_INFO)-1) WRITE (6,"(2I5,(9(1X,1PE12.5)))") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & (SYS_SOL_VCTR(BASE+PARAM), PARAM=1,NO_PARAMS_PER_NODE) ELSE WRITE (6,"(2I5,(2(1X,1PE12.5)),A)") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & " ------------------" END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO HEAD_PTR => ELEMENT_PTR % EDG_4 % NS_HDR IF (ELEMENT_PTR % CRNR_SEQ(4) == HEAD_PTR % RPTR % NPTR % ID) THEN RIGHT = .TRUE. TEMP_PTR => HEAD_PTR % RPTR % RPTR ELSE RIGHT = .FALSE. TEMP_PTR => HEAD_PTR % LPTR % LPTR END IF DO IF (TEMP_PTR % NPTR % ID == ELEMENT_PTR % CRNR_SEQ(1)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN BASE = NO_PARAMS_PER_NODE*(ABS(TEMP_PTR % NPTR % DOF_INFO)-1) WRITE (6,"(2I5,(9(1X,1PE12.5)))") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & (SYS_SOL_VCTR(BASE+PARAM), PARAM=1,NO_PARAMS_PER_NODE) ELSE WRITE (6,"(2I5,(2(1X,1PE12.5)),A)") ELEMENT_PTR % ID, & TEMP_PTR % NPTR % ID, TEMP_PTR % NPTR % X, TEMP_PTR % NPTR % Y, & " ------------------" END IF IF (RIGHT) THEN TEMP_PTR => TEMP_PTR % RPTR ELSE TEMP_PTR => TEMP_PTR % LPTR END IF END DO ELEMENT_PTR => ELEMENT_PTR % NEXT END DO END SUBROUTINE OUTPUT_BY_ELEMENTS SUBROUTINE OUTPUT_BY_NODES (NODE_LST, NO_SYS_DOFS, SYS_SOL_VCTR) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER, INTENT (IN) :: NO_SYS_DOFS REAL (DP), DIMENSION (NO_SYS_DOFS), INTENT (IN) :: SYS_SOL_VCTR TYPE (NODE_LST_NDE_TYPE), POINTER :: NODE_PTR INTEGER :: BASE INTEGER :: PARAM WRITE (6,"(//,A,//)") & " N O D A L R E S U L T S" SELECT CASE ( NO_PARAMS_PER_NODE ) CASE ( 1 ) WRITE (6,"(2A,/)") & " NODE ID X-COORDINATE Y-COORDINATE", & " PARAMETER" CASE ( 2 ) WRITE (6,"(/)") CASE DEFAULT WRITE (6,"(/)") END SELECT NODE_PTR => NODE_LST % HDR % NEXT_PTR DO IF (.NOT. ASSOCIATED (NODE_PTR)) EXIT IF (NODE_PTR % DOF_INFO /= 0) THEN BASE = NO_PARAMS_PER_NODE*(ABS(NODE_PTR % DOF_INFO)-1) SELECT CASE ( NO_PARAMS_PER_NODE ) CASE ( 1 ) WRITE (6,"(I15,9X,E12.5,6X,E12.5,4X,E12.5)") & NODE_PTR % ID, NODE_PTR % X, NODE_PTR % Y, & (SYS_SOL_VCTR(BASE+PARAM), PARAM=1,NO_PARAMS_PER_NODE) CASE ( 2 ) WRITE (6,"(/)") CASE DEFAULT WRITE (6,"(/)") END SELECT ELSE WRITE (6,"(I5,(2(1X,1PE12.5)),A)") NODE_PTR % ID, NODE_PTR % X, & NODE_PTR % Y," ------------------" END IF NODE_PTR => NODE_PTR % NEXT_PTR END DO END SUBROUTINE OUTPUT_BY_NODES SUBROUTINE OUTPUT_NDL_GRADS (NO_AN_NODES,NO_CMPNTS,NODE_LST,GRADS) USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_AN_NODES INTEGER, INTENT (IN) :: NO_CMPNTS TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST REAL (DP), DIMENSION (NO_AN_NODES,NO_CMPNTS), INTENT (IN) :: GRADS TYPE (NODE_LST_NDE_TYPE), POINTER :: NODE_PTR INTEGER :: INDX WRITE (6,"(//,A,//)") & " G L O B A L N O D E P O I N T D E R I V A T I V E S" WRITE (6,"(A,/)") & " NODE ID X-COORDINATE Y-COORDINATE DP/DX DP/DY" NODE_PTR => NODE_LST % HDR % NEXT_PTR DO IF (.NOT. ASSOCIATED (NODE_PTR)) EXIT INDX = ABS(NODE_PTR % DOF_INFO) WRITE (6,"(I9,3X,2(2X,1PE12.4),2(2X,1PE15.5))") & NODE_PTR % ID, NODE_PTR % X, NODE_PTR % Y, GRADS(INDX,1), GRADS(INDX,2) NODE_PTR => NODE_PTR % NEXT_PTR END DO END SUBROUTINE OUTPUT_NDL_GRADS SUBROUTINE P_ADAPT (EL_LST,EDG_LST,NODE_LST,EL_ENORM,CUTOFF) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDG_LST TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST REAL (DP), DIMENSION (NO_ELEMENTS), INTENT (IN) :: EL_ENORM REAL (DP), INTENT (IN) :: CUTOFF INTEGER, DIMENSION (NO_ELEMENTS*4) :: POLD INTEGER, DIMENSION (NO_ELEMENTS*4) :: PNEW REAL (DP), DIMENSION (NO_ELEMENTS*4) :: PDBL TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR REAL (DP) :: PTMP REAL (DP), DIMENSION (NO_ELEMENTS) :: XI INTEGER :: NO_MODS INTEGER :: MDFCTN INTEGER, DIMENSION (1) :: EL_TO_MDFY INTEGER :: EDG_INDX TYPE (EDGE_LST_NDE_TYPE), POINTER :: EDG_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR LOGICAL :: OK_TO_MODIFY INTEGER :: GET_NO_NODES_ON_EDGE EL_PTR => EL_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_PTR)) EXIT XI(EL_PTR % ID) = EL_ENORM(EL_PTR % ID)/CUTOFF print*,"xi(",el_ptr % id,") = ",xi(el_ptr % id) EL_PTR => EL_PTR % NEXT END DO POLD = 0 PNEW = 0 PDBL = 0.d0 NO_MODS = CEILING (NO_ELEMENTS*0.3) print*,"no_mods = ",no_mods DO MDFCTN = 1,NO_MODS EL_TO_MDFY = MAXLOC (XI) print*,"el_to_mdfy = ",el_to_mdfy EL_PTR => EL_LST % HDR % NEXT DO IF (EL_PTR % ID == EL_TO_MDFY(1)) EXIT EL_PTR => EL_PTR % NEXT END DO POLD(EL_PTR % EDG_1 % ID) = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_1) - 1 PTMP = POLD(EL_PTR % EDG_1 % ID) * & XI(EL_PTR % ID)**(1.0/POLD(EL_PTR % EDG_1 % ID)) IF (PTMP > PDBL(EL_PTR % EDG_1 % ID)) PDBL(EL_PTR % EDG_1 % ID) = PTMP POLD(EL_PTR % EDG_2 % ID) = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_2) - 1 PTMP = POLD(EL_PTR % EDG_2 % ID) * & XI(EL_PTR % ID)**(1.0/POLD(EL_PTR % EDG_2 % ID)) IF (PTMP > PDBL(EL_PTR % EDG_2 % ID)) PDBL(EL_PTR % EDG_2 % ID) = PTMP POLD(EL_PTR % EDG_3 % ID) = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_3) - 1 PTMP = POLD(EL_PTR % EDG_3 % ID) * & XI(EL_PTR % ID)**(1.0/POLD(EL_PTR % EDG_3 % ID)) IF (PTMP > PDBL(EL_PTR % EDG_3 % ID)) PDBL(EL_PTR % EDG_3 % ID) = PTMP POLD(EL_PTR % EDG_4 % ID) = GET_NO_NODES_ON_EDGE (EL_PTR % EDG_4) - 1 PTMP = POLD(EL_PTR % EDG_4 % ID) * & XI(EL_PTR % ID)**(1.0/POLD(EL_PTR % EDG_4 % ID)) IF (PTMP > PDBL(EL_PTR % EDG_4 % ID)) PDBL(EL_PTR % EDG_4 % ID) = PTMP XI(EL_TO_MDFY(1)) = 0.d0 END DO DO EDG_INDX = 1,NO_ELEMENTS*4 IF (PDBL(EDG_INDX) /= 0.d0) THEN IF (PDBL(EDG_INDX) < 1.d0) THEN PNEW(EDG_INDX) = 1 ELSE IF (PDBL(EDG_INDX) > POLD(EDG_INDX)) THEN PNEW(EDG_INDX) = CEILING (PDBL(EDG_INDX)) IF ((PNEW(EDG_INDX)-POLD(EDG_INDX)) > 2) & PNEW(EDG_INDX) = POLD(EDG_INDX) + 2 ELSE PNEW(EDG_INDX) = FLOOR (PDBL(EDG_INDX)) END IF END IF END DO EDG_PTR => EDG_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EDG_PTR)) EXIT IF (PNEW(EDG_PTR % ID) /= 0) THEN IF (PNEW(EDG_PTR % ID) /= POLD(EDG_PTR % ID)) THEN OK_TO_MODIFY = .FALSE. TEMP_PTR => EDG_PTR % NS_HDR % RPTR DO IF (ASSOCIATED (TEMP_PTR,EDG_PTR % NS_HDR)) EXIT IF (TEMP_PTR % NPTR % CNSTRNT_INDICATOR == 0) THEN OK_TO_MODIFY = .TRUE. EXIT END IF TEMP_PTR => TEMP_PTR % RPTR END DO IF (OK_TO_MODIFY) THEN print*,"modifying edge # ",EDG_PTR % ID print*,"pnew = ",PNEW(EDG_PTR % ID) CALL REINIT_EDGE (EDG_PTR) CALL MODIFY_EDGE (EDG_PTR,PNEW(EDG_PTR % ID),NODE_LST) END IF END IF END IF EDG_PTR => EDG_PTR % NEXT END DO END SUBROUTINE P_ADAPT SUBROUTINE PERFORM_EL_POST_PRCSSNG (ELEMENT_PTR,NO_AN_NODES,NO_EL_DOFS, & EL_SOL_VCTR,SPR_DATA_LST,BANNER) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_NDE_TYPE), TARGET :: ELEMENT_PTR INTEGER, INTENT (IN) :: NO_AN_NODES INTEGER, INTENT (IN) :: NO_EL_DOFS REAL (DP), DIMENSION (NO_EL_DOFS), INTENT (IN) :: EL_SOL_VCTR TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST LOGICAL, INTENT (INOUT) :: BANNER INTEGER :: NIP INTEGER :: QP REAL (DP), DIMENSION (SOL_SPACE_DIM) :: QP_XY REAL (DP), DIMENSION (SOL_SPACE_DIM,NO_AN_NODES) :: DGH REAL (DP), DIMENSION (SOL_SPACE_DIM) :: GRAD INTEGER :: I IF (BANNER) THEN WRITE (6,"(//,A,//)") & " G L O B A L G A U S S P O I N T D E R I V A T I V E S" WRITE (6,"(A,/)") & " ELEMENT ID IP X-COORDINATE Y-COORDINATE DP/DX DP/DY" BANNER = .FALSE. END IF READ (POST_PROC_IO_UNIT1) NIP DO QP = 1,NIP READ (POST_PROC_IO_UNIT1) QP_XY, DGH GRAD(1) = 0.d0 GRAD(2) = 0.d0 DO I = 1,NO_AN_NODES GRAD(1) = GRAD(1) + DGH(1,I)*EL_SOL_VCTR(I) GRAD(2) = GRAD(2) + DGH(2,I)*EL_SOL_VCTR(I) END DO CALL INSERT_SPR_SMPLNG_DATA (SPR_DATA_LST,ELEMENT_PTR % ID,QP_XY(1), & QP_XY(2),GRAD) WRITE (6,"(I8,I8,2(2X,1PE12.4),2(2X,1PE15.5))") & ELEMENT_PTR % ID, QP, QP_XY, GRAD END DO END SUBROUTINE PERFORM_EL_POST_PRCSSNG SUBROUTINE PERFORM_POST_PRCSSNG (ELEMENT_LST,NO_SYS_DOFS,SYS_SOL_VCTR, & SPR_DATA_LST) USE DBL_PRECISION_MODULE USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST INTEGER, INTENT (IN) :: NO_SYS_DOFS REAL (DP), DIMENSION (NO_SYS_DOFS), INTENT (IN) :: SYS_SOL_VCTR TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER :: NO_AN_NODES INTEGER :: NO_EL_DOFS INTEGER, DIMENSION (:), ALLOCATABLE :: DOF_MAP REAL (DP), DIMENSION (:,:), ALLOCATABLE :: NODE_CRDS REAL (DP), DIMENSION (:), ALLOCATABLE :: EL_SOL_VCTR TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ELEMENT_PTR INTEGER :: GET_NO_ANALYSIS_NODES INTEGER :: INDEX LOGICAL :: BANNER IF (POST_PROC_IO_UNIT1 > 0) THEN BANNER = .TRUE. REWIND POST_PROC_IO_UNIT1 ELEMENT_PTR => ELEMENT_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (ELEMENT_PTR)) EXIT NO_AN_NODES = GET_NO_ANALYSIS_NODES (ELEMENT_PTR) NO_EL_DOFS = NO_AN_NODES * NO_PARAMS_PER_NODE ALLOCATE (DOF_MAP(NO_EL_DOFS)) ALLOCATE (NODE_CRDS(NO_AN_NODES,SOL_SPACE_DIM)) ALLOCATE (EL_SOL_VCTR(NO_EL_DOFS)) CALL GET_SYS_DOFS_FOR_ELEMENT (ELEMENT_PTR,NO_EL_DOFS,DOF_MAP) CALL GET_NODAL_COORDINATES (ELEMENT_PTR,NO_AN_NODES,NODE_CRDS) DO INDEX = 1,NO_EL_DOFS EL_SOL_VCTR(INDEX) = SYS_SOL_VCTR(DOF_MAP(INDEX)) END DO CALL PERFORM_EL_POST_PRCSSNG (ELEMENT_PTR,NO_AN_NODES,NO_EL_DOFS, & EL_SOL_VCTR,SPR_DATA_LST,BANNER) DEALLOCATE (DOF_MAP) DEALLOCATE (NODE_CRDS) DEALLOCATE (EL_SOL_VCTR) ELEMENT_PTR => ELEMENT_PTR % NEXT END DO END IF END SUBROUTINE PERFORM_POST_PRCSSNG SUBROUTINE PERFORM_SKYLINE_ANALYSIS (NO_SYS_DOFS, ELEMENT_LST, & SYS_COL_HGHTS, DIAG_LCTNS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_DOFS TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST INTEGER, DIMENSION (NO_SYS_DOFS), INTENT (OUT) :: SYS_COL_HGHTS INTEGER, DIMENSION (NO_SYS_DOFS+1), INTENT (OUT) :: DIAG_LCTNS TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: ELEMENT_PTR INTEGER, DIMENSION (:), ALLOCATABLE :: DOF_MAP INTEGER :: NO_AN_NODES INTEGER :: NO_EL_DOFS INTEGER :: GET_NO_ANALYSIS_NODES INTEGER, DIMENSION (:), ALLOCATABLE :: EL_COL_HGHTS INTEGER :: EL_DOF INTEGER :: SYS_DOF SYS_COL_HGHTS = 0 DIAG_LCTNS = 0 ELEMENT_PTR => ELEMENT_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (ELEMENT_PTR)) EXIT NO_AN_NODES = GET_NO_ANALYSIS_NODES (ELEMENT_PTR) NO_EL_DOFS = NO_AN_NODES * NO_PARAMS_PER_NODE ALLOCATE (DOF_MAP(NO_EL_DOFS)) CALL GET_SYS_DOFS_FOR_ELEMENT (ELEMENT_PTR, NO_EL_DOFS, DOF_MAP) ALLOCATE (EL_COL_HGHTS(NO_EL_DOFS)) CALL GET_ELEMENT_COLUMN_HEIGHTS (NO_EL_DOFS, DOF_MAP, EL_COL_HGHTS) DO EL_DOF = 1,NO_EL_DOFS SYS_DOF = DOF_MAP(EL_DOF) SYS_COL_HGHTS(SYS_DOF) = MAX (SYS_COL_HGHTS(SYS_DOF), & EL_COL_HGHTS(EL_DOF)) END DO ELEMENT_PTR => ELEMENT_PTR % NEXT DEALLOCATE (DOF_MAP) DEALLOCATE (EL_COL_HGHTS) END DO CALL GET_SYS_DIAG_LCTNS (NO_SYS_DOFS, SYS_COL_HGHTS, DIAG_LCTNS) END SUBROUTINE PERFORM_SKYLINE_ANALYSIS SUBROUTINE PREPARE_COMPLETE_LS_PRBLM (M,N,DGR,SPR_DATA_LST,EL_LST, & MAX_PTCH_ELS,PTCH_IDS,P,SGMA_H) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: M INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: DGR TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS REAL (DP), DIMENSION (M,N), INTENT (OUT) :: P REAL (DP), DIMENSION (M,NO_ROWS_B), INTENT (OUT) :: SGMA_H INTEGER :: ROW INTEGER :: INDEX INTEGER :: CMPNT TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: EL_DATA_PTR TYPE (SMPLNG_DATA_LST_NDE_TYPE), POINTER :: SMPLNG_DATA_PTR REAL (DP) :: X REAL (DP) :: Y REAL (DP) :: XMIN REAL (DP) :: XMAX REAL (DP) :: YMIN REAL (DP) :: YMAX INTEGER :: I INTEGER :: J INTEGER :: TERM REAL (DP) :: XPRT REAL (DP) :: YPRT CALL DETERMINE_PTCH_BNDS (MAX_PTCH_ELS,PTCH_IDS,EL_LST,XMIN,XMAX,YMIN,YMAX) ROW = 0 DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN EL_DATA_PTR => SPR_DATA_LST % HDR % NEXT DO IF (PTCH_IDS(INDEX) == EL_DATA_PTR % EL_ID) EXIT EL_DATA_PTR => EL_DATA_PTR % NEXT END DO SMPLNG_DATA_PTR => EL_DATA_PTR % SD_HDR % NEXT DO IF (.NOT. ASSOCIATED (SMPLNG_DATA_PTR)) EXIT ROW = ROW + 1 DO CMPNT = 1,NO_ROWS_B SGMA_H(ROW,CMPNT) = SMPLNG_DATA_PTR % DERIVS(CMPNT) END DO X = -1.d0+2.d0*((SMPLNG_DATA_PTR % X - XMIN)/(XMAX-XMIN)) Y = -1.d0+2.d0*((SMPLNG_DATA_PTR % Y - YMIN)/(YMAX-YMIN)) TERM = 0 DO I = 0,DGR DO J = 0,I TERM = TERM + 1 IF ((I-J) == 0) THEN XPRT = 1.d0 ELSE IF (X == 0.d0) THEN XPRT = 0.d0 ELSE XPRT = X**(I-J) END IF IF ((J) == 0) THEN YPRT = 1.d0 ELSE IF (Y == 0.d0) THEN YPRT = 0.d0 ELSE YPRT = Y**J END IF P(ROW,TERM) = XPRT * YPRT END DO END DO SMPLNG_DATA_PTR => SMPLNG_DATA_PTR % NEXT END DO END IF END DO END SUBROUTINE PREPARE_COMPLETE_LS_PRBLM SUBROUTINE PREPARE_SHPFNCTN_LS_PRBLM (M,N,MAX_NODES_R,MAX_NODES_S, & EL_LST,SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS,P,SGMA_H) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: M INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: MAX_NODES_R INTEGER, INTENT (IN) :: MAX_NODES_S TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST INTEGER, INTENT (IN) :: MAX_PTCH_ELS INTEGER, DIMENSION (MAX_PTCH_ELS), INTENT (IN) :: PTCH_IDS REAL (DP), DIMENSION (M,N), INTENT (OUT) :: P REAL (DP), DIMENSION (M,NO_ROWS_B), INTENT (OUT) :: SGMA_H INTEGER :: ROW INTEGER :: INDEX INTEGER :: CMPNT TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: EL_DATA_PTR TYPE (SMPLNG_DATA_LST_NDE_TYPE), POINTER :: SMPLNG_DATA_PTR REAL (DP) :: X REAL (DP) :: Y REAL (DP) :: XMIN REAL (DP) :: XMAX REAL (DP) :: YMIN REAL (DP) :: YMAX INTEGER :: I INTEGER :: TERM CALL DETERMINE_PTCH_BNDS (MAX_PTCH_ELS,PTCH_IDS,EL_LST,XMIN,XMAX,YMIN,YMAX) ROW = 0 DO INDEX = 1,MAX_PTCH_ELS IF (PTCH_IDS(INDEX) /= 0) THEN EL_DATA_PTR => SPR_DATA_LST % HDR % NEXT DO IF (PTCH_IDS(INDEX) == EL_DATA_PTR % EL_ID) EXIT EL_DATA_PTR => EL_DATA_PTR % NEXT END DO SMPLNG_DATA_PTR => EL_DATA_PTR % SD_HDR % NEXT DO IF (.NOT. ASSOCIATED (SMPLNG_DATA_PTR)) EXIT ROW = ROW + 1 DO CMPNT = 1,NO_ROWS_B SGMA_H(ROW,CMPNT) = SMPLNG_DATA_PTR % DERIVS(CMPNT) END DO X = -1.d0 + 2.d0*((SMPLNG_DATA_PTR % X - XMIN)/(XMAX - XMIN)) Y = -1.d0 + 2.d0*((SMPLNG_DATA_PTR % Y - YMIN)/(YMAX - YMIN)) TERM = 0 DO I = 0,MAX_NODES_R-1 TERM = TERM + 1 IF (I == 0) THEN P(ROW,TERM) = 1.d0 ELSE IF (X == 0.d0) THEN P(ROW,TERM) = 0.d0 ELSE P(ROW,TERM) = X**I END IF END DO DO I = 0,MAX_NODES_R-1 TERM = TERM + 1 IF (I == 0) THEN P(ROW,TERM) = Y ELSE IF (X == 0.0 .OR. Y == 0) THEN P(ROW,TERM) = 0.d0 ELSE P(ROW,TERM) = X**I * Y END IF END DO DO I = 3,MAX_NODES_S TERM = TERM + 1 IF (Y == 0.d0) THEN P(ROW,TERM) = 0.d0 ELSE P(ROW,TERM) = Y**(I-1) END IF END DO DO I = 3,MAX_NODES_S TERM = TERM + 1 IF (Y == 0.0 .OR. X == 0) THEN P(ROW,TERM) = 0.d0 ELSE P(ROW,TERM) = X * Y**(I-1) END IF END DO SMPLNG_DATA_PTR => SMPLNG_DATA_PTR % NEXT END DO END IF END DO END SUBROUTINE PREPARE_SHPFNCTN_LS_PRBLM SUBROUTINE PRINT_EDGE_LIST (EDGE_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST TYPE (EDGE_LST_NDE_TYPE), POINTER :: TEMP_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: NS_TEMP_PTR TEMP_PTR => EDGE_LST % HDR % NEXT WRITE (6,"(/,A23,/)") "######## EDGE LIST DUMP" WRITE (6,"(A12)") "EDGE_LST.HDR" DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT WRITE (6,"(A1)") "|" WRITE (6,"(A1)") "|" WRITE (6,"(A1)") "V" WRITE (6,"(A5,I3)") "ID = ",TEMP_PTR % ID WRITE (6,"(A20)") "NS = NODE_SUBLST.HDR" NS_TEMP_PTR => TEMP_PTR % NS_HDR % RPTR DO IF (ASSOCIATED (NS_TEMP_PTR, TEMP_PTR % NS_HDR)) EXIT WRITE (6,"(A6)") " |" WRITE (6,"(A6)") " |" WRITE (6,"(A6)") " V" WRITE (6,"(A10,I3)") " ID = ",NS_TEMP_PTR % NPTR % ID NS_TEMP_PTR => NS_TEMP_PTR % RPTR END DO TEMP_PTR => TEMP_PTR % NEXT END DO WRITE (6,"(/,A27,/)") "######## END EDGE LIST DUMP" END SUBROUTINE PRINT_EDGE_LIST SUBROUTINE PRINT_ELEMENT_LIST (ELEMENT_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: TEMP_PTR TEMP_PTR => ELEMENT_LST % HDR % NEXT WRITE (6,"(/,A26,/)") "######## ELEMENT LIST DUMP" WRITE (6,"(A15)") "ELEMENT_LST.HDR" DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT WRITE (6,"(A1)") "|" WRITE (6,"(A1)") "|" WRITE (6,"(A1)") "V" WRITE (6,"(A5,I3)") "ID = ",TEMP_PTR % ID WRITE (6,"(A5,I3)") "C1 = ",TEMP_PTR % CRNR_SEQ (1) WRITE (6,"(A5,I3)") "C2 = ",TEMP_PTR % CRNR_SEQ (2) WRITE (6,"(A5,I3)") "C3 = ",TEMP_PTR % CRNR_SEQ (3) WRITE (6,"(A5,I3)") "C4 = ",TEMP_PTR % CRNR_SEQ (4) WRITE (6,"(A5,I3)") "E1 = ",TEMP_PTR % EDG_1 % ID WRITE (6,"(A5,I3)") "E2 = ",TEMP_PTR % EDG_2 % ID WRITE (6,"(A5,I3)") "E3 = ",TEMP_PTR % EDG_3 % ID WRITE (6,"(A5,I3)") "E4 = ",TEMP_PTR % EDG_4 % ID IF (.NOT.ASSOCIATED (TEMP_PTR % ADJ_1)) THEN WRITE (6,"(A)") "A1 = 0" ELSE WRITE (6,"(A5,I3)") "A1 = ",TEMP_PTR % ADJ_1 % ID END IF IF (.NOT.ASSOCIATED (TEMP_PTR % ADJ_2)) THEN WRITE (6,"(A)") "A2 = 0" ELSE WRITE (6,"(A5,I3)") "A2 = ",TEMP_PTR % ADJ_2 % ID END IF IF (.NOT.ASSOCIATED (TEMP_PTR % ADJ_3)) THEN WRITE (6,"(A)") "A3 = 0" ELSE WRITE (6,"(A5,I3)") "A3 = ",TEMP_PTR % ADJ_3 % ID END IF IF (.NOT.ASSOCIATED (TEMP_PTR % ADJ_4)) THEN WRITE (6,"(A)") "A4 = 0" ELSE WRITE (6,"(A5,I3)") "A4 = ",TEMP_PTR % ADJ_4 % ID END IF TEMP_PTR => TEMP_PTR % NEXT END DO WRITE (6,"(/,A30,/)") "######## END ELEMENT LIST DUMP" END SUBROUTINE PRINT_ELEMENT_LIST SUBROUTINE PRINT_LOGO () IMPLICIT NONE INTEGER, DIMENSION (8) :: STAMP CHARACTER (LEN = 3) :: MONTH CALL DATE_AND_TIME (VALUES=STAMP) SELECT CASE (STAMP(2)) CASE (1) MONTH = "JAN" CASE (2) MONTH = "FEB" CASE (3) MONTH = "MAR" CASE (4) MONTH = "APR" CASE (5) MONTH = "MAY" CASE (6) MONTH = "JUN" CASE (7) MONTH = "JUL" CASE (8) MONTH = "AUG" CASE (9) MONTH = "SEP" CASE (10) MONTH = "OCT" CASE (11) MONTH = "NOV" CASE (12) MONTH = "DEC" END SELECT WRITE (6,"(2A)") "****************************************", & "****************************************" WRITE (6,"(2A)") "****************************************", & "****************************************" WRITE (6,"(2A)") "*** ", & " ***" WRITE (6,"(3A)") "*** ", & "RP_ADAPT (Version 1.0)", & " ***" WRITE (6,"(2A)") "*** ", & " ***" WRITE (6,"(3A,I2,A,I4,A)") "*** ", & MONTH," ",STAMP(3),", ",STAMP(1), & " ***" WRITE (6,"(2A)") "*** ", & " ***" WRITE (6,"(A,3(I2,A))") "*** ", & STAMP(5),":",STAMP(6),":",STAMP(7), & " ***" WRITE (6,"(2A)") "*** ", & " ***" WRITE (6,"(2A)") "****************************************", & "****************************************" WRITE (6,"(2A)") "****************************************", & "****************************************" END SUBROUTINE PRINT_LOGO SUBROUTINE PRINT_NODE_LIST (NODE_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST TYPE (NODE_LST_NDE_TYPE), POINTER :: TEMP_PTR ! TYPE (EDGE_SUBLST_NDE_TYPE), POINTER :: ES_TEMP_PTR TEMP_PTR => NODE_LST % HDR % NEXT_PTR WRITE (6,"(/,A23,/)") "######## NODE LIST DUMP" WRITE (6,"(A12)") "NODE_LST.HDR" DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT WRITE (6,"(A)") "|" WRITE (6,"(A)") "|" WRITE (6,"(A)") "V" WRITE (6,"(A5,I3)") "ID = ",TEMP_PTR % ID WRITE (6,"(A5,F6.2)") "X = ",TEMP_PTR % X WRITE (6,"(A5,F6.2)") "Y = ",TEMP_PTR % Y ! WRITE (6,"(A20)") "ES = EDGE_SUBLST.HDR" ! ES_TEMP_PTR => TEMP_PTR % EDGE_SUBLST % HDR % NEXT ! DO ! ! IF (.NOT. ASSOCIATED (ES_TEMP_PTR)) EXIT ! WRITE (6,"(A6)") " |" ! WRITE (6,"(A6)") " |" ! WRITE (6,"(A6)") " V" ! WRITE (6,"(A10,I3)") " ID = ",ES_TEMP_PTR % ID ! ES_TEMP_PTR => ES_TEMP_PTR % NEXT ! ! END DO WRITE (6,"(A5,I3)") "IN = ",TEMP_PTR % DOF_INFO WRITE (6,"(A5,I3)") "CI = ",TEMP_PTR % CNSTRNT_INDICATOR TEMP_PTR => TEMP_PTR % NEXT_PTR END DO WRITE (6,"(/,A27,/)") "######## END NODE LIST DUMP" END SUBROUTINE PRINT_NODE_LIST SUBROUTINE PRINT_SPR_DATA_LIST (SPR_DATA_LST) USE CONTROL_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: EL_DATA_TEMP_PTR TYPE (SMPLNG_DATA_LST_NDE_TYPE), POINTER :: SMPLNG_DATA_TEMP_PTR INTEGER :: I EL_DATA_TEMP_PTR => SPR_DATA_LST % HDR % NEXT WRITE (6,"(/,A,/)") "######## SPR DATA LIST DUMP" WRITE (6,"(A)") "SPR_DATA_LST.HDR" DO IF (.NOT. ASSOCIATED (EL_DATA_TEMP_PTR)) EXIT WRITE (6,"(A)") "|" WRITE (6,"(A)") "|" WRITE (6,"(A)") "V" WRITE (6,"(A,I3)") "EL_ID = ",EL_DATA_TEMP_PTR % EL_ID WRITE (6,"(A,I3)") "MAX_NDS_R = ",EL_DATA_TEMP_PTR % MAX_NDS_R WRITE (6,"(A,I3)") "MAX_NDS_S = ",EL_DATA_TEMP_PTR % MAX_NDS_S WRITE (6,"(A,I3)") "NO_SP = ",EL_DATA_TEMP_PTR % NO_SP WRITE (6,"(A)") "SD_HDR = SMPLNG_DATA.HDR" SMPLNG_DATA_TEMP_PTR => EL_DATA_TEMP_PTR % SD_HDR % NEXT DO IF (.NOT. ASSOCIATED (SMPLNG_DATA_TEMP_PTR)) EXIT WRITE (6,"(A)") " |" WRITE (6,"(A)") " |" WRITE (6,"(A)") " V" WRITE (6,"(A,E12.5)") " X = ",SMPLNG_DATA_TEMP_PTR % X WRITE (6,"(A,E12.5)") " Y = ",SMPLNG_DATA_TEMP_PTR % Y DO I = 1,NO_ROWS_B WRITE (6,"(A,I2,A,E12.5)") & " ",I," = ",SMPLNG_DATA_TEMP_PTR % DERIVS(I) END DO SMPLNG_DATA_TEMP_PTR => SMPLNG_DATA_TEMP_PTR % NEXT END DO EL_DATA_TEMP_PTR => EL_DATA_TEMP_PTR % NEXT END DO WRITE (6,"(/,A,/)") "######## END SPR DATA LIST DUMP" END SUBROUTINE PRINT_SPR_DATA_LIST MODULE PROPERTIES_MODULE USE DBL_PRECISION_MODULE INTEGER, DIMENSION (:,:), ALLOCATABLE :: NODE_PRPS_INT REAL (DP), DIMENSION (:,:), ALLOCATABLE :: NODE_PRPS_REAL INTEGER, DIMENSION (:,:), ALLOCATABLE :: EL_PRPS_INT REAL (DP), DIMENSION (:,:), ALLOCATABLE :: EL_PRPS_REAL INTEGER, DIMENSION (:,:), ALLOCATABLE :: SEG_PRPS_INT REAL (DP), DIMENSION (:,:), ALLOCATABLE :: SEG_PRPS_REAL INTEGER, DIMENSION (:), ALLOCATABLE :: MISC_PRPS_INT REAL (DP), DIMENSION (:), ALLOCATABLE :: MISC_PRPS_REAL END MODULE PROPERTIES_MODULE SUBROUTINE PURGE_SPR_DATA_LIST (SPR_DATA_LST) USE SPR_TYPES_MODULE IMPLICIT NONE TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (INOUT) :: SPR_DATA_LST TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: EL_DATA_TEMP_PTR TYPE (SMPLNG_DATA_LST_NDE_TYPE), POINTER :: SMPLNG_DATA_TEMP_PTR EL_DATA_TEMP_PTR => SPR_DATA_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_DATA_TEMP_PTR)) EXIT SMPLNG_DATA_TEMP_PTR => EL_DATA_TEMP_PTR % SD_HDR % NEXT DO IF (.NOT. ASSOCIATED (SMPLNG_DATA_TEMP_PTR)) EXIT EL_DATA_TEMP_PTR % SD_HDR % NEXT => SMPLNG_DATA_TEMP_PTR % NEXT ! DEALLOCATE (SMPLNG_DATA_TEMP_PTR % DERIVS) DEALLOCATE (SMPLNG_DATA_TEMP_PTR) SMPLNG_DATA_TEMP_PTR => EL_DATA_TEMP_PTR % SD_HDR % NEXT END DO SPR_DATA_LST % HDR % NEXT => EL_DATA_TEMP_PTR % NEXT DEALLOCATE (EL_DATA_TEMP_PTR) EL_DATA_TEMP_PTR => SPR_DATA_LST % HDR % NEXT END DO END SUBROUTINE PURGE_SPR_DATA_LIST SUBROUTINE READ_CONTROL_DATA () USE CONTROL_MODULE IMPLICIT NONE INTEGER :: LINE_CNT ! Loop control CHARACTER (LEN = 80) :: REMARK_LINE ! Dummy used to read user remarks ! Print the logo, run date, and run time CALL PRINT_LOGO () ! Read and echo the title line READ (5,"(A80)") TITLE WRITE (6,"(/,A)") TITLE ! Read the first line of problem control data READ (5,"(16I5)") NO_SYS_GEOM_NODES, & NO_ELEMENTS, & NO_PARAMS_PER_NODE, & UNUSED_1, & SOL_SPACE_DIM, & UNUSED_2, & UNUSED_3, & UNUSED_4, & UNUSED_5, & RHS_FLAG, & NO_USER_REMARKS, & NO_ROWS_B, & NO_QUAD_PTS, & UNUSED_6, & UNUSED_7, & DEBUG_MODE ! Echo this data back to the user WRITE (6,"(//,2A,//)") & " P R O B L E M P A R A M ", & "E T E R E C H O" WRITE (6,"(A,I5)") & "NUMBER OF NODAL POINTS IN SYSTEM.............", & NO_SYS_GEOM_NODES WRITE (6,"(A,I5)") & "NUMBER OF ELEMENTS IN SYSTEM.................", & NO_ELEMENTS WRITE (6,"(A,I5)") & "NUMBER OF PARAMETERS PER NODE (1)............", & NO_PARAMS_PER_NODE WRITE (6,"(A,I5)") & "DIMENSION OF SPACE (1).......................", & SOL_SPACE_DIM WRITE (6,"(A,I5)") & "INITIAL FORCING VECTOR (0-OMIT, 1-READ)......", & RHS_FLAG WRITE (6,"(A,I5)") & "NUMBER OF USER REMARK LINES..................", & NO_USER_REMARKS WRITE (6,"(A,I5)") & "NUMBER OF ROWS IN B MATRIX (1)...............", & NO_ROWS_B WRITE (6,"(A,I5)") & "NUMBER OF QUADRATURE POINTS..................", & NO_QUAD_PTS WRITE (6,"(A,I5)") & "DEBUG MODE (0-OFF, 1-ON).....................", & DEBUG_MODE ! Read the second line of control data READ (5,"(16I5)") NO_INT_PRPS_PER_NODE, & NO_REAL_PRPS_PER_NODE, & NO_INT_PRPS_PER_EL, & NO_REAL_PRPS_PER_EL, & NO_INT_PRPS_MISC, & NO_REAL_PRPS_MISC, & NODAL_PRPS_FLAG, & ELEMENTAL_PRPS_FLAG, & NODAL_PRINT_FLAG, & ELEMENTAL_PRINT_FLAG, & POST_PROC_IO_UNIT1, & POST_PROC_IO_UNIT2, & NO_FLUX_COMPS, & EL_COL_FLAG, & NO_INT_PRPS_PER_SEG, & NO_REAL_PRPS_PER_SEG ! Echo selected items back to the user WRITE (6,"(A,I5)") & "NUMBER OF FLUX COMPONENTS PER NODE...........", & NO_FLUX_COMPS WRITE (6,"(A,I5)") & "NUMBER OF INTEGER PROPERTIES PER NODE........", & NO_INT_PRPS_PER_NODE WRITE (6,"(A,I5)") & "NUMBER OF REAL PROPERTIES PER NODE...........", & NO_REAL_PRPS_PER_NODE WRITE (6,"(A,I5)") & "NUMBER OF INTEGER PROPERTIES PER ELEMENT.....", & NO_INT_PRPS_PER_EL WRITE (6,"(A,I5)") & "NUMBER OF REAL PROPERTIES PER ELEMENT........", & NO_REAL_PRPS_PER_EL WRITE (6,"(A,I5)") & "NUMBER OF INTEGER PROPERTIES PER SEGMENT.....", & NO_INT_PRPS_PER_SEG WRITE (6,"(A,I5)") & "NUMBER OF REAL PROPERTIES PER SEGMENT........", & NO_REAL_PRPS_PER_SEG WRITE (6,"(A,I5)") & "NUMBER OF INTEGER MISCELLANEOUS PROPERTIES...", & NO_INT_PRPS_MISC WRITE (6,"(A,I5)") & "NUMBER OF REAL MISCELLANEOUS PROPERTIES......", & NO_REAL_PRPS_MISC ! Read the third line of control data READ (5,"(4I5,F4.1)") ADAPT_FLAG, & PTCH_DFN_FLAG, & PLYNML_FLAG, & MAX_ITERS, & PRCNT_CUTOFF ! Echo selected items back to the user WRITE (6,"(A,I5)") & "ADAPTIVITY FLAG (0=NO ADAPTIVITY)............", & ADAPT_FLAG IF (ADAPT_FLAG > 0) THEN WRITE (6,"(A,I5)") & "PATCH DEFINITION METHOD......................", & PTCH_DFN_FLAG WRITE (6,"(A,I5)") & "SPR LEAST SQUARES POLYNOMIAL TYPE............", & PLYNML_FLAG WRITE (6,"(A,I5)") & "MAXIMUM ADAPTIVE ITERATIONS..................", & MAX_ITERS WRITE (6,"(A,F4.1)") & "PERCENT CUTOFF...............................", & PRCNT_CUTOFF END IF ! Process the user remark lines WRITE (6,"(/,A,I5,A)") & "THE NEXT ", NO_USER_REMARKS, & " LINES ARE USER SUPPLIED" DO LINE_CNT = 1, NO_USER_REMARKS READ (5,"(A80)") REMARK_LINE WRITE (6,"(A80)") REMARK_LINE END DO ! LINE_CNT NEXT_UNUSED_NODE = NO_SYS_GEOM_NODES + 1 END SUBROUTINE READ_CONTROL_DATA SUBROUTINE RECOVER_SPRCNVRGNT_DERIVS (NO_SYS_AN_NODES,EL_LST,SPR_DATA_LST, & SPRCNVRGNT_NDL_DRVS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE USE SPR_TYPES_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_AN_NODES TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: EL_LST TYPE (SPR_DATA_LST_PTR_TYPE), INTENT (IN) :: SPR_DATA_LST REAL (DP), DIMENSION (NO_SYS_AN_NODES,NO_ROWS_B) :: SPRCNVRGNT_NDL_DRVS TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR INTEGER, PARAMETER :: MAX_PTCH_ELS = 9 INTEGER, DIMENSION (MAX_PTCH_ELS) :: PTCH_IDS INTEGER :: M INTEGER :: DGR INTEGER :: N INTEGER :: MAX_NODES_R INTEGER :: MAX_NODES_S REAL (DP), DIMENSION (:,:), ALLOCATABLE :: P REAL (DP), DIMENSION (:,:), ALLOCATABLE :: SGMA_H REAL (DP), DIMENSION (:,:), ALLOCATABLE :: A REAL (DP), DIMENSION (:), ALLOCATABLE :: W REAL (DP), DIMENSION (:,:), ALLOCATABLE :: V INTEGER :: INDEX INTEGER :: CMPNT INTEGER, DIMENSION (NO_SYS_AN_NODES) :: CNTS REAL (DP), DIMENSION (NO_SYS_AN_NODES,NO_ROWS_B) :: SUMS INTEGER :: GET_NO_PTCH_GAUSS_PTS INTEGER :: GET_HGHST_PTCH_DGR REAL (DP) :: EPSILON = 1.0E-5 CNTS = 0 SUMS = 0.d0 EL_PTR => EL_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_PTR)) EXIT print*,"ELEMENT = ",EL_PTR % ID,"************" CALL GET_ELS_IN_PATCH (EL_LST,EL_PTR,PTCH_DFN_FLAG,MAX_PTCH_ELS,PTCH_IDS) M = GET_NO_PTCH_GAUSS_PTS (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS) print*,"m = ",m IF (PLYNML_FLAG == 0) THEN DGR = GET_HGHST_PTCH_DGR (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS) N = INT ((DGR + 1.0)*(DGR/2.0 + 1.0)) print*,"n = ",n ALLOCATE (P(M,N)) ALLOCATE (SGMA_H(M,NO_ROWS_B)) CALL PREPARE_COMPLETE_LS_PRBLM (M,N,DGR,SPR_DATA_LST,EL_LST, & MAX_PTCH_ELS,PTCH_IDS,P,SGMA_H) ELSE IF (PLYNML_FLAG == 1) THEN CALL GET_MAX_NODES_IN_R_S (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS, & MAX_NODES_R,MAX_NODES_S) N = MAX_NODES_R*2 + MAX_NODES_S*2 - 4 ALLOCATE (P(M,N)) ALLOCATE (SGMA_H(M,NO_ROWS_B)) CALL PREPARE_SHPFNCTN_LS_PRBLM (M,N,MAX_NODES_R,MAX_NODES_S, & EL_LST,SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS,P,SGMA_H) ELSE DGR = GET_HGHST_PTCH_DGR (SPR_DATA_LST,MAX_PTCH_ELS,PTCH_IDS) + 1 N = INT ((DGR + 1.0)*(DGR/2.0 + 1.0)) ALLOCATE (P(M,N)) ALLOCATE (SGMA_H(M,NO_ROWS_B)) CALL PREPARE_COMPLETE_LS_PRBLM (M,N,DGR,SPR_DATA_LST,EL_LST, & MAX_PTCH_ELS,PTCH_IDS,P,SGMA_H) END IF ALLOCATE (A(N,NO_ROWS_B)) ALLOCATE (W(N)) ALLOCATE (V(N,N)) CALL SVDCMP (P,M,N,M,N,W,V) DO INDEX = 1,N IF (W(INDEX) < EPSILON) W(INDEX) = 0.d0 END DO DO CMPNT = 1,NO_ROWS_B CALL SVBKSB (P,W,V,M,N,M,N,SGMA_H(1,CMPNT),A(1,CMPNT)) END DO IF (PLYNML_FLAG == 0 .OR. PLYNML_FLAG == 2) THEN CALL EVALUATE_COMPLETE_NDL_DRVS (N,DGR,MAX_PTCH_ELS,NO_SYS_AN_NODES, & EL_LST,PTCH_IDS,A,CNTS,SUMS) ELSE CALL EVALUATE_SHPFNCTN_NDL_DRVS (N,MAX_NODES_R,MAX_NODES_S, & MAX_PTCH_ELS,NO_SYS_AN_NODES,EL_LST,PTCH_IDS,A,CNTS,SUMS) END IF DEALLOCATE (P) DEALLOCATE (SGMA_H) DEALLOCATE (A) DEALLOCATE (W) DEALLOCATE (V) EL_PTR => EL_PTR % NEXT END DO DO INDEX = 1,NO_SYS_AN_NODES DO CMPNT = 1,NO_ROWS_B SPRCNVRGNT_NDL_DRVS(INDEX,CMPNT) = SUMS(INDEX,CMPNT)/CNTS(INDEX) END DO END DO END SUBROUTINE RECOVER_SPRCNVRGNT_DERIVS SUBROUTINE REINIT_EDGE (EDG_PTR) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_NDE_TYPE), TARGET :: EDG_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: TEMP_PTR TYPE (NODE_SUBLST_NDE_TYPE), POINTER :: FREE_PTR TEMP_PTR => EDG_PTR % NS_HDR % RPTR DO NULLIFY (FREE_PTR) IF (ASSOCIATED (TEMP_PTR,EDG_PTR % NS_HDR)) EXIT IF (TEMP_PTR % NPTR % DOF_INFO /= 0) THEN IF (TEMP_PTR % NPTR % DOF_INFO < 0) THEN TEMP_PTR % NPTR % DOF_INFO = 0 ELSE ! Delete the node from the node list IF (.NOT. ASSOCIATED (TEMP_PTR % NPTR % NEXT_PTR)) THEN NULLIFY (TEMP_PTR % NPTR % PREV_PTR % NEXT_PTR) NULLIFY (TEMP_PTR % NPTR % PREV_PTR) DEALLOCATE (TEMP_PTR % NPTR) NULLIFY (TEMP_PTR % NPTR) ELSE TEMP_PTR % NPTR % PREV_PTR % NEXT_PTR => & TEMP_PTR % NPTR % NEXT_PTR TEMP_PTR % NPTR % NEXT_PTR % PREV_PTR => & TEMP_PTR % NPTR % PREV_PTR NULLIFY (TEMP_PTR % NPTR % PREV_PTR) NULLIFY (TEMP_PTR % NPTR % NEXT_PTR) DEALLOCATE (TEMP_PTR % NPTR) NULLIFY (TEMP_PTR % NPTR) END IF ! Mark sublist node for deletion FREE_PTR => TEMP_PTR END IF END IF TEMP_PTR => TEMP_PTR % RPTR IF (ASSOCIATED (FREE_PTR)) THEN FREE_PTR % LPTR % RPTR => FREE_PTR % RPTR FREE_PTR % RPTR % LPTR => FREE_PTR % LPTR NULLIFY (FREE_PTR % LPTR) NULLIFY (FREE_PTR % RPTR) DEALLOCATE (FREE_PTR) END IF END DO END SUBROUTINE REINIT_EDGE FUNCTION RESET_SYS_DOF_INFO (NODE_LST) RESULT (NO_SYS_DOFS) USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER :: NO_SYS_DOFS INTEGER :: NO_AN_NODES TYPE (NODE_LST_NDE_TYPE), POINTER :: TEMP_PTR NO_AN_NODES = 0 TEMP_PTR => NODE_LST % HDR % NEXT_PTR DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT ! The node is purely an analysis node IF (TEMP_PTR % DOF_INFO > 0) THEN NO_AN_NODES = NO_AN_NODES + 1 TEMP_PTR % DOF_INFO = NO_AN_NODES END IF ! The node is both a geometry node and an analysis node IF (TEMP_PTR % DOF_INFO < 0) THEN NO_AN_NODES = NO_AN_NODES + 1 TEMP_PTR % DOF_INFO = -NO_AN_NODES END IF TEMP_PTR => TEMP_PTR % NEXT_PTR END DO NO_SYS_DOFS = NO_AN_NODES * NO_PARAMS_PER_NODE END FUNCTION RESET_SYS_DOF_INFO SUBROUTINE SET_ADJACENT_EL_PTRS (ELEMENT_LST) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (ELEMENT_LST_PTR_TYPE), INTENT (IN) :: ELEMENT_LST TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: EL_PTR TYPE (ELEMENT_LST_NDE_TYPE), POINTER :: TEMP_PTR EL_PTR => ELEMENT_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (EL_PTR)) EXIT TEMP_PTR => ELEMENT_LST % HDR % NEXT DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF (EL_PTR % ID /= TEMP_PTR % ID) THEN IF (ASSOCIATED (EL_PTR % EDG_1,TEMP_PTR % EDG_1) .OR. & ASSOCIATED (EL_PTR % EDG_1,TEMP_PTR % EDG_2) .OR. & ASSOCIATED (EL_PTR % EDG_1,TEMP_PTR % EDG_3) .OR. & ASSOCIATED (EL_PTR % EDG_1,TEMP_PTR % EDG_4)) THEN EL_PTR % ADJ_1 => TEMP_PTR ELSE IF (ASSOCIATED (EL_PTR % EDG_2,TEMP_PTR % EDG_1) .OR. & ASSOCIATED (EL_PTR % EDG_2,TEMP_PTR % EDG_2) .OR. & ASSOCIATED (EL_PTR % EDG_2,TEMP_PTR % EDG_3) .OR. & ASSOCIATED (EL_PTR % EDG_2,TEMP_PTR % EDG_4)) THEN EL_PTR % ADJ_2 => TEMP_PTR ELSE IF (ASSOCIATED (EL_PTR % EDG_3,TEMP_PTR % EDG_1) .OR. & ASSOCIATED (EL_PTR % EDG_3,TEMP_PTR % EDG_2) .OR. & ASSOCIATED (EL_PTR % EDG_3,TEMP_PTR % EDG_3) .OR. & ASSOCIATED (EL_PTR % EDG_3,TEMP_PTR % EDG_4)) THEN EL_PTR % ADJ_3 => TEMP_PTR ELSE IF (ASSOCIATED (EL_PTR % EDG_4,TEMP_PTR % EDG_1) .OR. & ASSOCIATED (EL_PTR % EDG_4,TEMP_PTR % EDG_2) .OR. & ASSOCIATED (EL_PTR % EDG_4,TEMP_PTR % EDG_3) .OR. & ASSOCIATED (EL_PTR % EDG_4,TEMP_PTR % EDG_4)) THEN EL_PTR % ADJ_4 => TEMP_PTR END IF END IF TEMP_PTR => TEMP_PTR % NEXT END DO EL_PTR => EL_PTR % NEXT END DO END SUBROUTINE SET_ADJACENT_EL_PTRS SUBROUTINE SKYFAC (N, LD, SINGAB, PDCHEK, ICOND, A, ACOND, DETCF, & IDETEX, NEGEIG, IFAIL) USE DBL_PRECISION_MODULE IMPLICIT NONE INTERFACE FUNCTION DOTPRD (VCTR_1, VCTR_2, LNGTH) RESULT (VALUE) USE DBL_PRECISION_MODULE IMPLICIT NONE REAL (DP), DIMENSION (:), INTENT (IN) :: VCTR_1 REAL (DP), DIMENSION (:), INTENT (IN) :: VCTR_2 INTEGER, INTENT (IN) :: LNGTH REAL (DP) :: VALUE END FUNCTION DOTPRD END INTERFACE INTEGER, INTENT (IN) :: N INTEGER, DIMENSION (N+1), INTENT (IN) :: LD LOGICAL, INTENT (IN) :: SINGAB LOGICAL, INTENT (IN) :: PDCHEK INTEGER, INTENT (IN) :: ICOND REAL (DP), DIMENSION (LD(N+1)), INTENT (INOUT) :: A REAL (DP), INTENT (OUT) :: ACOND REAL (DP), INTENT (OUT) :: DETCF INTEGER, INTENT (OUT) :: IDETEX INTEGER, INTENT (OUT) :: NEGEIG INTEGER, INTENT (OUT) :: IFAIL REAL (DP), DIMENSION (:), ALLOCATABLE :: V INTEGER :: NBEG INTEGER :: NEND LOGICAL :: SINGULAR LOGICAL :: INDEF REAL (DP) :: EPSMAC INTEGER :: NBEGP1 INTEGER :: II INTEGER :: M INTEGER :: L REAL (DP) :: D INTEGER :: JK INTEGER :: JMJ INTEGER :: KU INTEGER :: IJ REAL (DP) :: TOLROW REAL (DP) :: DMAX REAL (DP) :: AIJ2 REAL (DP) :: NWA INTEGER :: JJ INTEGER :: I INTEGER :: J INTEGER :: K REAL (DP) :: DELTA ! Allocate the scratch array IF (ICOND == 0) THEN ALLOCATE (V(N)) ELSE ALLOCATE (V(4*N)) END IF ! Initialization NBEG = 0 NEND = N IFAIL = 0 IDETEX = 0 NEGEIG = 0 SINGULAR = .FALSE. INDEF = .FALSE. EPSMAC = TINY (EPSMAC) NWA = ABS (LD(N+1)) ! Compute squared lengths of unconstrained rows NEGB+1 thru N NBEGP1 = NBEG + 1 DO I = NBEGP1,N II = LD(I+1) IF (II > 0) THEN V(I) = A(II)**2 M = II - I K = MAX (NBEGP1,ABS(LD(I))-M+1) L = MIN (NEND,I) - 1 IF (K-L <= 0) THEN DO J = K,L IF (LD(J+1) > 0) THEN AIJ2 = A(M+J)**2 V(I) = V(I) + AIJ2 V(J) = V(J) + AIJ2 END IF END DO ! J END IF END IF END DO ! I ! Factorization DO J = NBEGP1,NEND JJ = LD(J+1) IF (JJ > 0) THEN D = A(JJ) JMJ = ABS(LD(J)) JK = JJ - JMJ KU = JK - 1 IF (KU /= 0) THEN DO K = 1,KU I = J - JK + K V(K) = 0.d0 II = LD(I+1) IF (II > 0) THEN M = MIN (II-ABS(LD(I)),K) - 1 IJ = JMJ + K V(K) = A(IJ) - DOTPRD (A(II-M:II-1),V(K-M:K-1),M) A(IJ) = V(K) * A(II) END IF END DO ! K D = D - DOTPRD (A(JMJ+1:JMJ+KU),V(1:KU),KU) END IF ! Singularity test TOLROW = 8.0 * EPSMAC * SQRT(V(J)) IF (ABS(D) <= TOLROW) THEN IF (SINGAB) THEN SINGULAR = .TRUE. EXIT ELSE D = TOLROW ENDIF END IF A(JJ) = 1.d0/D ! Update determinant if DETCF is nonzero IF (DETCF /= 0.d0) THEN DETCF = DETCF * D DO IF (ABS(DETCF) < 1.d0) EXIT DETCF = DETCF * 0.0625 IDETEX = IDETEX + 4 END DO DO IF (ABS(DETCF) >= 0.0625) EXIT DETCF = DETCF * 16.d0 IDETEX = IDETEX - 4 END DO END IF ! Positive definitness check if PDCHEK is .TRUE. IF (D <= 0.d0) THEN NEGEIG = NEGEIG + 1 IF (PDCHEK) THEN INDEF = .TRUE. EXIT END IF END IF END IF END DO ! J IF (.NOT.SINGULAR .AND. .NOT.INDEF) THEN IF (ICOND /= 0) THEN K = 0 DMAX = 0.d0 DO I = 1,ICOND ! CALL SKYRDM (V, NWA, 0.0, K) ! CALL SKYSOL DMAX = MAX (DELTA, DMAX) END DO ACOND = DMAX / ((1.d0+DMAX)*EPSMAC) END IF ELSE IF (SINGULAR) THEN IFAIL = J DETCF = 0.d0 ELSE IFAIL = -J END IF END IF DEALLOCATE (V) END SUBROUTINE SKYFAC SUBROUTINE SKYRDM (V, N, VMEAN, IRDM) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: N REAL (DP), DIMENSION (4*N), INTENT (OUT) :: V REAL (DP), INTENT (IN) :: VMEAN INTEGER, INTENT (INOUT) :: IRDM REAL (DP) :: C INTEGER :: IX REAL (DP) :: X INTEGER :: I C = VMEAN - 0.5 IF (IRDM == 0) THEN IX = 26903 IRDM = 1 END IF DO I = 1,N IX = MOD (27181*IX+13849,65536) X = IX V(I) = 0.15258789E-4 * X + C END DO END SUBROUTINE SKYRDM SUBROUTINE SKYSOL (N, LD, A, B, X) USE DBL_PRECISION_MODULE IMPLICIT NONE INTERFACE FUNCTION DOTPRD (VCTR_1, VCTR_2, LNGTH) RESULT (VALUE) USE DBL_PRECISION_MODULE IMPLICIT NONE REAL (DP), DIMENSION (:), INTENT (IN) :: VCTR_1 REAL (DP), DIMENSION (:), INTENT (IN) :: VCTR_2 INTEGER, INTENT (IN) :: LNGTH REAL (DP) :: VALUE END FUNCTION DOTPRD END INTERFACE INTEGER, INTENT (IN) :: N INTEGER, DIMENSION (N+1), INTENT (IN) :: LD REAL (DP), DIMENSION (LD(N+1)), INTENT (IN) :: A REAL (DP), DIMENSION (N), INTENT (IN) :: B REAL (DP), DIMENSION (N), INTENT (OUT) :: X REAL (DP) :: BXFAC INTEGER :: I INTEGER :: II INTEGER :: IMI INTEGER :: M INTEGER :: K INTEGER :: J ! Initialization BXFAC = 1.d0 DO I = 1,N X(I) = B(I) END DO ! Forward substitution pass DO I = 1,N II = LD(I+1) IF (II <= 0) THEN X(I) = 0.d0 ELSE IMI = ABS (LD(I)) M = II - IMI - 1 X(I) = X(I) - DOTPRD (A(IMI+1:IMI+M),X(I-M:I-1),M) END IF END DO ! Scaling pass DO I = 1,N II = ABS (LD(I+1)) X(I) = A(II) * X(I) END DO ! Backsubstitution pass I = N DO K = 1,N II = LD(I+1) IF (II <= 0) THEN X(I) = BXFAC * B(I) ELSE M = II - ABS (LD(I)) - 1 IF (M /= 0) THEN DO J = 1,M X(I-J) = X(I-J) - A(II-J) * X(I) END DO END IF END IF I = I - 1 END DO END SUBROUTINE SKYSOL MODULE SPR_TYPES_MODULE USE CONTROL_MODULE USE DBL_PRECISION_MODULE INTEGER, PARAMETER :: TEMP = 2 TYPE SMPLNG_DATA_LST_NDE_TYPE REAL (DP) :: X REAL (DP) :: Y REAL (DP), DIMENSION (TEMP) :: DERIVS TYPE (SMPLNG_DATA_LST_NDE_TYPE), POINTER :: NEXT END TYPE SMPLNG_DATA_LST_NDE_TYPE TYPE EL_DATA_LST_NDE_TYPE INTEGER :: EL_ID INTEGER :: MAX_NDS_R INTEGER :: MAX_NDS_S INTEGER :: NO_SP TYPE (SMPLNG_DATA_LST_NDE_TYPE) :: SD_HDR TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: NEXT END TYPE EL_DATA_LST_NDE_TYPE TYPE SPR_DATA_LST_PTR_TYPE TYPE (EL_DATA_LST_NDE_TYPE), POINTER :: HDR END TYPE SPR_DATA_LST_PTR_TYPE END MODULE SPR_TYPES_MODULE SUBROUTINE STORE_FOR_POST_PRCSSNG () USE CONTROL_MODULE USE DATA_STRUCTURES_MODULE USE ELMNT_ASSMBLY_MODULE IMPLICIT NONE END SUBROUTINE STORE_FOR_POST_PRCSSNG SUBROUTINE SVBKSB (U, W, V, M, N, MP, NP, B, X) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: M INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: MP INTEGER, INTENT (IN) :: NP REAL (DP), DIMENSION (MP,NP), INTENT (IN) :: U REAL (DP), DIMENSION (NP), INTENT (IN) :: W REAL (DP), DIMENSION (NP,NP), INTENT (IN) :: V REAL (DP), DIMENSION (MP), INTENT (IN) :: B REAL (DP), DIMENSION (NP), INTENT (OUT) :: X INTEGER, PARAMETER :: NMAX = 100 REAL (DP), DIMENSION (NMAX) :: TMP REAL (DP) :: S INTEGER :: I INTEGER :: J INTEGER :: JJ DO J = 1,N S = 0.d0 IF (W(J) /= 0.d0) THEN DO I = 1,M S = S+U(I,J)*B(I) END DO S = S/W(J) END IF TMP(J) = S END DO DO J = 1,N S = 0.d0 DO JJ = 1,N S = S+V(J,JJ)*TMP(JJ) END DO X(J) = S END DO END SUBROUTINE SVBKSB SUBROUTINE SVDCMP (A, M, N, MP, NP, W, V) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: M INTEGER, INTENT (IN) :: N INTEGER, INTENT (IN) :: MP INTEGER, INTENT (IN) :: NP REAL (DP), DIMENSION (MP,NP), INTENT (INOUT) :: A REAL (DP), DIMENSION (NP), INTENT (OUT) :: W REAL (DP), DIMENSION (NP,NP), INTENT (OUT) :: V INTEGER, PARAMETER :: NMAX = 100 REAL (DP), DIMENSION (NMAX) :: RV1 REAL (DP) :: G REAL (DP) :: SCALE REAL (DP) :: ANORM REAL (DP) :: S REAL (DP) :: F REAL (DP) :: H REAL (DP) :: C REAL (DP) :: X REAL (DP) :: Y REAL (DP) :: Z INTEGER :: ITS INTEGER :: I INTEGER :: J INTEGER :: JJ INTEGER :: K INTEGER :: L INTEGER :: NM IF (M < N) THEN WRITE (6,"(//,A)") "FATAL ERROR" WRITE (6,"(A)") "Fewer equations than unknowns in routine SVDCMP" STOP END IF ! Householder reduction to bidiagonal form G = 0.d0 SCALE = 0.d0 ANORM = 0.d0 DO I = 1,N L = I+1 RV1(I) = SCALE*G G = 0.d0 S = 0.d0 SCALE = 0.d0 IF (I <= M) THEN DO K = I,M SCALE = SCALE+ABS(A(K,I)) END DO IF (SCALE /= 0.d0) THEN DO K = I,M A(K,I) = A(K,I)/SCALE S = S+A(K,I)*A(K,I) END DO F = A(I,I) G = -SIGN(SQRT(S),F) H = F*G-S A(I,I) = F-G IF (I /= N) THEN DO J = L,N S = 0.d0 DO K = I,M S = S+A(K,I)*A(K,J) END DO F = S/H DO K = I,M A(K,J) = A(K,J)+F*A(K,I) END DO END DO END IF DO K = I,M A(K,I) = SCALE*A(K,I) END DO END IF END IF W(I) = SCALE*G G = 0.d0 S = 0.d0 SCALE = 0.d0 IF ((I <= M) .AND. (I /= N)) THEN DO K = L,N SCALE = SCALE+ABS(A(I,K)) END DO IF (SCALE /= 0.d0) THEN DO K = L,N A(I,K) = A(I,K)/SCALE S = S+A(I,K)*A(I,K) END DO F = A(I,L) G = -SIGN(SQRT(S),F) H = F*G-S A(I,L) = F-G DO K = L,N RV1(K) = A(I,K)/H END DO IF (I /= M) THEN DO J = L,M S = 0.d0 DO K = L,N S = S+A(J,K)*A(I,K) END DO DO K = L,N A(J,K) = A(J,K)+S*RV1(K) END DO END DO END IF DO K = L,N A(I,K) = SCALE*A(I,K) END DO END IF END IF ANORM = MAX(ANORM,(ABS(W(I))+ABS(RV1(I)))) END DO ! Accumulation of right-hand transformations DO I = N,1,-1 IF (I < N) THEN IF (G /= 0.d0) THEN DO J = L,N V(J,I) = (A(I,J)/A(I,L))/G END DO DO J = L,N S = 0.d0 DO K = L,N S = S+A(I,K)*V(K,J) END DO DO K = L,N V(K,J) = V(K,J)+S*V(K,I) END DO END DO END IF DO J = L,N V(I,J) = 0.d0 V(J,I) = 0.d0 END DO END IF V(I,I) = 1.d0 G = RV1(I) L = I END DO ! Accumulation of left-hand transformations DO I = N,1,-1 L = I+1 G = W(I) IF (I < N) THEN DO J = L,N A(I,J) = 0.d0 END DO END IF IF (G /= 0.d0) THEN G = 1.d0/G IF (I /= N) THEN DO J = L,N S = 0.d0 DO K = L,M S = S+A(K,I)*A(K,J) END DO F = (S/A(I,I))*G DO K = I,M A(K,J) = A(K,J)+F*A(K,I) END DO END DO END IF DO J = I,M A(J,I) = A(J,I)*G END DO ELSE DO J = I,M A(J,I) = 0.d0 END DO END IF A(I,I) = A(I,I)+1.d0 END DO ! Diagonalization of the bidiagonal form DO K = N,1,-1 DO ITS = 1,30 DO L = K,1,-1 NM = L-1 IF ((ABS(RV1(L))+ANORM) == ANORM) GOTO 2 IF ((ABS(W(NM))+ANORM) == ANORM) GOTO 1 END DO 1 CONTINUE C = 0.d0 S = 1.d0 DO I = L,K F = S*RV1(I) RV1(I) = C*RV1(I) IF ((ABS(F)+ANORM) == ANORM) GOTO 2 G = W(I) H = SQRT(F*F+G*G) W(I) = H H = 1.d0/H C = G*H S = -(F*H) DO J = 1,M Y = A(J,NM) Z = A(J,I) A(J,NM) = (Y*C)+(Z*S) A(J,I) = -(Y*S)+(Z*C) END DO END DO 2 CONTINUE Z = W(K) IF (L == K) THEN IF (Z < 0.d0) THEN W(K) = -Z DO J = 1,N V(J,K) = -V(J,K) END DO END IF GOTO 3 END IF IF (ITS == 30) WRITE (6,"(A)") "No convergence in 30 iterations" X = W(L) NM = K-1 Y = W(NM) G = RV1(NM) H = RV1(K) F = ((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y) G = SQRT(F*F+1.d0) F = ((X-Z)*(X+Z)+H*((Y/(F+SIGN(G,F)))-H))/X ! Next QR transformation C = 1.d0 S = 1.d0 DO J = L,NM I = J+1 G = RV1(I) Y = W(I) H = S*G G = C*G Z = SQRT(F*F+H*H) RV1(J) = Z C = F/Z S = H/Z F = (X*C)+(G*S) G = -(X*S)+(G*C) H = Y*S Y = Y*C DO JJ = 1,N X = V(JJ,J) Z = V(JJ,I) V(JJ,J) = (X*C)+(Z*S) V(JJ,I) = -(X*S)+(Z*C) END DO Z = SQRT(F*F+H*H) W(J) = Z IF (Z /= 0.d0) THEN Z = 1.d0/Z C = F*Z S = H*Z END IF F = (C*G)+(S*Y) X = -(S*G)+(C*Y) DO JJ = 1,M Y = A(JJ,J) Z = A(JJ,I) A(JJ,J) = (Y*C)+(Z*S) A(JJ,I) = -(Y*S)+(Z*C) END DO END DO RV1(L) = 0.d0 RV1(K) = F W(K) = X END DO 3 CONTINUE END DO END SUBROUTINE SVDCMP MODULE SYSTEM_EQNS_MODULE USE DBL_PRECISION_MODULE REAL (DP), DIMENSION (:), ALLOCATABLE :: SYS_SQR_MTRX INTEGER, DIMENSION (:), ALLOCATABLE :: SYS_COL_HGHTS INTEGER, DIMENSION (:), ALLOCATABLE :: DIAG_LCTNS REAL (DP), DIMENSION (:), ALLOCATABLE :: SYS_SRC_VCTR REAL (DP), DIMENSION (:), ALLOCATABLE :: SYS_SOL_VCTR END MODULE SYSTEM_EQNS_MODULE SUBROUTINE TAG_ANALYSIS_NODES (NODE_LST, CORNERS) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER, DIMENSION (1:4), INTENT (IN) :: CORNERS TYPE (NODE_LST_NDE_TYPE), POINTER :: TEMP_PTR TEMP_PTR => NODE_LST % HDR % NEXT_PTR DO IF (.NOT. ASSOCIATED (TEMP_PTR)) EXIT IF ((TEMP_PTR % ID .EQ. CORNERS(1)) .OR. & (TEMP_PTR % ID .EQ. CORNERS(2)) .OR. & (TEMP_PTR % ID .EQ. CORNERS(3)) .OR. & (TEMP_PTR % ID .EQ. CORNERS(4))) & TEMP_PTR % DOF_INFO = -1 TEMP_PTR => TEMP_PTR % NEXT_PTR END DO END SUBROUTINE TAG_ANALYSIS_NODES SUBROUTINE UPDATE_EDGE_LIST (NODE_LST, EDGE_LST, CORNERS, SIDES, & NEXT_UNUSED_EDGE_ID) USE DATA_STRUCTURES_MODULE IMPLICIT NONE INTERFACE SUBROUTINE INSERT_EDGE (EDGE_LST, EDGE_ID, END_1, MID, END_2) USE DATA_STRUCTURES_MODULE IMPLICIT NONE TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST INTEGER, INTENT (IN) :: EDGE_ID TYPE (NODE_LST_NDE_TYPE), TARGET :: END_1 TYPE (NODE_LST_NDE_TYPE), POINTER :: MID TYPE (NODE_LST_NDE_TYPE), TARGET :: END_2 END SUBROUTINE INSERT_EDGE END INTERFACE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST TYPE (EDGE_LST_PTR_TYPE), INTENT (IN) :: EDGE_LST INTEGER, DIMENSION (1:4), INTENT (IN) :: CORNERS INTEGER, DIMENSION (1:4), INTENT (IN) :: SIDES INTEGER, INTENT (INOUT) :: NEXT_UNUSED_EDGE_ID INTEGER :: EDGE_ID TYPE (NODE_LST_NDE_TYPE), POINTER :: END_1 TYPE (NODE_LST_NDE_TYPE), POINTER :: END_2 TYPE (NODE_LST_NDE_TYPE), POINTER :: MID INTEGER :: GET_NEW_EDGE_ID TYPE (NODE_LST_NDE_TYPE), POINTER :: GET_NODE_PTR ! Process the first edge EDGE_ID = GET_NEW_EDGE_ID (EDGE_LST, CORNERS(1), CORNERS(2), & NEXT_UNUSED_EDGE_ID) IF (EDGE_ID .NE. 0) THEN END_1 => GET_NODE_PTR (NODE_LST, CORNERS(1)) END_2 => GET_NODE_PTR (NODE_LST, CORNERS(2)) IF (SIDES(1) .NE. 0) THEN MID => GET_NODE_PTR (NODE_LST, SIDES(1)) ELSE NULLIFY (MID) ENDIF CALL INSERT_EDGE (EDGE_LST, EDGE_ID, END_1, MID, END_2) ENDIF ! Process the second edge EDGE_ID = GET_NEW_EDGE_ID (EDGE_LST, CORNERS(2), CORNERS(3), & NEXT_UNUSED_EDGE_ID) IF (EDGE_ID .NE. 0) THEN END_1 => GET_NODE_PTR (NODE_LST, CORNERS(2)) END_2 => GET_NODE_PTR (NODE_LST, CORNERS(3)) IF (SIDES(2) .NE. 0) THEN MID => GET_NODE_PTR (NODE_LST, SIDES(2)) ELSE NULLIFY (MID) ENDIF CALL INSERT_EDGE (EDGE_LST, EDGE_ID, END_1, MID, END_2) ENDIF ! Process the third edge EDGE_ID = GET_NEW_EDGE_ID (EDGE_LST, CORNERS(3), CORNERS(4), & NEXT_UNUSED_EDGE_ID) IF (EDGE_ID .NE. 0) THEN END_1 => GET_NODE_PTR (NODE_LST, CORNERS(3)) END_2 => GET_NODE_PTR (NODE_LST, CORNERS(4)) IF (SIDES(3) .NE. 0) THEN MID => GET_NODE_PTR (NODE_LST, SIDES(3)) ELSE NULLIFY (MID) ENDIF CALL INSERT_EDGE (EDGE_LST, EDGE_ID, END_1, MID, END_2) ENDIF ! Process the fourth edge EDGE_ID = GET_NEW_EDGE_ID (EDGE_LST, CORNERS(4), CORNERS(1), & NEXT_UNUSED_EDGE_ID) IF (EDGE_ID .NE. 0) THEN END_1 => GET_NODE_PTR (NODE_LST, CORNERS(4)) END_2 => GET_NODE_PTR (NODE_LST, CORNERS(1)) IF (SIDES(4) .NE. 0) THEN MID => GET_NODE_PTR (NODE_LST, SIDES(4)) ELSE NULLIFY (MID) ENDIF CALL INSERT_EDGE (EDGE_LST, EDGE_ID, END_1, MID, END_2) ENDIF END SUBROUTINE UPDATE_EDGE_LIST SUBROUTINE UPDATE_NDL_DRV_AVRG_DATA (NO_SYS_AN_NODES,N,P,A,INDX,CNTS,SUMS) USE CONTROL_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_SYS_AN_NODES INTEGER, INTENT (IN) :: N REAL (DP), DIMENSION (N), INTENT (IN) :: P REAL (DP), DIMENSION (N,NO_ROWS_B), INTENT (IN) :: A INTEGER, INTENT (IN) :: INDX INTEGER, DIMENSION (NO_SYS_AN_NODES),INTENT (OUT) :: CNTS REAL (DP), DIMENSION (NO_SYS_AN_NODES,NO_ROWS_B),INTENT (OUT) :: SUMS INTEGER :: CMPNT REAL (DP) :: DRV INTEGER :: PLYNML_TERM DO CMPNT = 1,NO_ROWS_B DRV = 0.d0 DO PLYNML_TERM = 1,N DRV = DRV + P(PLYNML_TERM) * A(PLYNML_TERM,CMPNT) END DO SUMS(INDX,CMPNT) = SUMS(INDX,CMPNT) + DRV END DO CNTS(INDX) = CNTS(INDX) + 1 END SUBROUTINE UPDATE_NDL_DRV_AVRG_DATA SUBROUTINE UPDATE_SYS_SQR_MTRX (NO_EL_DOFS, NO_SYS_DOFS, DOF_MAP, & DIAG_LCTNS, EL_SQR_MTRX, SYS_SQR_MTRX) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_EL_DOFS INTEGER, INTENT (IN) :: NO_SYS_DOFS INTEGER, DIMENSION (NO_EL_DOFS), INTENT (IN) :: DOF_MAP INTEGER, DIMENSION (NO_SYS_DOFS+1), INTENT (IN) :: DIAG_LCTNS REAL (DP), DIMENSION (NO_EL_DOFS,NO_EL_DOFS), INTENT (IN) :: EL_SQR_MTRX REAL (DP), DIMENSION (DIAG_LCTNS(NO_SYS_DOFS+1)), INTENT (OUT) :: SYS_SQR_MTRX INTEGER :: EL_COL INTEGER :: SYS_COL INTEGER :: BASE INTEGER :: EL_ROW INTEGER :: SYS_ROW INTEGER :: INDX DO EL_COL = 1,NO_EL_DOFS SYS_COL = DOF_MAP(EL_COL) BASE = DIAG_LCTNS(SYS_COL+1) - SYS_COL DO EL_ROW = 1,NO_EL_DOFS SYS_ROW = DOF_MAP(EL_ROW) IF (SYS_ROW <= SYS_COL) THEN INDX = BASE + SYS_ROW SYS_SQR_MTRX(INDX) = SYS_SQR_MTRX(INDX) + & EL_SQR_MTRX(EL_ROW,EL_COL) ENDIF END DO END DO END SUBROUTINE UPDATE_SYS_SQR_MTRX SUBROUTINE UPDATE_SYS_SRC_VCTR (NO_EL_DOFS, NO_SYS_DOFS, DOF_MAP, & EL_SRC_VCTR, SYS_SRC_VCTR) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_EL_DOFS INTEGER, INTENT (IN) :: NO_SYS_DOFS INTEGER, DIMENSION (NO_EL_DOFS), INTENT (IN) :: DOF_MAP REAL (DP), DIMENSION (NO_EL_DOFS), INTENT (IN) :: EL_SRC_VCTR REAL (DP), DIMENSION (NO_SYS_DOFS), INTENT (OUT) :: SYS_SRC_VCTR INTEGER :: EL_ROW INTEGER :: SYS_ROW DO EL_ROW = 1,NO_EL_DOFS SYS_ROW = DOF_MAP(EL_ROW) SYS_SRC_VCTR(SYS_ROW) = SYS_SRC_VCTR(SYS_ROW) + & EL_SRC_VCTR(EL_ROW) END DO END SUBROUTINE UPDATE_SYS_SRC_VCTR SUBROUTINE WRITE_IDEAS_NODE_UNV (NODE_LST,ITERATION) USE DATA_STRUCTURES_MODULE USE DBL_PRECISION_MODULE IMPLICIT NONE TYPE (NODE_LST_PTR_TYPE), INTENT (IN) :: NODE_LST INTEGER, INTENT (IN) :: ITERATION CHARACTER (LEN = 12) :: UNV_NAME = "nodes_00.unv" INTEGER :: UNV_UNIT = 30 TYPE (NODE_LST_NDE_TYPE), POINTER :: NODE_PTR INTEGER :: GLOBAL_CS = 0 INTEGER :: WHITE = 15 REAL (DP) :: Z_VALUE = 0.d0 IF (ITERATION < 10) THEN WRITE (UNV_NAME(8:8),"(I1)") ITERATION ELSE WRITE (UNV_NAME(7:8),"(I2)") ITERATION END IF OPEN (UNV_UNIT,FILE=UNV_NAME,STATUS="UNKNOWN") WRITE (UNV_UNIT,"(A)") " -1" WRITE (UNV_UNIT,"(A)") " 15" NODE_PTR => NODE_LST % HDR % NEXT_PTR DO IF (.NOT. ASSOCIATED (NODE_PTR)) EXIT WRITE (UNV_UNIT,"(4I10,1P3E13.5)") NODE_PTR % ID, GLOBAL_CS, & GLOBAL_CS, WHITE, NODE_PTR % X, NODE_PTR % Y, Z_VALUE NODE_PTR => NODE_PTR % NEXT_PTR END DO WRITE (UNV_UNIT,"(A)") " -1" CLOSE (UNV_UNIT,STATUS="KEEP") END SUBROUTINE WRITE_IDEAS_NODE_UNV SUBROUTINE ZERO_OUT_MATRIX (ROWS, COLS, MATRIX) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: ROWS INTEGER, INTENT (IN) :: COLS REAL (DP), DIMENSION (ROWS,COLS), INTENT (OUT) :: MATRIX INTEGER :: ROW_INDEX INTEGER :: COL_INDEX DO ROW_INDEX = 1,ROWS DO COL_INDEX = 1,COLS MATRIX(ROW_INDEX,COL_INDEX) = 0.d0 END DO END DO END SUBROUTINE ZERO_OUT_MATRIX SUBROUTINE ZERO_OUT_VECTOR (NO_OF_ENTRIES, VECTOR) USE DBL_PRECISION_MODULE IMPLICIT NONE INTEGER, INTENT (IN) :: NO_OF_ENTRIES REAL (DP), DIMENSION (1:NO_OF_ENTRIES), INTENT (OUT) :: VECTOR INTEGER :: INDEX DO INDEX = 1, NO_OF_ENTRIES VECTOR(INDEX) = 0.d0 END DO END SUBROUTINE ZERO_OUT_VECTOR