! copyright 2005, J. E. Akin, all rights reserved. ! Begin demo_select_lib.f (also see select_source_lib.f) ! NOTE: This is the set of programs that are valid for one example ! application. Each example is identified in its data file by the ! keyword "example" followed by a number, say nnn. Given nnn the ! corresponding files here are executed by the code in the ! select_source_lib.f using SELECT CASE features. Thus both ! libraries need to be edited when a new example is added to ! this source library. ! Only ELEM_SQ_MATRIX is always required by every example problem. ! It must be edited to contain the problem dependent source, or the ! include path to the problem dependent source. Then this library ! must be compiled and linked to the other archived executables. ! NOTE: Routines B_MATRIX_EX_nnn and E_MATRIX_EX_nnn are required only ! if the SCP error erstimator is invoked. (With nnn replaced with the ! three digit EXAMPLE number.) !============= Begin Files for EXAMPLE number nnn ============= SUBROUTINE DESCRIBE_EXAMPLE_nnn ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE nnn ***' END SUBROUTINE DESCRIBE_EXAMPLE_nnn FUNCTION GET_B_MATRIX_EX_nnn (DGH, XYZ) RESULT (B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT B MATRIX FOR STRAIN OR GRADIENTS ! (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for LT_N, LT_FREE, & H (LT_N) IF NEEDED IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP) :: B (N_R_B, LT_FREE) B = 0.d0 ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *,DGH (1, 1), XYZ (1), B (1, 1) ! STOP 'Edit GET_B_MATRIX_EX_nnn to use my_b_matrix_inc' STOP 'ERROR: NO SOURCE AT GET_B_MATRIX_EX_nnn' END FUNCTION GET_B_MATRIX_EX_nnn SUBROUTINE B_MATRIX_EX_nnn (DGH, XYZ, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT B MATRIX FOR STRAIN OR GRADIENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for LT_N, LT_FREE, & H (LT_N) IF NEEDED IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP), INTENT(OUT) :: B (N_R_B, LT_FREE) B = 0.d0 ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, DGH (1, 1), XYZ (1), B (1, 1) STOP 'ERROR: NO SOURCE AT B_MATRIX_EX_nnn' END SUBROUTINE B_MATRIX_EX_nnn FUNCTION GET_E_MATRIX_EX_nnn (IE, XYZ) RESULT (E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT (IN) :: IE REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP) :: E (N_R_B, N_R_B) E = 0.d0 ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, IE, XYZ (1), E (1, 1) STOP 'ERROR: NO SOURCE AT GET_E_MATRIX_EX_nnn' END FUNCTION GET_E_MATRIX_EX_nnn SUBROUTINE E_MATRIX_EX_nnn (IE, XYZ, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW ! (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT (IN) :: IE REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP), INTENT(OUT) :: E (N_R_B, N_R_B) E = 0.d0 ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, IE, XYZ (1), E (1, 1) STOP 'ERROR: NO SOURCE AT E_MATRIX_EX_nnn' END SUBROUTINE E_MATRIX_EX_nnn SUBROUTINE ELEM_SQ_EX_nnn (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_* IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) STOP 'ERROR: NO SOURCE AT ELEM_SQ_EX_nnn' END SUBROUTINE ELEM_SQ_EX_nnn SUBROUTINE MIXED_SQ_EX_nnn (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE MIXED_BC SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** MIXED_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! suppress compiler warnings (never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) STOP 'ERROR: NO SOURCE AT MIXED_SQ_EX_nnn' END SUBROUTINE MIXED_SQ_EX_nnn SUBROUTINE ELEM_COL_EX_nnn (E, H_INTG, PRT_L_PT, PRT_MAT,& L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE OPTIONAL ELEMENT COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & XYZ (N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! ..................................................... ! *** ELEM_COL_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! ..................................................... ! suppress compiler warnings (never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) STOP 'ERROR: NO SOURCE AT ELEM_COL_EX_nnn' END SUBROUTINE ELEM_COL_EX_nnn SUBROUTINE SEG_COL_EX_nnn (E, H_INTG, PRT_L_PT, PRT_MAT,& L_PT_PROP, IE, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE FLUX SEGMENT COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Geometric_Properties Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! OPTIONAL FLUX VALUES IF BOUNDARY SEGMENT ELEMENT REAL(DP), INTENT(IN) :: FLUX (L_B_N, N_G_FLUX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & XYZ (N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! ..................................................... ! *** SEG_COL_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! ..................................................... ! suppress compiler warnings (never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, FLUX STOP 'ERROR: NO SOURCE AT SEG_COL_EX_nnn' END SUBROUTINE SEG_COL_EX_nnn SUBROUTINE ELEM_POST_EX_nnn (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, FLUX) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE DATA FOR ELEMENT POST-SOLUTION USE IN POST_PROCESS_ELEM ! IF THAT WAS NOT DONE IN ELEM_SQ_MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! OPTIONAL FLUX VALUES IF BOUNDARY SEGMENT ELEMENT REAL(DP), INTENT(IN) :: FLUX (L_B_N, N_G_FLUX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & XYZ (N_SPACE) REAL(DP) :: H_INTG (LT_N ), & B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2), BODY (N_SPACE) REAL(DP) :: DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! ..................................................... ! *** ELEM_POST_EX_nnn PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! ..................................................... ! suppress compiler warnings (never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, FLUX STOP 'ERROR: NO SOURCE AT ELEM_POST_EX_nnn' END SUBROUTINE ELEM_POST_EX_nnn SUBROUTINE POST_PROCESS_ELEM_EX_nnn (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_nnn ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_ELEM_EX_nnn' END SUBROUTINE POST_PROCESS_ELEM_EX_nnn SUBROUTINE POST_PROCESS_MIXED_EX_nnn (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! MIXED SEGMENT POST-SOLUTION CALCULATIONS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, MIXEDENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_MIXED PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_MIXED_EX_nnn' END SUBROUTINE POST_PROCESS_MIXED_EX_nnn FUNCTION START_DOF_VALUE_EX_nnn (IG, XYZ) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! DEFINE STARTING VALUE OF PARAMETER IG IN TERMS OF ! COORDINATES OF THE NODE (FOR ITERATIVE SOLUTIONS) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! A PROBLEM DEPENDENT ROUTINE Use System_Constants ! For EXAMPLE IMPLICIT NONE INTEGER, INTENT(IN) :: IG ! local dof number REAL(DP), INTENT(IN) :: XYZ (N_SPACE) ! position REAL(DP) :: START_DOF_VALUE_EX_nnn ! result ! IG = LOCAL PARAMETER NUMBER AT NODE ! .................................................... ! ** PROBLEM DEPENDENT START STATEMENTS FOLLOW ** ! .................................................... START_DOF_VALUE_EX_nnn = 0.d0 ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *,'START_DOF_VALUE:', IG, XYZ STOP 'ERROR: NO SOURCE AT START_DOF_VALUE_EX_nnn' END FUNCTION START_DOF_VALUE_EX_nnn !SUBROUTINE GET_CONSTANT_E_nnn (IE, E) !b turned off 4/8/02 !! * * * * * * * * * * * * * * * * * * * * * * * * * * * * !! COMPUTE A CONSTANT CONSTITUTIVE MATRIX, E !! * * * * * * * * * * * * * * * * * * * * * * * * * * * * !Use System_Constants ! For EXAMPLE !Use Elem_Type_Data ! for: LT_N, LT_PARM, COORD (LT_N, N_SPACE), & ! ! DLH (LT_PARM, LT_N), H (LT_N), etc !Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX ! IMPLICIT NONE ! INTEGER, INTENT(IN) :: IE ! ELEMENT NUMBER ! REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) ! CONSTITUTIVE MATRIX !! suppress compiler warnings (touch is never true) ! IF ( TOUCH ) PRINT *, IE, E (1, 1) ! STOP 'ERROR: NO SOURCE AT GET_CONSTANT_E_nnn' ! !b STOP 'Edit GET_CONSTANT_E to use my_e_matrix_inc' !END SUBROUTINE GET_CONSTANT_E_nnn ! ============= End Files for EXAMPLE number nnn ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 101 ============= SUBROUTINE DESCRIBE_EXAMPLE_101 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *,'*** DESCRIPTIONS OF EXAMPLE 101 ***' PRINT *,'Heat conduction through, convection from a bar:' PRINT *,'K*A*U,XX - h*P*(U-U_ext) = 0, U(0)=U_0, dU/dx(L)=0' PRINT *,'For globally constant data the analytic solution is:' PRINT *,'U(x) = U_ext - (U_0-U_ext) * Cosh [m*(L-x)] / Cosh [mL]' PRINT *,'Reaction is equal and opposite to the total convection ' PRINT *,'loss = h_e*P_e*(U_0-U_ext)*Sinh(m*L)/(m*Cosh_mL)' PRINT *,'m^2 = h_e*P_e/(K_e*A_e)' PRINT *,'U_ext = GET_REAL_MISC (1) for external temperature' PRINT *,'L = GET_REAL_MISC (2) for exact length' PRINT *,'U_0 = GET_REAL_MISC (3) for essential bc at x = 0' PRINT *,'K_e = GET_REAL_LP (1) thermal conductivity' PRINT *,'A_e = GET_REAL_LP (2) area of bar' PRINT *,'h_e = GET_REAL_LP (3) convection coefficent' PRINT *,'P_e = GET_REAL_LP (4) perimeter of area A_e' PRINT *,' ' PRINT *,'OR' PRINT *,'' PRINT *,'Structural pile (axial bar) in an elastic foundation:' PRINT *,'E*A*U,XX - k*P*U = 0, U(0)=U_0, dU/dx(L)=0' PRINT *,'For globally constant data the analytic solution is:' PRINT *,'U(x) = U_0 * Cosh [m*(L-x)] / Cosh [mL]' PRINT *,'Soil reaction is equal and opposite to total load:' PRINT *,'Soil = k_e*P_e*U_0*Sinh(m*L)/(m*Cosh_mL)' PRINT *,'m^2 = k_e*P_e/(E_e*A_e)' PRINT *,'Miscellaneous data used ONLY for analytic solution:' PRINT *,'U_ext = GET_REAL_MISC (1) external soil settlement, 0' PRINT *,'L = GET_REAL_MISC (2) for exact length' PRINT *,'U_0 = GET_REAL_MISC (3) for essential bc at x = 0' PRINT *,'Real FE problem properties are:' PRINT *,'E_e = GET_REAL_LP (1) elastic modulus, kN/cm^2' PRINT *,'A_e = GET_REAL_LP (2) area of bar, cm^2' PRINT *,'k_e = GET_REAL_LP (3) foundation stiffness, kN/cm^2' PRINT *,'P_e = GET_REAL_LP (4) perimeter of area A_e, cm' END SUBROUTINE DESCRIBE_EXAMPLE_101 SUBROUTINE ELEM_SQ_EX_101 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! ! Combined heat conduction through, convection from a bar: ! K*A*U,XX - h*P*(U-U_ext) = 0, U(0)=U_0, dU/dx(L)=0 ! For globally constant data the analytic soultion is: ! U(x) = U_ext - (U_0-U_ext) * Cosh [m*(L-x)] / Cosh [mL] ! Reaction is equal and opposite to the total convection ! loss = h_e*P_e*(U_0-U_ext)*Sinh(m*L)/(m*Cosh_mL) ! m^2 = h_e*P_e/(K_e*A_e) ! U_ext = GET_REAL_MISC (1) ! for external reference temperature ! L = GET_REAL_MISC (2) ! for exact length ! U_0 = GET_REAL_MISC (3) ! for essential bc at x = 0 ! K_e = GET_REAL_LP (1) ! thermal conductivity ! A_e = GET_REAL_LP (2) ! area of bar ! h_e = GET_REAL_LP (3) ! convection coefficent on perimeter ! P_e = GET_REAL_LP (4) ! perimeter of area A_e REAL(DP) :: DL, DX_DR ! Length, Jacobian REAL(DP) :: K_e, A_e, h_e, P_e, U_ext ! properties INTEGER :: IQ ! Loops IF ( DEBUG_PROPERTY ) PRINT *, 'Entering ELEM_SQ_EX_101 ' DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH DX_DR = DL / 2. ! CONSTANT JACOBIAN U_ext = GET_REAL_MISC (1) ! external temperature K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e E = A_e * K_e ! constitutive CALL STORE_FLUX_POINT_COUNT ! Save LT_QP !b ! S, C, H_INTG already zeroed DO IQ = 1, LT_QP ! LOOP OVER QUADRATURES ! GET INTERPOLATION FUNCTIONS, AND X-COORD H = GET_H_AT_QP (IQ) XYZ = MATMUL (H, COORD) ! ISOPARAMETRIC ! LOCAL AND GLOBAL DERIVATIVES DLH = GET_DLH_AT_QP (IQ) DGH = DLH / DX_DR ! CONVECTION SOURCE C = C + h_e * P_e * U_ext * H * WT (IQ) * DX_DR ! SQUARE MATRIX, CONDUCTION & CONVECTION S = S + ( K_e * A_e * MATMUL (TRANSPOSE(DGH), DGH) & + h_e * P_e * OUTER_PRODUCT (H, H) ) * WT (IQ) * DX_DR ! INTEGRATING FOR CONVECTION LOSS, POST PROCESSING H_INTG = H_INTG + h_e * P_e * H * WT (IQ) * DX_DR ! SAVE FOR FLUX AVERAGING OR POST PROCESSING, B == DGH CALL STORE_FLUX_POINT_DATA (XYZ, E, DGH) !b !b IF ( N_FILE1 > 0) WRITE (N_FILE1) DGH ! FOR GRADIENT POST_PROC END DO ! QUADRATURE IF ( N_FILE2 > 0) WRITE (N_FILE2) H_INTG ! FOR HEAT LOSS ! End of application dependent code ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_101 SUBROUTINE POST_PROCESS_ELEM_EX_101 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_101 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! ! Combined heat conduction through, convection from a bar: ! K*A*U,XX - h*P*(U-U_ext) = 0, U(0)=U_0, dU/dx(L)=0 ! For globally constant data the analytic soultion is: ! U(x) = U_ext - (U_0-U_ext) * Cosh [m*(L-x)] / Cosh [mL] ! Reaction is equal and opposite to the total convection ! loss = h_e*P_e*(U_0-U_ext)*Sinh(m*L)/(m*Cosh_mL) ! m^2 = h_e*P_e/(K_e*A_e) ! U_ext = GET_REAL_MISC (1) ! for external reference temperature ! L = GET_REAL_MISC (2) ! for exact length ! U_0 = GET_REAL_MISC (3) ! for essential bc at x = 0 ! K_e = GET_REAL_LP (1) ! thermal conductivity ! A_e = GET_REAL_LP (2) ! area of bar ! h_e = GET_REAL_LP (3) ! convection coefficent on perimeter ! P_e = GET_REAL_LP (4) ! perimeter of area A_e REAL(DP) :: EXACT, GRAD, VALUE, X_IQ, P_IQ (1) REAL(DP) :: DL, DX_DR REAL(DP) :: K_e, A_e, h_e, P_e ! properties REAL(DP), SAVE :: Q_E_LOSS, TOTAL_LOSS, EXACT_LOSS ! convection losses REAL(DP), SAVE :: L, m, COSH_ML, U_0, U_ext ! in exact sol INTEGER :: IQ, N_IP ! loops IF ( DEBUG_PROPERTY ) PRINT *, 'Entering POST_PROCESS_ELEM_EX_101 ' K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e IF ( IE == 1 ) THEN !--> WRITE GRADIENT HEADINGS WRITE (6, 10) 10 FORMAT ( /, '** ELEMENT GAUSS POINT RESULTS **',/, & & 'ELEMENT X EXACT FEA ', & & 'GRADIENT FE_GRADIENT' ) TOTAL_LOSS = 0.d0 ! INITIALIZE m = sqrt (h_e*P_e/(K_e*A_e)) ! for exact sol U_ext = GET_REAL_MISC (1) ! external ref temp L = GET_REAL_MISC (2) ! for exact U_0 = GET_REAL_MISC (3) ! essential bc at x = 0 COSH_ML = COSH (m*L) ! for exact END IF DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH DX_DR = DL / 2. ! CONSTANT JACOBIAN ! GET RESULTS AT FIRST NODE OF ELEMENT P_IQ (1) = -1.d0 ; X_IQ = COORD (1, 1) ! First node CALL SCALAR_SHAPES (P_IQ (1), H) ! shapes VALUE = DOT_PRODUCT ( H, D) ! fe value EXACT = U_ext + (U_0-U_ext)*COSH (m*(L-X_IQ)) / COSH_ML ! exact value CALL SCALAR_DERIVS (P_IQ (1), DLH) ! local grad DGH = DLH / DX_DR ! global gradient STRESS = MATMUL (DGH, D) ! FE GRADIENT GRAD = -(U_0-U_ext)*m*SINH (m*(L-X_IQ)) / COSH_ML ! exact grad WRITE (6, 11) IE, x_iq, exact, value, grad, STRESS (1) !b 11 FORMAT (I6, 2X, f6.3, 4(2x, 1PE12.5)) 11 FORMAT (I3, 1X, g9.4, 4(2x, 1PE12.5)) CALL READ_FLUX_POINT_COUNT (N_IP) !b !--> LOOP OVER QUADRATURE POINTS (ISOPARAMETRIC ELEMENT) !b DO IQ = 1, LT_QP DO IQ = 1, N_IP !b H = GET_H_AT_QP (IQ) ! shapes X_IQ = DOT_PRODUCT ( H, COORD (1:LT_N, 1)) ! x point VALUE = DOT_PRODUCT ( H, D) ! fe value EXACT = U_ext + (U_0-U_ext)*COSH(m*(L-X_IQ)) / COSH_ML ! exact value ! ---> CALCULATE GRADIENT, STRESS = DGH*D !b READ (N_FILE1) DGH ! Gradient matrix CALL READ_FLUX_POINT_DATA (XYZ, E, DGH) ! B=DGH ! Gradient matrix STRESS = MATMUL (DGH, D) ! FE GRADIENT GRAD = -(U_0-U_ext)*m*SINH (m*(L-X_IQ)) / COSH_ML ! exact grad WRITE (6, 11) IE, x_iq, exact, value, grad, STRESS (1) END DO ! GET RESULTS AT LAST NODE OF ELEMENT P_IQ (1) = 1.d0 ; X_IQ = COORD (LT_N, 1) ! Last node CALL SCALAR_SHAPES (P_IQ (1), H) ! shapes VALUE = DOT_PRODUCT ( H, D) ! fe value EXACT = U_ext + (U_0-U_ext)*COSH (m*(L-X_IQ)) / COSH_ML ! exact value CALL SCALAR_DERIVS (P_IQ (1), DLH) ! local grad DGH = DLH / DX_DR ! global gradient STRESS = MATMUL (DGH, D) ! FE GRADIENT GRAD = -(U_0-U_ext)*m*SINH (m*(L-X_IQ)) / COSH_ML ! exact grad WRITE (6, 11) IE, x_iq, exact, value, grad, STRESS (1) ! GET CONVECTIVE HEAT LOSS IF ( N_FILE2 > 0 ) THEN READ (N_FILE2) H_INTG Q_E_LOSS = DOT_PRODUCT ((D-U_ext), H_INTG) ! Elem loss WRITE (6, *) 'ELEMENT CONVECTION HEAT LOSS = ', Q_E_LOSS TOTAL_LOSS = TOTAL_LOSS + Q_E_LOSS IF ( IE == N_ELEMS ) THEN WRITE (6, *) 'TOTAL HEAT LOSS = ', TOTAL_LOSS EXACT_LOSS = H_e*P_e*(U_0-U_ext)*Sinh(m*L)/(m*Cosh_mL) WRITE (6, *) 'EXACT HEAT LOSS = ', EXACT_LOSS END IF ! LAST ELEMENT END IF ! LOSS DESIRED ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER END SUBROUTINE POST_PROCESS_ELEM_EX_101 ! =============== End Files for EXAMPLE number 101 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 102 ============= SUBROUTINE DESCRIBE_EXAMPLE_102 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 102 ***' PRINT *,'APPLICATION: LINEAR SLIDER BEARING' PRINT *,'[P,x h^3 / 6 nu],x = [U h],x + h,t' PRINT *,'UNKNOWN IS FLUID PRESSURE, P' PRINT *,'N_SPACE = 1, NOD_PER_EL = 2, N_G_DOF = 1' PRINT *,'N_EL_FRE = 2, MISC_FL = 2' PRINT *,'FLT_MISC(1) = VISCOSITY nu, FLT_MISC(2) = VELOCITY U' PRINT *,'EL_PROP(1) OR PRT_L_PT(K,1) = FILM THICKNESS, h' PRINT *,'EL_PROP(2) OR PRT_L_PT(K,2) = FILM SQUEEZE, h,t' END SUBROUTINE DESCRIBE_EXAMPLE_102 SUBROUTINE ELEM_SQ_EX_102 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! APPLICATION: LINEAR SLIDER BEARING, Example 102 ! Equation: [P,x h^3 / 6 nu],x = [U h],x + h,t ! N_SPACE = 1, NOD_PER_EL = 2, N_G_DOF = 1, N_EL_FRE = 2, ! MISC_FL = 2, either el_real or pt_real = 0, 1, or 2 ! GET_REAL_MISC (1) = VISCOSITY, GET_REAL_MISC (2) = VELOCITY ! GET_REAL_LP (1) OR PRT_L_PT (K,1) = FILM THICKNESS ! GET_REAL_LP (2) OR PRT_L_PT (K,2) = FILM SQUEEZE INTEGER, SAVE :: KALL = 1 REAL (DP), SAVE :: VIS, VEL, DL REAL (DP) :: THICK, CONST, SQUEEZE IF ( KALL == 1 ) THEN ! GET GLOBAL REAL CONSTANTS KALL = 0 VIS = GET_REAL_MISC (1) ; VEL = GET_REAL_MISC (2) END IF ! FIRST CALL !--> DEFINE ELEMENT LENGTH AND ELEMENT THICKNESS DL = COORD (2, 1) - COORD (1, 1) THICK = 0.d0 IF ( EL_REAL > 0 ) THICK = GET_REAL_LP (1) IF ( EL_REAL > 1 ) SQUEEZE = GET_REAL_LP (2) ! CHECK FOR ALTERNATE AVERAGE NODE THICKNESS IF ( THICK == 0.d0 ) THEN ! USE NODAL PROPERTY IF ( N_NP_FLO > 0 ) THEN ! DATA EXISTS THICK = 0.5d0 * (PRT_L_PT (1, 1) + PRT_L_PT (2, 1) ) SQUEEZE = 0.5d0 * (PRT_L_PT (1, 2) + PRT_L_PT (2, 2) ) ELSE STOP 'NO SLIDER BEARING THICKNESS DATA' END IF END IF ! NODAL THICKNESS DATA !--> GENERATE ELEMENT SQUARE MATRIX & COLUMN MATRIX CONST = THICK**3 / (6.0_DP * VIS * DL) S (1, 1) = CONST ; S (2, 2) = CONST S (1, 2) = -CONST ; S (2, 1) = -CONST C (1) = VEL * THICK + DL * SQUEEZE *0.5d0 !B check sign C (2) = -VEL * THICK + DL * SQUEEZE *0.5d0 !--> GENERATE DATA FOR LOAD CALCULATIONS AND STORE H_INTG (1) = 0.5_DP * DL ; H_INTG (2) = 0.5_DP * DL IF ( N_FILE1 > 0 ) WRITE (N_FILE1) H_INTG ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_102 SUBROUTINE POST_PROCESS_ELEM_EX_102 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_102 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! Define any new array or variable types, then give statements !APPLICATION: LINEAR SLIDER BEARING ! N_SPACE = 1, NOD_PER_EL = 2, N_G_DOF = 1 ! N_EL_FRE = 2, MISC_FL = 2 INTEGER, SAVE :: KALL = 1 REAL(DP), SAVE :: FORCE, TOTAL = 0.0_DP IF (KALL == 1) THEN !--> PRINT TITLES ON THE FIRST CALL KALL = 0 ; WRITE (6, 5) 5 FORMAT (/, '*** E L E M E N T L O A D S ***',/, & 'ELEMENT LOAD TOTAL') ENDIF !--> CALCULATE LOADS ON THE ELEMENTS, F = H_INTG*D READ (N_FILE1) H_INTG FORCE = DOT_PRODUCT (H_INTG, D) TOTAL = TOTAL + FORCE WRITE (6, 10) IE, FORCE, TOTAL 10 FORMAT (I5, 1PE16.5, 3X, 1PE16.5) ! *** END POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER END SUBROUTINE POST_PROCESS_ELEM_EX_102 ! ============= End Files for EXAMPLE number 102 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 103 ============= SUBROUTINE DESCRIBE_EXAMPLE_103 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 103 ***' PRINT *, 'Least squares solution of dy/dx + A y = F with y(0) = 0' PRINT *, 'Analytically integrated linear line element' PRINT *, 'A, F are misc real properties 1, 2 respectively' PRINT *, 'Exact solution: y(x) = (1 - e (-Ax)) * F / A' END SUBROUTINE DESCRIBE_EXAMPLE_103 SUBROUTINE ELEM_SQ_EX_103 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays !b REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & !b AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) !b REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & !b STRESS (N_R_B + 2) !b REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & !b DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! APPLICATION: LEAST SQUARES SOLUTION OF Y' + A * Y = F ! N_SPACE = 1, NOD_PER_EL = 2, N_G_DOF = 1 REAL(DP), SAVE :: A, F, DL ! DL = LENGTH OF ELEMENT ! RECOVER MISC PROBLEM COEFFICIENTS, A AND F (ON FIRST CALL) IF ( IE == 1 ) THEN ! Get coeffients A = GET_REAL_MISC (1) ; F = GET_REAL_MISC (2) END IF DL = COORD (2, 1) - COORD (1, 1) S (1, 1) = (3.d0 - 3.d0 * A * DL + A * A * DL * DL) / 3.d0 / DL S (2, 2) = (3.d0 + 3.d0 * A * DL + A * A * DL * DL) / 3.d0 / DL S (1, 2) = (A * A * DL * DL - 6.d0) / 6.d0 / DL S (2, 1) = S (1, 2) C (1) = 0.5d0 * F * (A * DL - 2.d0) C (2) = 0.5d0 * F * (A * DL + 2.d0) !b ! suppress compiler warnings (touch is never true) !b IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & !b DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & !b PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_103 SUBROUTINE POST_PROCESS_ELEM_EX_103 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_103 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_ELEM_EX_103' END SUBROUTINE POST_PROCESS_ELEM_EX_103 ! ============= End Files for EXAMPLE number 103 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 104 ============= SUBROUTINE DESCRIBE_EXAMPLE_104 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 104 ***' PRINT *, 'Solution of the ODE u" + u + x = 0 for x in ]0,1[ ' PRINT *, 'with various boundary conditions, by Galerkin method.' PRINT *, 'Optional source to post-process invoked if keywords "post_1"' PRINT *, 'and/or "post_2" appear in the keywords control for' PRINT *, 'gradient comparisons and solution norms' END SUBROUTINE DESCRIBE_EXAMPLE_104 SUBROUTINE ELEM_SQ_EX_104 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! NOTE: IS_ELEMENT SHOULD BE TRUE HERE ! .............................................................. ! Define any new array or variable types, then give statements ! ! APPLICATION DEPENDENT Galerkin MWR FOR ODE ! U,XX + U + X = 0, with boundary conditions like ! U(0)=0=U(1), so U = Sin(x)/Sin(1) - x or ! U(0)=0,U'(1)=0, so U = Sin(x)/Cos(1) - x etc REAL(DP) :: DL, DX_DR ! Length, Jacobian REAL(DP) :: NORM_SQ (LT_FREE, LT_FREE) ! elem norm INTEGER :: IQ ! Loops DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH IF ( DL < 0.d0 ) THEN PRINT *,'WARNING: CORRECTED NEGATIVE LENGTH ELEMENT ', THIS_EL N_WARN = N_WARN + 1 ; DL = ABS (DL) END IF DX_DR = DL / 2. ! CONSTANT JACOBIAN E = 1.d0 ! CONSTANT E NORM_SQ = 0.d0 ! S, C already zeroed IF ( U_FLUX > 0 ) WRITE (U_FLUX) LT_QP ! for SCP or post-process DO IQ = 1, LT_QP ! LOOP OVER QUADRATURES ! GET INTERPOLATION FUNCTIONS, AND X-COORD H = GET_H_AT_QP (IQ) XYZ = MATMUL (H, COORD) ! ISOPARAMETRIC ! LOCAL AND GLOBAL DERIVATIVES, B = DGH DLH = GET_DLH_AT_QP (IQ) DGH = DLH / DX_DR C = C + H * XYZ (1) * WT (IQ) * DX_DR ! SOURCE VECTOR ! SQUARE MATRIX S = S + ( MATMUL (TRANSPOSE(DGH), DGH) & - OUTER_PRODUCT (H, H) ) * WT (IQ) * DX_DR ! NORM MATRIX NORM_SQ = NORM_SQ + ( MATMUL (TRANSPOSE(DGH), DGH) & + OUTER_PRODUCT (H, H) ) * WT (IQ) * DX_DR IF ( U_FLUX > 0 ) WRITE (U_FLUX) XYZ, E, DGH ! for SCP or post-proc END DO ! QUADRATURE IF ( N_FILE2 > 0) WRITE (N_FILE2) NORM_SQ ! FOR H_1_NORM POST_PROC ! End of application dependent code ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_104 SUBROUTINE POST_PROCESS_ELEM_EX_104 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_104 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Define any new array or variable types, then give statements ! ! APPLICATION DEPENDENT Galerkin MWR FOR ODE ! U,XX + U + X = 0, U(0)=0=U(1), U = Sin(x)/Sin(1) - x REAL(DP) :: EXACT, GRAD, VALUE, X_IQ, P_IQ (1) REAL(DP) :: DL, DX_DR REAL(DP) :: NORM_SQ (LT_FREE, LT_FREE) ! elem norm matrix REAL(DP), SAVE :: H_1_ELEM, H_1_TOTAL = 0.d0 ! elem norms INTEGER :: IQ ! Do not duplicate exact flux results if given elsewhere !b IF ( N_FILE1 > 0 .AND. .NOT. USE_EXACT_FLUX ) THEN ! GET GRADIENTS IF ( N_FILE1 > 0 ) THEN ! GET GRADIENTS IF ( IE == 1 ) THEN !--> WRITE GRADIENT HEADINGS IF ( EXACT_CASE /= 9 .AND. EXACT_CASE /= 10 ) THEN PRINT *, 'WARNING: POST_PROCESS_ELEM_EX_104 "post_1" INVALID' PRINT *, 'EXCEPT FOR EXACT_CASE = 9 AND 10' N_WARN = N_WARN + 1 END IF WRITE (6, 10) 10 FORMAT ( /, '** ELEMENT FLUX RESULTS **',/, & & 'ELEMENT X EXACT FEA ', & & 'GRADIENT FE_GRADIENT' ) H_1_TOTAL = 0.d0 ! INITIALIZE END IF DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH DX_DR = DL / 2. ! CONSTANT JACOBIAN ! GET RESULTS AT FIRST NODE OF ELEMENT P_IQ (1) = -1.d0 ; X_IQ = COORD (1, 1) ! First node CALL SCALAR_SHAPES (P_IQ (1), H) ! shapes CALL SCALAR_DERIVS (P_IQ (1), DLH) ! local grad VALUE = DOT_PRODUCT ( H, D) ! fe value DGH = DLH / DX_DR ! physical gradient STRESS (1) = DOT_PRODUCT (DGH(1,1:LT_N), D(1:LT_N)) ! FE GRADIENT SELECT CASE (EXACT_CASE) CASE (9) exact = sin (x_iq) / sin (1.d0) - x_iq grad = cos (x_iq) / sin (1.d0) - 1.d0 CASE (10) exact = sin (x_iq) / cos (1.d0) - x_iq grad = cos (x_iq) / cos (1.d0) - 1.d0 CASE DEFAULT ; exact = 0.d0 ; grad = 0.d0 END SELECT WRITE (6, 11) IE, x_iq, exact, value, grad, STRESS (1) 11 FORMAT (I6, 2X, f6.3, 4(2x, 1PE12.5)) IF ( U_FLUX > 0 ) THEN READ (U_FLUX) LT_QP !--> LOOP OVER QUADRATURE POINTS (ISOPARAMETRIC ELEMENT) DO IQ = 1, LT_QP H = GET_H_AT_QP (IQ) ! shapes X_IQ = DOT_PRODUCT ( H, COORD (1:LT_N, 1)) ! x point VALUE = DOT_PRODUCT ( H, D) ! fe value ! ---> CALCULATE GRADIENT, STRESS = DGH*D READ (U_FLUX) XYZ, E, DGH STRESS (1) = DOT_PRODUCT (DGH(1,1:LT_N), D(1:LT_N)) ! FE GRAD SELECT CASE (EXACT_CASE) CASE (9) exact = sin (x_iq) / sin (1.d0) - x_iq grad = cos (x_iq) / sin (1.d0) - 1.d0 CASE (10) exact = sin (x_iq) / cos (1.d0) - x_iq grad = cos (x_iq) / cos (1.d0) - 1.d0 CASE DEFAULT ; exact = 0.d0 ; grad = 0.d0 END SELECT WRITE (6, 11) IE, x_iq, exact, value, grad, STRESS (1) END DO END IF ! saved qp ! GET RESULTS AT LAST NODE OF ELEMENT P_IQ (1) = 1.d0 ; X_IQ = COORD (LT_N, 1) ! Last node CALL SCALAR_SHAPES (P_IQ (1), H) ! shapes CALL SCALAR_DERIVS (P_IQ (1), DLH) ! local grad VALUE = DOT_PRODUCT ( H, D) ! fe value DGH = DLH / DX_DR ! physical gradient STRESS (1) = DOT_PRODUCT (DGH(1,1:LT_N), D(1:LT_N)) ! FE GRAD SELECT CASE (EXACT_CASE) CASE (9) exact = sin (x_iq) / sin (1.d0) - x_iq grad = cos (x_iq) / sin (1.d0) - 1.d0 CASE (10) exact = sin (x_iq) / cos (1.d0) - x_iq grad = cos (x_iq) / cos (1.d0) - 1.d0 CASE DEFAULT ; exact = 0.d0 ; grad = 0.d0 END SELECT WRITE (6, 11) IE, x_iq, exact, value, grad, STRESS (1) END IF ! gradients ! GET NORM DATA IF ( N_FILE2 > 0) THEN READ (N_FILE2) NORM_SQ ! FOR H_1_NORM POST_PROC H_1_ELEM = DOT_PRODUCT (D, (MATMUL(NORM_SQ, D))) ! Elem value WRITE (6, *) 'SQUARE OF ELEMENT H_1 NORM = ', H_1_ELEM H_1_TOTAL = H_1_TOTAL + H_1_ELEM IF ( IE == N_ELEMS ) WRITE (6, *) & 'SQUARE OF TOTAL H_1 NORM = ', H_1_TOTAL END IF ! NORMS DESIRED ! *** END POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER END SUBROUTINE POST_PROCESS_ELEM_EX_104 ! ============= End Files for EXAMPLE number 104 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 105 ============= SUBROUTINE DESCRIBE_EXAMPLE_105 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 105 ***' PRINT *, 'U" + Q(x) = 0, Usually with U(0)=0=U(1)' PRINT *, 'Q = X^L gives U(x) = (x-x^(L+2))/((L+1)(L+2))' PRINT *, 'L is first misc integer property. This is EXACT_CASE 11' PRINT *, 'If constant in elem Q is first real property.' PRINT *, 'Q_STEP = 1 for X<=1/2 else = 0 gives exact' PRINT *, 'U = X(3-4X)/8, X <= 1/2, else U = (1-X)/8, EXACT_CASE = 14' END SUBROUTINE DESCRIBE_EXAMPLE_105 SUBROUTINE ELEM_SQ_EX_105 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_SOURCE IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! NOTE: IS_ELEMENT SHOULD BE TRUE HERE ! .............................................................. ! Define any new array or variable types, then give statements ! ! APPLICATION DEPENDENT Galerkin MWR FOR ODE ! U,XX + Q = 0, U(0)=0=U(1) with Q in data file or ! Q(X) given by SELECT_EXACT_SOURCE (defaults to my_source_inc) REAL(DP) :: DL, DX_DR ! Length, Jacobian REAL(DP) :: SOURCE ! Q source term !b REAL(DP) :: NORM_SQ (LT_FREE, LT_FREE) ! elem norm INTEGER :: IQ, L ! Loops, optional DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH DX_DR = DL / 2. ! CONSTANT JACOBIAN ! SET CONSTANT PROPERTIES SOURCE = 0.d0 ; L = 0 ! DEFAULTS IF ( EL_REAL > 0 ) SOURCE = GET_REAL_LP (1) ! FOR Q CONST IF ( INTEGERS > 0 ) L = GET_INTEGER_MISC (1) ! FOR Q=X^L E(1, 1) = 1 ! identity matrix !b NORM_SQ = 0.d0 ! S, C already zeroed CALL STORE_FLUX_POINT_COUNT ! Save LT_QP !b DO IQ = 1, LT_QP ! LOOP OVER QUADRATURES ! GET INTERPOLATION FUNCTIONS, AND X-COORD H = GET_H_AT_QP (IQ) XYZ = MATMUL (H, COORD) ! ISOPARAMETRIC ! LOCAL AND GLOBAL DERIVATIVES, B = DGH DLH = GET_DLH_AT_QP (IQ) DGH = DLH / DX_DR ! Variable source Q(X) ? IF ( INTEGERS > 0 ) SOURCE = XYZ (1) **L ! Q = X^L ! Override via EXACT_CASE value or my_source_inc IF ( USE_EXACT_SOURCE ) CALL SELECT_EXACT_SOURCE (XYZ, SOURCE) ! RESULTANT SOURCE VECTOR C = C + SOURCE * H * WT (IQ) * DX_DR ! SQUARE MATRIX S = S + MATMUL (TRANSPOSE(DGH), DGH) * WT (IQ) * DX_DR !b ! NORM MATRIX !b NORM_SQ = NORM_SQ + ( MATMUL (TRANSPOSE(DGH), DGH) & !b + OUTER_PRODUCT (H, H) ) * WT (IQ) * DX_DR CALL STORE_FLUX_POINT_DATA (XYZ, E, DGH) !b END DO ! QUADRATURE !b IF ( N_FILE2 > 0) WRITE (N_FILE2) NORM_SQ ! FOR H_1_NORM POST_PROC ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_105 SUBROUTINE POST_PROCESS_ELEM_EX_105 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_105 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_ELEM_EX_105' END SUBROUTINE POST_PROCESS_ELEM_EX_105 ! ============= End Files for EXAMPLE number 105 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 106 ============= SUBROUTINE DESCRIBE_EXAMPLE_106 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 106 ***' PRINT *, "Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 482" PRINT *, "An unsymmetric Galerkin solution" PRINT *, "Y(0)=2, Y(1)=5/3, Y=X^4/6 - 3X^2/2 + X + 2, EXACT_CASE 15" PRINT *, "Y(0)=1, Y'(1)+Y(1)=0, Y=(X^4 -3X^2 -X +6)/6, EXACT_CASE 16" PRINT *, "Y'(0)+Y(0)=0, Y'(1)-Y(1)=3, Y=X^4/6 +3X^2/2 +X -1, EXACT_CASE 17" PRINT *, ' ' END SUBROUTINE DESCRIBE_EXAMPLE_106 SUBROUTINE ELEM_SQ_EX_106 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_SOURCE IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) REAL(DP) :: SOURCE ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! NOTE: IS_ELEMENT SHOULD BE TRUE HERE ! .............................................................. ! Define any new array or variable types, then give statements ! ! APPLICATION DEPENDENT Galerkin MWR FOR ODE ! Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1), Fausett p. 482 ! Y'' + f_1 Y' + f_2 Y = f_3 , Fausett p. 482 ! -f_3 given in SELECT_EXACT_SOURCE or my_source_inc REAL(DP) :: DL, DX_DR ! Length, Jacobian REAL(DP) :: f_1, f_2, f_3 ! Length, Jacobian !b REAL(DP) :: NORM_SQ (LT_FREE, LT_FREE) ! elem norm INTEGER :: IQ ! Loops DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH DX_DR = DL / 2. ! CONSTANT JACOBIAN E(1, 1) = 1 ! identity matrix !b NORM_SQ = 0.d0 ! S, C already zeroed CALL STORE_FLUX_POINT_COUNT ! Save LT_QP !b DO IQ = 1, LT_QP ! LOOP OVER QUADRATURES ! GET INTERPOLATION FUNCTIONS, AND X-COORD H = GET_H_AT_QP (IQ) XYZ = MATMUL (H, COORD) ! ISOPARAMETRIC ! GET VARIABLE COEFFICIENTS f_3 = 1.d0 + XYZ (1) **2 ; SOURCE = -f_3 f_2 = 2.d0 / f_3 f_1 = -XYZ (1) * f_2 ! LOCAL AND GLOBAL DERIVATIVES DLH = GET_DLH_AT_QP (IQ) DGH = DLH / DX_DR ! RESULTANT SOURCE VECTOR !b via EXACT_CASE value or my_source_inc !b CALL SELECT_EXACT_SOURCE (XYZ, SOURCE) C = C + SOURCE * H * WT (IQ) * DX_DR ! SQUARE MATRIX S = S + MATMUL (TRANSPOSE(DGH), DGH) * WT (IQ) * DX_DR & - f_1 * OUTER_PRODUCT (H, DGH(1, :)) * WT (IQ) * DX_DR & - f_2 * OUTER_PRODUCT (H, H) * WT (IQ) * DX_DR !b ! NORM MATRIX !b NORM_SQ = NORM_SQ + ( MATMUL (TRANSPOSE(DGH), DGH) & !b + OUTER_PRODUCT (H, H) ) * WT (IQ) * DX_DR CALL STORE_FLUX_POINT_DATA (XYZ, E, DGH) !b END DO ! QUADRATURE !b IF ( N_FILE2 > 0) WRITE (N_FILE2) NORM_SQ ! FOR H_1_NORM POST_PROC ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_106 SUBROUTINE MIXED_SQ_EX_106 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE MIXED_BC SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_ROBIN_DATA IMPLICIT NONE INTEGER, INTENT(IN) :: IE ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), H_INTG (LT_N), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** MIXED_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Mixed or Robin boundary condition, Standard form: ! K_n * U,n + ROBIN_1_SEG * U + ROBIN_2_SEG = 0 ! GET ROBIN_DATA COMPONENTS, These are point values !b XYZ (:) = COORD (1, :) ! A point "element" !b CALL SELECT_EXACT_ROBIN_DATA (XYZ, ROBIN_1_SEG, ROBIN_2_SEG) !b S (1, 1) = ROBIN_1_SEG !b sign ?? xxx !b C (1) = -ROBIN_2_SEG !b sign ?? xxx ! GET ROBIN_DATA COMPONENTS, These are point values IF ( MIXED_REAL > 1 ) THEN !bb ROBIN_1_SEG = GET_REAL_MX (1) ! recover mixed seg data ROBIN_2_SEG = GET_REAL_MX (2) ! recover mixed seg data ELSE !bb XYZ (:) = COORD (1, :) ! A point "element" !bb CALL SELECT_EXACT_ROBIN_DATA (XYZ, ROBIN_1_SEG, ROBIN_2_SEG) !bb END IF !bb S (1, 1) = ROBIN_1_SEG ! square "matrix" C (1) = -ROBIN_2_SEG ! column "vector" IF ( DEBUG_MIX_SQ ) THEN print *, 'MIXED_REAL ', MIXED_REAL print *, 'S (1, 1) ', S (1, 1) print *, 'C (1) ', C (1) END IF ! END APPLICATION DEPENDENT MIXED_SQ_MATRIX STATEMENTS ! suppress compiler warnings (never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE MIXED_SQ_EX_106 SUBROUTINE POST_PROCESS_ELEM_EX_106 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_106 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_ELEM_EX_106' END SUBROUTINE POST_PROCESS_ELEM_EX_106 ! ============= End Files for EXAMPLE number 106 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 107 ============= SUBROUTINE DESCRIBE_EXAMPLE_107 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 107 ***' PRINT *, 'LEAST SQUARES METHOD SOLUTION OF' PRINT *, "Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1)" PRINT *, 'USING 3-RD ORDER HERMITE ELEMENTS IN UNIT COORDINATES' END SUBROUTINE DESCRIBE_EXAMPLE_107 SUBROUTINE B_MATRIX_EX_107 (DGH, XYZ, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT B MATRIX FOR STRAIN OR GRADIENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for LT_N, LT_FREE, & H (LT_N) IF NEEDED IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP), INTENT(OUT) :: B (N_R_B, LT_FREE) ! ................................................................... ! ** GET_APPLICATION_B_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ** ! For REAL (DP) :: B (N_R_B, LT_FREE) ! Given REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N), XYZ (N_SPACE) ! ................................................................... ! LEAST SQ. SOL. OF: Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1) ! Define any new array or variable types, then give statements REAL(DP) :: DGV (LT_FREE) ! v' NOT GLOBAL as in my_el_sq_inc REAL(DP) :: D2GV (LT_FREE) ! v'' REAL(DP) :: PT_UNIT ! unit coord pt REAL(DP) :: DL ! physical length ! DGV = DERIVATIVE OF GENERALIZED (VECTOR) INTERPOLATION ! D2GV = SECOND DERIVATIVE OF GENERALIZED (VECTOR) INTERPOLATION ! used only in the error estimator IF ( DEBUG_B .OR. DEBUG_INCLUDE ) & WRITE (N_BUG, *) 'Entering my_b_matrix_inc' DL = COORD (2, 1) - COORD (1, 1) PT_UNIT = (XYZ (1) - COORD (1, 1)) / DL CALL DERIV_C1_L (PT_UNIT, DL, DGV) ! DV / DX CALL DERIV2_C1_L (PT_UNIT, DL, D2GV) ! D^2 V / DX^2 B (1, :) = DGV (:) ! NOT Global B (2, :) = D2GV (:) ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, DGH (1, 1), XYZ END SUBROUTINE B_MATRIX_EX_107 SUBROUTINE ELEM_SQ_EX_107 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Define any new array or variable types, then give statements ! LEAST SQ. SOL. OF: Y'' - 2XY'/(X^2+1) + 2Y/(X^2+1) = (X^2+1) ! a) WITH Y(0)=2, Y(1)=5/3, Y(X) = X^4/6 - 3X^2/2 + X + 2 (p482) ! b) WITH Y(0)=1 AND MPC Y'(1)+Y(1)=0 (p484) ! Y(X) = (x^4 -3X^2 -X +6)/6 ! c) WITH MPCs: Y'(0)+Y(0)=0, Y'(1)-Y(1)=3 ! SO Y(X) = X^4/6 + 3X^2/2 + X - 1, Fausett p. 485 ! USING 3-RD ORDER HERMITE ELEMENTS IN UNIT COORDINATES ! >>> NOTE NOT SCALAR INTERPOLATION, C_1 ELEMENTS <<< ! ------------------------------------------- REAL(DP) :: D2GV (LT_FREE) ! v'' REAL(DP) :: F (LT_FREE) ! Nonlinear Workspace REAL(DP) :: PT_UNIT, WT_UNIT ! unit coord pt, weight REAL(DP) :: DL, X_PT, X_SQ ! physical length INTEGER :: IP ! loops ! DGV = DERIVATIVE OF GENERALIZED (VECTOR) INTERPOLATION ! D2GV = SECOND DERIVATIVE OF GENERALIZED (VECTOR) INTERPOLATION ! V = GENERALIZED (VECTOR) INTERPOLATION FUNCTIONS IF ( DEBUG_EL_SQ .OR. DEBUG_INCLUDE ) & WRITE (N_BUG, *) 'Entering my_b_matrix_inc' E = GET_REAL_IDENTITY (N_R_B) ! DUMMY CONSTITUTIVE DL = COORD (2, 1) - COORD (1, 1) ! GET THE LENGTH CALL STORE_FLUX_POINT_COUNT ! Save LT_QP ! NUMERICAL INTEGRATION LOOP DO IP = 1, LT_QP !--> FIND UNIT COORDINATES AND WEIGHT FOR INTEGRATION PT_UNIT = (1.d0 + PT(1,IP)) / 2.d0 WT_UNIT = WT (IP) / 2.0d0 XYZ (1) = COORD (1, 1) + PT_UNIT * DL X_SQ = XYZ (1) * XYZ (1) !--> EVALUATE HERMITE SHAPE FUNCTIONS AND DERIVATIVES CALL SHAPE_C1_L (PT_UNIT, DL, V) ! V (NOT H) CALL DERIV_C1_L (PT_UNIT, DL, DGV) ! DV / DX CALL DERIV2_C1_L (PT_UNIT, DL, D2GV) ! D^2 V / DX^2 ! WORK VECTOR, F = Y'' - 2XY'/(X^2+1) + 2Y/(X^2-1) F = D2GV - 2.d0 * XYZ (1) * DGV(1,:) / (X_SQ + 1.d0) & + 2.d0 * V / (X_SQ + 1.d0) ! COMPLETE THE SQUARE MATRIX AND SOURCE VECTOR ! S_IJ = S_IJ + WT_UNIT * F_I * F_J * DL S = S + WT_UNIT * OUTER_PRODUCT (F, F) * DL C = C + WT_UNIT * (X_SQ + 1.d0) * F * DL ! STORE DATA FOR SCP OR POST PROCESSING B (1, :) = DGV (1, :) ; B (2, :) = D2GV (:) CALL STORE_FLUX_POINT_DATA (XYZ, E, B) END DO ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_107 SUBROUTINE POST_PROCESS_ELEM_EX_107 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_107 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_ELEM_EX_107' END SUBROUTINE POST_PROCESS_ELEM_EX_107 ! ============= End Files for EXAMPLE number 107 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 108 ============= SUBROUTINE DESCRIBE_EXAMPLE_108 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 108 ***' PRINT *, 'A conducting, convection truss member with heat generation' PRINT *, 'Equation: K*A U,xx - H*P (U - U_ext) - Q_e= 0,' PRINT *, 'U = temperature, K = conductivity,' PRINT *, 'A = area, H = convection coeff,' PRINT *, 'P = perimeter, Q_e = heat source per unit length.' END SUBROUTINE DESCRIBE_EXAMPLE_108 SUBROUTINE ELEM_SQ_EX_108 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Define any new array or variable types, then give statements ! 108.my_el_sq_inc ! A conducting, convection truss member with heat generation ! Equation: K*A U,xx - H*P (U - U_e) + Q_e = 0, ! U = temperature, K = conductivity, A = area, H = convection coeff, ! P = perimeter, Q_e = heat source per unit length ! 1 *---(K_e, A_e, h_e, P_e, U_e, Q_e)---* 2, Element in xyz REAL(DP) :: L_BAR ! Length REAL(DP) :: K_e, A_e, h_e, P_e, U_e, Q_e ! properties IF ( debug_el_sq .or. debug_include ) & WRITE (N_BUG, *) 'Entering my_el_sq_inc' L_BAR = SQRT( SUM( (COORD (2, 1:N_SPACE) & - COORD (1, 1:N_SPACE)) **2) ) K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e Q_e = GET_REAL_LP (5) ! source per unit length, BTU/ hr ft U_e = GET_REAL_LP (6) ! convecting temperature, F S (1, 1) = K_e * A_e / L_BAR + h_e * P_e * L_BAR / 3.d0 S (2, 1) = -K_e * A_e / L_BAR + h_e * P_e * L_BAR / 6.d0 S (1, 2) = -K_e * A_e / L_BAR + h_e * P_e * L_BAR / 6.d0 S (2, 2) = K_e * A_e / L_BAR + h_e * P_e * L_BAR / 3.d0 C (1) = (h_e * P_e * U_e + Q_e) * L_BAR / 2.d0 C (2) = (h_e * P_e * U_e + Q_e) * L_BAR / 2.d0 ! end file: my_el_sq_inc ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_108 ! ============= End Files for EXAMPLE number 108 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 109 ============= SUBROUTINE DESCRIBE_EXAMPLE_109 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 109 ***' PRINT *, 'Cylindrical Heat Transfer' PRINT *, '1/R * d[R K_RR dT/dR]/dR + Q = 0,' PRINT *, 'K_RR = GET_REAL_LP (1) ; Q = GET_REAL_LP (2' END SUBROUTINE DESCRIBE_EXAMPLE_109 SUBROUTINE E_MATRIX_EX_109 (IE, XYZ, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW ! (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT (IN) :: IE REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP), INTENT(OUT) :: E (N_R_B, N_R_B) ! ................................................................... ! ** GET_APPLICATION_E_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ** ! For required REAL (DP) :: E (N_R_B, N_R_B) ! Given INTEGER, INTENT (IN) :: N_R_B ! ................................................................... ! used only in the error estimator IF ( DEBUG_INCLUDE .OR. DEBUG_E ) & WRITE (N_BUG, *) 'Entering my_e_matrix_inc' E (1, 1) = GET_REAL_LP (1) ! K_RR ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, IE, XYZ END SUBROUTINE E_MATRIX_EX_109 SUBROUTINE ELEM_SQ_EX_109 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Define any new array or variable types, then give statements ! 109.my_el_sq_inc ! CYLINDRICAL HEAT CONDUCTION, (see Section 6.2) ! Global constant TWO_PI = 6.2831853072d0 REAL(DP) :: CONST, DET REAL(DP) :: K_RR, SOURCE INTEGER :: IP ! 1/R * d[R K_RR dT/dR]/dR + Q = 0, Example 109 ! PROP_1 = CONDUCTIVITY K_RR ! PROP_2 = SOURCE PER UNIT VOLUME, Q, or use my_exact_source_inc !--> DEFINE ELEMENT PROPERTIES K_RR = GET_REAL_LP (1) ; SOURCE = GET_REAL_LP (2) E (1, 1) = K_RR ! CONSTITUTIVE ! STORE NUMBER OF POINTS FOR FLUX CALCULATIONS CALL STORE_FLUX_POINT_COUNT ! Save LT_QP !--> NUMERICAL INTEGRATION LOOP DO IP = 1, LT_QP H = GET_H_AT_QP (IP) ! EVALUATE INTERPOLATION FUNCTIONS XYZ = MATMUL (H, COORD) ! FIND RADIUS (R), ISOPARAMETRIC DLH = GET_DLH_AT_QP (IP) ! FIND LOCAL DERIVATIVES AJ = MATMUL (DLH, COORD) ! FIND JACOBIAN AT THE PT ! FORM INVERSE AND DETERMINATE OF JACOBIAN CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) CONST = TWO_PI * DET * WT(IP) * XYZ (1) ! TWO_PI*|J|*w*R ! EVALUATE GLOBAL DERIVATIVES, DGH == B DGH = MATMUL (AJ_INV, DLH) B = COPY_DGH_INTO_B_MATRIX (DGH) ! B = DGH C = C + CONST * SOURCE * H ! VOLUMETRIC SOURCE OPTION ! CONDUCTION SQUARE MATRIX S = S + CONST * MATMUL ((MATMUL (TRANSPOSE (B), E)), B) !--> SAVE COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING CALL STORE_FLUX_POINT_DATA (XYZ, E, B) END DO ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_109 ! ============= End Files for EXAMPLE number 109 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 110 ============= SUBROUTINE DESCRIBE_EXAMPLE_110 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 110 ***' PRINT *, 'Cylindrical Stress Analysis, (see Section 6.3)' PRINT *, "Elem prop: Young's modulus, Poisson's ratio, Mass density" PRINT *, 'Misc prop: Spin about z-axis, radians per second' END SUBROUTINE DESCRIBE_EXAMPLE_110 SUBROUTINE B_MATRIX_EX_110 (DGH, XYZ, B) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT B MATRIX FOR STRAIN OR GRADIENTS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for LT_N, LT_FREE, & H (LT_N) IF NEEDED Use Interface_Header ! for COPY_DGH_INTO_B_MATRIX IMPLICIT NONE REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N) REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP), INTENT(OUT) :: B (N_R_B, LT_FREE) ! ................................................................... ! ** GET_APPLICATION_B_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ** ! For REAL (DP) :: B (N_R_B, LT_FREE) ! Given REAL (DP), INTENT(IN) :: DGH (N_SPACE, LT_N), XYZ (N_SPACE) ! ................................................................... B (1, :) = DGH (1, :) ! DU/DR radial strain B (2, :) = H (:) / XYZ (1) ! U/R hoop strain ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, DGH (1, 1), XYZ END SUBROUTINE B_MATRIX_EX_110 SUBROUTINE E_MATRIX_EX_110 (IE, XYZ, E) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET APPLICATION DEPENDENT E MATRIX FOR CONSTITUTIVE LAW ! (USED IF SUPERCONVERGENT PATCH GRADIENTS ARE ACTIVE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Sys_Properties_Data IMPLICIT NONE INTEGER, INTENT (IN) :: IE REAL (DP), INTENT(IN) :: XYZ (N_SPACE) REAL (DP), INTENT(OUT) :: E (N_R_B, N_R_B) ! ................................................................... ! ** GET_APPLICATION_E_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW ** ! For required REAL (DP) :: E (N_R_B, N_R_B) ! Given INTEGER, INTENT (IN) :: N_R_B ! ................................................................... REAL(DP) :: E_mod, P_ratio, Rho, Spin ! used only in the error estimator IF ( DEBUG_INCLUDE .OR. DEBUG_E ) & WRITE (N_BUG, *) 'Entering my_e_matrix_inc' E_mod = GET_REAL_LP (1) ; P_ratio = GET_REAL_LP (2) Rho = GET_REAL_LP (3) ; Spin = GET_REAL_MISC (1) ! CONSTITUTIVE LAW E(1,1) = E_mod * (1 - P_ratio)/((1 + P_ratio)*(1 - 2*P_ratio)) E(2,1) = E_mod * P_ratio/((1 + P_ratio)*(1 - 2*P_ratio)) E(1,2) = E(2,1) ; E(2,2) = E(1,1) ! suppress compiler warnings (touch is never true) IF ( TOUCH ) PRINT *, IE, XYZ END SUBROUTINE E_MATRIX_EX_110 SUBROUTINE ELEM_SQ_EX_110 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE, TWO_PI Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Define any new array or variable types, then give statements ! 110.my_el_sq_inc ! Cylindrical Stress Analysis, (see Section 6.3) REAL(DP) :: CONST, DET REAL(DP) :: E_mod, P_ratio, Rho, Spin INTEGER :: IP ! Elem real prop: Young's modulus, Poisson's ratio, Mass density ! Misc real prop: Spin about z-axis, radians per second E_mod = GET_REAL_LP (1) ; P_ratio = GET_REAL_LP (2) Rho = GET_REAL_LP (3) ; Spin = GET_REAL_MISC (1) ! CONSTITUTIVE LAW E(1,1) = E_mod * (1 - P_ratio)/((1 + P_ratio)*(1 - 2*P_ratio)) E(2,1) = E_mod * P_ratio/((1 + P_ratio)*(1 - 2*P_ratio)) E(1,2) = E(2,1) ; E(2,2) = E(1,1) ! STORE NUMBER OF POINTS FOR FLUX CALCULATIONS CALL STORE_FLUX_POINT_COUNT ! Save LT_QP !--> NUMERICAL INTEGRATION LOOP DO IP = 1, LT_QP H = GET_H_AT_QP (IP) ! EVALUATE INTERPOLATION FUNCTIONS XYZ = MATMUL (H, COORD) ! FIND RADIUS (R), ISOPARAMETRIC DLH = GET_DLH_AT_QP (IP) ! FIND LOCAL DERIVATIVES AJ = MATMUL (DLH, COORD) ! FIND JACOBIAN AT THE PT ! FORM INVERSE AND DETERMINATE OF JACOBIAN CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) CONST = TWO_PI * DET * WT(IP) * XYZ (1) ! TWO_PI*|J|*w*R ! EVALUATE GLOBAL DERIVATIVES & STRAIN-DISPLACEMENT DGH = MATMUL (AJ_INV, DLH) B (1, :) = DGH (1, :) ! DU/DR radial strain B (2, :) = H (:) / XYZ (1) ! U/R hoop strain ! STIFFNESS MATRIX S = S + CONST * MATMUL ((MATMUL (TRANSPOSE (B), E)), B) BODY = - Rho * XYZ (1) * Spin **2 ! -Rho R Omega^2 C = C + CONST * BODY (1) * H ! CENTRIFUGAL FORCE RESULT !--> SAVE COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING CALL STORE_FLUX_POINT_DATA (XYZ, E, B) END DO ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_110 ! ============= End Files for EXAMPLE number 110 ============= ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 111 ============= SUBROUTINE DESCRIBE_EXAMPLE_111 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 111 ***' PRINT *, 'Direct Current Network' PRINT *, 'Elem prop: resistance' PRINT *, 'Results: Voltage at nodes' PRINT *, 'ELEMENT REACTIONS ARE THE IN & OUT CURRENT FLOWS' END SUBROUTINE DESCRIBE_EXAMPLE_111 SUBROUTINE ELEM_SQ_EX_111 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE, TWO_PI Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! Define any new array or variable types, then give statements ! 111.my_el_sq_inc ! NOTE: ELEMENT REACTIONS ARE THE IN & OUT CURRENT FLOWS REAL(DP) :: resistance, conductance ! R, 1/R ! Get resistance resistance = GET_REAL_LP (1) ! real element property IF ( resistance > 0.d0 ) THEN conductance = 1.d0 / resistance ELSE PRINT *, 'WARNING: Invalid resistance, element ', IE N_WARN = N_WARN + 1 ; conductance = 1.d0 END IF ! valid data ! Conductance matrix S (1,1) = conductance ; S (2,1) = -conductance S (1,2) = -conductance ; S (2,2) = conductance ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_111 ! ============= End Files for EXAMPLE number 111 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 112 ============= SUBROUTINE DESCRIBE_EXAMPLE_112 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 112 ***' PRINT *, 'SUPG Advection-Diffusion in 1-D for L2, L3, or L4' PRINT *, 'u * dp/dx - d(k * dp/dx)/dx = Q; u, k, Q constan' PRINT *, 'u = GET_REAL_LP (1) ! velocity' PRINT *, 'k = GET_REAL_LP (2) ! conductivity' PRINT *, 'Q = GET_REAL_LP (3) ! source per unit length' PRINT *, 'Misc integer: Galerkin option = 0, 1=SUPG' PRINT *, 'Misc real: fake Pe overwritten for exact solution' END SUBROUTINE DESCRIBE_EXAMPLE_112 SUBROUTINE ELEM_SQ_EX_112 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE, X) !b L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE) ! Global coordinates ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! ..................................................... ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! ..................................................... ! Define any new array or variable types, then give statements ! ! Galerkin or SUPG 1-d Advection-Diffusion Problem ! u * dp/dx - d(k * dp/dx)/dx = Q, assume u, k, Q constant ! ! u = GET_REAL_LP (1) ! velocity ! k = GET_REAL_LP (2) ! conductivity ! Q = GET_REAL_LP (3) ! source per unit length ! OPTION = GET_INTEGER_MISC (1) ! 0=Galerkin (default), 1=SUPG ! Pe_sys = GET_REAL_MISC (1) ! System Peclet, for exact solution ! LT_N = number of nodes for this element type ! MISC_FL = number of real miscellaneous properties ! MISC_FX = number of integer miscellaneous properties ! N_LP_FLO = number of real element properties ! W = Petrov weight, DGW its global derivative REAL(DP) :: W (LT_N), DGW (1, LT_N) ! SUPG & deriv REAL(DP) :: D2GH (1, LT_N) ! SUPG, zero ? REAL(DP) :: DL, DX_DR, DL_A ! Length, Jacobian REAL(DP) :: u, k, Q ! speed, diffusion, source REAL(DP) :: Pe, A, COTH ! Peclet data, L2 INTEGER :: IQ ! Loops REAL(DP), SAVE :: Pe_max ! debugging INTEGER, SAVE :: OPTION = 0 ! Method of solution DL = COORD (LT_N, 1) - COORD (1, 1) ! LENGTH ELEM DX_DR = DL / 2. ! CONSTANT JACOBIAN DL_A = DL / (LT_N - 1) ! SUPG length ! DATA READS AND SAVES u = GET_REAL_LP (1) ! velocity k = GET_REAL_LP (2) ; E = k ! conductivity Q = 0.d0 ; IF ( N_LP_FLO > 2 ) Q = GET_REAL_LP (3) ! source E = k ! constitutive IF ( IE == 1 ) THEN ! FIRST ELEMENT, ONE TIME ACTIONS Pe_max = 0.d0 ! initialize IF ( MISC_FX > 0 ) OPTION = GET_INTEGER_MISC (1) ! 0=Galerkin, 1=SUPG IF ( OPTION == 0 ) THEN ! echo choice PRINT *,'NOTE: Galerkin method' ! default ELSE ; PRINT *,'NOTE: SUPG method' ; END IF ! supg IF ( (USE_EXACT .OR. USE_EXACT_FLUX) .AND. MISC_FL == 0 ) THEN PRINT *, 'WARNING: ELEM_SQ_MATRIX, to use Pe_sys in exact ' PRINT *, 'solution set keyword "reals 1" in control' N_WARN = N_WARN + 1 ; END IF ! exact data given END IF ! FIRST ELEMENT ! SUPG TERMS (ASSUMING L2 ELEMENT), ? Bias for L3, L4 ? Pe = u * DL / k ! Grid Peclet IF ( Pe > Pe_max ) Pe_max = Pe ! for debug COTH = COSH (Pe/2) / SINH (Pe/2) ! Optimal SUPG A = ABS (COTH) - 1.d0 / ABS (Pe/2) ! Optimal SUPG L2 A = SIGN (A, u) ! abs(A) * sign of u ! ------- ELEMENT MATRICES FORMATION ---------- CALL STORE_FLUX_POINT_COUNT ! Save LT_QP FOR SCP DO IQ = 1, LT_QP ! LOOP OVER QUADRATURES, S, C already zeroed ! GET TRIAL INTERPOLATION FUNCTIONS, AND X-COORD H = GET_H_AT_QP (IQ) ! SOLUTION INTERPOLATION XYZ = MATMUL (H, COORD) ! ISOPARAMETRIC ! LOCAL AND GLOBAL FIRST DERIVATIVES DLH = GET_DLH_AT_QP (IQ) ! LOCAL DERIVATIVE DGH = DLH / DX_DR ! PHYSICAL DERIVATIVE ! *** SELECT STANDARD GALERKIN OR SUPG *** IF ( OPTION == 0 ) THEN ! Galerkin W = H ; DGW (1, :) = DGH (1, :) ELSE ! SUPG Method ! LOCAL AND GLOBAL SECOND DERIVATIVES (FOR N_SPACE == 1) SELECT CASE (LT_N) ! ELEMENT LIBRARY CASE (2) ; D2LH = 0.d0 ! L2 CASE (3) ; CALL DERIV2_3_L (PT (1, IQ), D2LH (1, :)) ! L3 CASE (4) ; CALL DERIV2_4_L (PT (1, IQ), D2LH (1, :)) ! L4 CASE DEFAULT ; STOP 'SUPG, NO SECOND DERIVATIVE IN LIBRARY' END SELECT D2GH = D2LH / DX_DR**2 ! PHYSICAL SECOND DERIVATIVE ! SUPG WEIGHTINGS, NOTE SECOND DERIVATIVE IN DGW W = H + A * DGH (1, :) * DL_A * 0.5d0 DGW (1, :) = DGH (1, :) + A * D2GH (1, :) * DL_A * 0.5d0 ! PRE-INSERT SECOND DERIVATIVE RESIDUAL, IF ANY IF ( LT_N > 2 ) S = S + k * A * DL_A * WT (IQ) * DX_DR & * OUTER_PRODUCT (DGH (1, :), D2GH (1, :)) END IF ! Method option ! MATRICES: SOURCE, CONDUCTION & ADVECTION C = C + Q * W * WT (IQ) * DX_DR ! SOURCE S = S + ( k * MATMUL (TRANSPOSE(DGH), DGH) & + u * OUTER_PRODUCT (W, DGH(1,:) )) * WT (IQ) * DX_DR !--> SAVE COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING CALL STORE_FLUX_POINT_DATA (XYZ, E, DGH) END DO ! QUADRATURE IF ( IE == N_ELEMS ) PRINT *, 'Maximum element Pe = ', PE_max ! *** END ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS *** END SUBROUTINE ELEM_SQ_EX_112 SUBROUTINE POST_PROCESS_ELEM_EX_112 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_112 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER STOP 'ERROR: NO SOURCE AT POST_PROCESS_ELEM_EX_112' END SUBROUTINE POST_PROCESS_ELEM_EX_112 ! ============= End Files for EXAMPLE number 112 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 113 ============= SUBROUTINE DESCRIBE_EXAMPLE_113 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 113 ***' PRINT *,'An axial elastic bar using closed form linear elements' PRINT *,'Element properties: 1-area, 2-elastic modulus, 3-temperature' PRINT *,'rise, 4-coeff. of thermal expansion, 5- weight per unit volume' PRINT *,'Element post-processing lists the mechanical strain, thermal' PRINT *,'strain, and mechanical stress, at the element center' END SUBROUTINE DESCRIBE_EXAMPLE_113 SUBROUTINE ELEM_SQ_EX_113 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_* IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .............................................................. ! AN AXIAL BAR BY DIRECT ENERGY APPROACH ! ELEMENT REAL PROPERTIES: (1) = AREA, (2) = ELASTIC MODULUS ! (3) = TEMP RISE, (4) = COEFF EXPANSION, (5) = WEIGHT DENSITY REAL(DP) :: BAR_L ! length REAL(DP) :: DELTA_T, ALPHA ! temp rise, expansion REAL(DP) :: AREA, GAMMA ! area, wt. density REAL(DP) :: M_E, THERMAL ! modulus, thermal strain ! Get properties for this element AREA = GET_REAL_LP (1); M_E = GET_REAL_LP (2) DELTA_T = GET_REAL_LP (3); ALPHA = GET_REAL_LP (4) GAMMA = GET_REAL_LP (5) ! Find bar length and direction cosines BAR_L = COORD (2, 1) - COORD (1, 1) ! length ! Form global strain-displacement matrix B (1, :) = (/ -1, 1 /) / BAR_L ! Form global stiffness, S = B' EAL B S = M_E * AREA * BAR_L * MATMUL ( TRANSPOSE (B), B ) ! Initial (thermal) strain loading THERMAL = ALPHA * DELTA_T ! strain C = B (1, :) * M_E * THERMAL * AREA * BAR_L ! force ! Weight load, in positive X-direction (wt density * volume) C = C + (/ 0.5d0, 0.5d0 /) * GAMMA * AREA * BAR_L ! total weight ! Save for stress post-processing (set post_1 in keywords) IF ( N_FILE1 > 0 ) WRITE (N_FILE1) M_E, B, THERMAL ! End of application dependent code ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_113 SUBROUTINE POST_PROCESS_ELEM_EX_113 (PRT_L_PT, PRT_MAT, & L_PT_PROP, ITER, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! ELEMENT LEVEL POST-SOLUTION CALCULATIONS, FROM ELEM_POST_EX_113 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions IMPLICIT NONE INTEGER, INTENT(IN) :: ITER, IE ! ITERATION, ELEMENT NUMBER ! OPTIONAL PROPERTY AND SOLUTION VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays (Usually used.) REAL(DP) :: B (N_R_B, LT_FREE), DGH (N_SPACE, LT_N), & DGV (N_SPACE, LT_FREE), & E (N_R_B, N_R_B), EB (N_R_B, LT_FREE), & STRAIN (N_R_B + 2), STRAIN_0 (N_R_B), & STRESS (N_R_B + 2), BODY (N_SPACE) ! OPTIONAL FOR NUMERICAL INTEGRATION REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), & H_INTG (LT_N), XYZ (N_SPACE) ! VARIABLES: See file NOTATION.f ! .................................................... ! *** POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! .................................................... ! AN AXIAL BAR BY DIRECT ENERGY APPROACH ! ELEMENT REAL PROPERTIES: (1) = AREA, (2) = ELASTIC MODULUS ! (3) = TEMP RISE, (4) = COEFF EXPANSION, (5) = WEIGHT DENSITY ! STRESS = M_E * (MECHANICAL STRAIN - INITIAL STRAIN) REAL(DP) :: THERMAL, M_E ! initial strain, modulus LOGICAL, SAVE :: FIRST = .TRUE. ! printing IF ( FIRST ) THEN ! first call FIRST = .FALSE. ; WRITE (6, 5) ! print headings 5 FORMAT (' E L E M E N T S T R E S S E S', /, & & ' ELEMENT STRESS MECH. STRAIN THERMAL STRAIN') END IF ! first call !--> Read stress strain data from N_FILE1 (set by post_1) READ (N_FILE1) M_E, B, STRAIN_0 (1) ! THERMAL = STRAIN_0 !--> Calculate mechanical strain, STRAIN = B * D STRAIN (1) = DOT_PRODUCT ( B(1, :), D ) !--> Generalized Hooke's Law STRESS (1) = M_E * (STRAIN (1) - STRAIN_0 (1)) WRITE (6, 1) IE, STRESS (1), STRAIN (1), STRAIN_0 (1) 1 FORMAT (I5, 3ES15.5) ! *** END POST_PROCESS_ELEM PROBLEM DEPENDENT STATEMENTS *** ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) IF ( TOUCH ) PRINT *, ITER END SUBROUTINE POST_PROCESS_ELEM_EX_113 !============= End Files for EXAMPLE number 113 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 114 ============= SUBROUTINE DESCRIBE_EXAMPLE_114 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE PRINT *, '*** DESCRIPTIONS OF EXAMPLE 114 ***' PRINT *,'Heat & Mass Transfer in 1-D: (unsymmetric system) L2 ONLY' PRINT *,'Equation: K*A U,xx - H*P (U - U_e) - m_dot*c_p*U,x + Q_e = 0,' PRINT *,'1 *---(Properties)---* 2, Properties are:' PRINT *,'1: K_e = thermal conductivity' PRINT *,'2: A_e = area of bar' PRINT *,'3: h_e = convection coefficent on perimeter' PRINT *,'4: P_e = perimeter of area A_e' PRINT *,'5: Q_e = source per unit length, BTU/ hr ft' PRINT *,'6: U_e = convecting temperature, F' PRINT *,'7: m_dot = mass flow rate' PRINT *,'8: c_p = specific heat, constant pressure' END SUBROUTINE DESCRIBE_EXAMPLE_114 SUBROUTINE ELEM_SQ_EX_114 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_* IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! For required REAL (DP) :: S (LT_FREE, LT_FREE) ! and optional REAL (DP) :: C (LT_FREE) ! .............................................................. ! Define any new array or variable types, then give statements ! Heat and Mass Transfer in 1-D: (unsymmetric system) L2 ONLY ! Equation: K*A U,xx - H*P (U - U_e) - m_dot*c_p*U,x + Q_e = 0, ! U = temperature, K = conductivity, A = area, H = convection coeff, ! P = perimeter, m_dot = mass flow rate, c_p = specific heat, ! Q_e = heat source per unit length ! 1 *---(K_e, A_e, h_e, P_e, U_e, m_dot, c_p, Q_e)---* 2, Element in x REAL(DP) :: L_BAR ! Length REAL(DP) :: K_e, A_e, h_e, P_e, U_e, Q_e, m_dot, c_p ! properties L_BAR = SQRT( SUM( (COORD (2, 1:N_SPACE) & - COORD (1, 1:N_SPACE)) **2) ) IF ( EL_REAL < 8 ) PRINT *, 'STOP: NEED el_real 8 FOR MASS TRANSFER' K_e = GET_REAL_LP (1) ! thermal conductivity A_e = GET_REAL_LP (2) ! area of bar h_e = GET_REAL_LP (3) ! convection coefficent on perimeter P_e = GET_REAL_LP (4) ! perimeter of area A_e Q_e = GET_REAL_LP (5) ! source per unit length, BTU/ hr ft U_e = GET_REAL_LP (6) ! convecting temperature, F m_dot = GET_REAL_LP (7) ! mass flow rate c_p = GET_REAL_LP (8) ! specific heat, constant pressure ! conduction and convection effects S (1, 1) = K_e * A_e / L_BAR + h_e * P_e * L_BAR / 3.d0 S (2, 1) = -K_e * A_e / L_BAR + h_e * P_e * L_BAR / 6.d0 S (1, 2) = -K_e * A_e / L_BAR + h_e * P_e * L_BAR / 6.d0 S (2, 2) = K_e * A_e / L_BAR + h_e * P_e * L_BAR / 3.d0 ! mass transfer effects S (1, 1) = S (1, 1) - m_dot * c_p / 2 S (2, 1) = S (2, 1) - m_dot * c_p / 2 S (1, 2) = S (1, 2) + m_dot * c_p / 2 S (2, 2) = S (2, 2) + m_dot * c_p / 2 ! convection and source effects C (1) = (h_e * P_e * U_e + Q_e) * L_BAR / 2.d0 C (2) = (h_e * P_e * U_e + Q_e) * L_BAR / 2.d0 ! end file: my_el_sq_inc ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_114 !============= End Files for EXAMPLE number 114 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 115 ============= !============= End Files for EXAMPLE number 115 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 116 ============= !============= End Files for EXAMPLE number 116 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 117 ============= !============= End Files for EXAMPLE number 117 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 118 ============= !============= End Files for EXAMPLE number 118 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 119 ============= !============= End Files for EXAMPLE number 119 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 120 ============= !============= End Files for EXAMPLE number 120 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 121 ============= SUBROUTINE DESCRIBE_EXAMPLE_121 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE = 121 PRINT *, "*** DESCRIPTIONS OF EXAMPLE ", EXAMPLE, " ***" PRINT *, "Galerkin FOR TRANSIENT DE VIA ITERATION " PRINT *, "K U,xx - R U,t = 0, with U(x,0), U(0,t), U(L,t)" PRINT *, "----------------------------------------" PRINT *, "Uses iteration loops to manually treat time integration" PRINT *, "NOTE: This is very inefficient since the matrices" PRINT *, "generation, assembly, boundary conditions, and " PRINT *, "factorization are repeated ever time step (iteration)" PRINT *, "Need key relaxation 1.0 so D_OLD replaced by D_NEW" PRINT *, "Keyword time_step 1.d-3, etc sets time step size" PRINT *, "Keyword time_method sets method: " PRINT *, " 1-Euler, 2-Crank-Nicolson, 3-Galerkin, 4-Least Sq" END SUBROUTINE DESCRIBE_EXAMPLE_121 SUBROUTINE ELEM_SQ_EX_121 (E, H_INTG, PRT_L_PT, PRT_MAT, & L_PT_PROP, IE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GENERATE ELEMENT SQUARE MATRIX, OPTIONAL COLUMN MATRIX ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE Use Geometric_Properties Use Elem_Type_Data ! for: Use Sys_Properties_Data ! for FLT_MISC, MISC_FIX Use Interface_Header ! for functions Use Select_Source ! for SELECT_EXACT_* IMPLICIT NONE INTEGER, INTENT(IN) :: IE REAL(DP), INTENT(INOUT) :: E (N_R_B, N_R_B) REAL(DP), INTENT(INOUT) :: H_INTG (LT_N) !b make optional !uu ! OPTIONAL PROPERTY VALUES REAL(DP), INTENT(IN) :: PRT_L_PT (LT_N, N_NP_FLO), & PRT_MAT (MISC_FL) INTEGER, INTENT(IN) :: L_PT_PROP (N_NP_FIX) ! Automatic Arrays REAL(DP) :: AJ (N_SPACE, N_SPACE), XYZ (N_SPACE), & AJ_INV (N_SPACE, N_SPACE), BODY (N_SPACE) REAL(DP) :: STRAIN_0 (N_R_B), STRAIN (N_R_B + 2), & STRESS (N_R_B + 2) REAL(DP) :: B (N_R_B, LT_FREE), EB (N_R_B, LT_FREE), & DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE) ! VARIABLES: See file NOTATION.f ! .............................................................. ! begin file 121.my_el_sq_inc ! *** ELEM_SQ_MATRIX PROBLEM DEPENDENT STATEMENTS FOLLOW *** ! For required REAL (DP) :: S (LT_FREE, LT_FREE) ! and optional REAL (DP) :: C (LT_FREE) ! using previous solution D (LT_FREE) with ! D = D_OLD + RELAXATION (D_NEW - D_OLD) via key relaxation 1.0 ! Globals START_TIME, TIME_STEP, TIME are available ! DIAG_M (LT_FREE) is global ! .............................................................. ! APPLICATION DEPENDENT Galerkin FOR TRANSIENT DE ! via element level assembly of semi-discrete form(s) ! K_e U,xx + Q - Rho_e U,t = 0, with U(x,0) given by keyword ! start_value, U(0,t) and U(L,t) given by EBC constant in time ! Keywords time_step, time_method, scalar_source, ! diagonal_mass are also available ! Spatial items REAL(DP) :: DL, DX_DR, Q = 0.d0 ! Length, Jacobian, source INTEGER :: IQ, J ! Loops ! Work items for time integration INTEGER, SAVE :: METHOD = 1 ! defaults REAL(DP), SAVE :: Del_t = 1.d0 ! defaults REAL(DP), SAVE :: K_e = 1.d0, Rho_e = 1.d0 ! defaults REAL(DP) :: K (LT_FREE, LT_FREE), M (LT_FREE, LT_FREE) REAL(DP) :: F (LT_FREE), WORK (LT_FREE, LT_FREE) IF ( debug_el_sq .or. debug_include ) & WRITE (N_BUG, *) 'Entering my_el_sq_inc' IF ( THIS_EL == 1 ) THEN METHOD = TIME_METHOD ! keyword option time_method IF ( INTEGERS > 0 ) METHOD = GET_INTEGER_MISC (1) ! alt IF ( TIME_STEP /= 1.d0 ) Del_t = TIME_STEP TIME = START_TIME + THIS_ITER * TIME_STEP PRINT *, 'TIME = ', TIME END IF K = 0 ; M = 0 ; F = 0 ! Initialize spatial arrays IF ( N_LP_FLO > 1 ) THEN ! use non-default properties K_e = GET_REAL_LP (1) ! thermal conductivity Rho_e = GET_REAL_LP (2) ! rho, c_p END IF E = K_e ! Diffusion DL = COORD (LT_N, 1) - COORD (1, 1) ! Length DX_DR = DL / 2. ! Jacobian IF ( SCALAR_SOURCE /= 0.d0 ) Q = SCALAR_SOURCE ! source CALL STORE_FLUX_POINT_COUNT ! Save LT_QP for post-processing DO IQ = 1, LT_QP ! LOOP OVER QUADRATURES ! GET INTERPOLATION FUNCTIONS, AND X-COORD H = GET_H_AT_QP (IQ) XYZ = MATMUL (H, COORD) ! ISOPARAMETRIC ! LOCAL AND GLOBAL DERIVATIVES, B = DGH DLH = GET_DLH_AT_QP (IQ) ! dH / dr DGH = DLH / DX_DR ! dH / dx ! SQUARE MATRICES (STIFFNESS & MASS) & SOURCE VECTOR K = K + MATMUL (TRANSPOSE(DGH), DGH)* WT (IQ) * DX_DR * K_e M = M + OUTER_PRODUCT (H, H) * WT (IQ) * DX_DR * Rho_e F = F + H * Q * WT (IQ) * DX_DR CALL STORE_FLUX_POINT_DATA (XYZ, E, DGH) ! for post-processing END DO ! QUADRATURE ! Assemble for semi-discrete one-step time rule ! NOTE: This is very inefficient since the assembly, boundary ! conditions, and factorization are repeated ever time step IF ( DIAGONAL_MASS ) THEN ! use the scaled diagonal mass form CALL DIAGONALIZE_SQ_MATRIX (LT_FREE, M, DIAG_M) ; M = 0.d0 DO J = 1, LT_FREE M (J, J) = DIAG_M (J) END DO END IF SELECT CASE ( METHOD ) ! for one-step time integration rule CASE ( 2 ) ! Crank-Nicolson WORK = M / Del_t - K / 2 C = F + MATMUL (WORK, D) S = M / Del_t + K / 2 CASE ( 3 ) ! Galerkin continuous in time WORK = M / Del_t - K / 3 C = F + MATMUL (WORK, D) S = 2 * K / 3.d0 + M / Del_t CASE ( 4 ) ! Least Squares in time, F constant in time WORK = MATMUL (TRANSPOSE(M), M) / Del_t & + MATMUL (TRANSPOSE(K), M) / 2 & - MATMUL (TRANSPOSE(M), K) / 2 & - MATMUL (TRANSPOSE(K), K) * Del_t / 6 C = -MATMUL (TRANSPOSE(K), F) * Del_t / 2 & - MATMUL (TRANSPOSE(M), F) & + MATMUL (WORK, D) S = MATMUL (TRANSPOSE(K), K) * Del_t / 3 & + MATMUL (TRANSPOSE(K), M) / 2 & + MATMUL (TRANSPOSE(M), K) / 2 & + MATMUL (TRANSPOSE(M), M) / Del_t CASE DEFAULT ! Method 1, Euler forward difference WORK = M / Del_t C = F + MATMUL (WORK, D) S = K + M / Del_t END SELECT IF ( debug_el_sq .AND. THIS_EL == 1 ) THEN PRINT *, 'K, M, F, WORK, D, S, C for method = ', METHOD CALL RPRINT (K, LT_FREE, LT_FREE, 1) CALL RPRINT (M, LT_FREE, LT_FREE, 1) CALL RPRINT (F, 1, LT_FREE, 1) CALL RPRINT (D, 1, LT_FREE, 1) CALL RPRINT (WORK, LT_FREE, LT_FREE, 1) CALL RPRINT (S, LT_FREE, LT_FREE, 1) CALL RPRINT (C, 1, LT_FREE, 1) END IF ! end file: 121.my_el_sq_inc ! suppress compiler warnings (touch is never true) IF ( TOUCH ) CALL TOUCH_UNUSED_ITEMS_1 (AJ, AJ_INV, B, BODY, DGH, & DGV, E, EB, G, H_INTG, IE, L_PT_PROP, PRT_L_PT, & PRT_MAT, STRAIN, STRAIN_0, STRESS, XYZ) END SUBROUTINE ELEM_SQ_EX_121 !FUNCTION START_DOF_VALUE_EX_121 (IG, XYZ) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! DEFINE STARTING VALUE OF PARAMETER IG IN TERMS OF !! COORDINATES OF THE NODE (FOR ITERATIVE SOLUTIONS) !! * * * * * * * * * * * * * * * * * * * * * * * * * * !! A PROBLEM DEPENDENT ROUTINE !Use System_Constants ! For EXAMPLE !IMPLICIT NONE ! INTEGER, INTENT(IN) :: IG ! local dof number ! REAL(DP), INTENT(IN) :: XYZ (N_SPACE) ! position ! REAL(DP) :: START_DOF_VALUE_EX_121 ! result !! IG = LOCAL PARAMETER NUMBER AT NODE !! ................................................................... !! ** START_DOF_VALUE PROBLEM DEPENDENT STATEMENTS FOLLOW ** !! For REAL (DP) :: START_DOF_VALUE ! result !! Given INTEGER, INTENT(IN) :: IG ! local dof number !! REAL(DP), INTENT(IN) :: XYZ (N_SPACE) ! position !! Given modules: Select_source, System_Constants1 !! ................................................................... ! INTEGER, SAVE :: NOTE = 0 ! IF ( NOTE == 0 ) THEN ; NOTE = 1; ! IF ( DEBUG_INCLUDE ) PRINT *, 'NOTE: USED START_DOF_VALUE_EX_121' ! PRINT *,'NOTE: DEFAULT START VALUE IS ', INITIAL_VALUE ! END IF !!--> ASSIGN UNIT STARTING VALUE ONLY ! START_DOF_VALUE_EX_121 = INITIAL_VALUE !! suppress compiler warnings (touch is never true) ! IF ( TOUCH ) PRINT *,'START_DOF_VALUE:', IG, XYZ !END FUNCTION START_DOF_VALUE_EX_121 !============= End Files for EXAMPLE number 121 ============= ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| ! ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||| !============= Begin Files for EXAMPLE number 122 ============= SUBROUTINE DESCRIBE_EXAMPLE_122 ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! LIST A DESCRIPTION OF THE EXAMPLE APPLICATION,ITS ! TYPICAL DATA AND ANY EXACT SOLUTIONS AVAILABLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants ! For EXAMPLE = 122 PRINT *, "*** DESCRIPTIONS OF EXAMPLE ", EXAMPLE, " ***" PRINT *,