! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE DERIV_3_L (X, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIND LOCAL DERIVATIVES FOR A 3 NODE LINE ELEMENT ! IN NATURAL COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: DH (3) ! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS (SHAPE_3_L) ! X = LOCAL COORDINATE OF POINT, -1 TO +1 ! LOCAL NODE COORD. ARE -1,0,+1 1----2----3 DH (1) = X - 0.5d0 DH (2) = - 2.d0 * X DH (3) = X + 0.5d0 END SUBROUTINE DERIV_3_L SUBROUTINE DERIV2_3_L (X, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIND LOCAL SECOND DERIVATIVES FOR A 3 NODE LINE ELEMENT ! IN NATURAL COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: D2H (3) ! D2H = LOCAL DERIVATIVES OF SHAPE FUNCTIONS (SHAPE_3_L) ! X = LOCAL COORDINATE OF POINT, -1 TO +1 ! LOCAL NODE COORD. ARE -1,0,+1. 1----2----3 D2H (1) = 1.d0 D2H (2) = - 2.d0 D2H (3) = 1.d0 END SUBROUTINE DERIV2_3_L SUBROUTINE DERIV_6_T (S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES FOR A SIX NODE UNIT TRIANGLE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: S, T REAL(DP), INTENT(OUT) :: DH (2, 6) ! S,T = LOCAL COORDINATES, SEE SHAPE_6_T ! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS ! DH(1,K) = DH(K)/DS, DH(2,K)=DH(K)/DT ! NODAL COORDS : 1-(0,0) 2-(1,0) 3-(0,1) ! 4-(0.5,0) 5-(0.5,0.5) 6-(0,0.5) DH (1, 1) = - 3.d0 + 4.d0 * S + 4.d0 * T DH (1, 2) = - 1.d0 + 4.d0 * S DH (1, 3) = 0.d0 DH (1, 4) = 4.d0 - 8.d0 * S - 4.d0 * T DH (1, 5) = 4.d0 * T DH (1, 6) = - 4.d0 * T DH (2, 1) = - 3.d0 + 4.d0 * S + 4.d0 * T DH (2, 2) = 0.d0 DH (2, 3) = - 1.d0 + 4.d0 * T DH (2, 4) = - 4.d0 * S DH (2, 5) = 4.d0 * S DH (2, 6) = 4.d0 - 4.d0 * S - 8.d0 * T END SUBROUTINE DERIV_6_T SUBROUTINE DERIV2_6_T (S, T, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL SECOND DERIVATIVES FOR A SIX NODE UNIT TRIANGLE ! H,SS H,ST H,TT ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: S, T REAL(DP), INTENT(OUT) :: D2H (3, 6) ! S,T = LOCAL COORDINATES, SEE SHAPE_6_T ! D2H = LOCAL DERIVATIVES OF SHAPE FUNCTIONS ! D2H(1,K) = D2H(K)/DS_DS, D2H(2,K)=D2H(K)/DS_DT, ETC ! NODAL COORDS : 1-(0,0) 2-(1,0) 3-(0,1) ! 4-(0.5,0) 5-(0.5,0.5) 6-(0,0.5) ! | V,AA | | V,II | A,B IMPLY LOCAL COORD ! | V,AB | = (SDJ) | V,IJ | I,J IMPLY GLOBAL COORD ! | V,BB | | V,JJ | D2H (1, 1) = 4.d0 D2H (1, 2) = 4.d0 D2H (1, 3) = 0.d0 D2H (1, 4) = -8.d0 D2H (1, 5) = 0.d0 D2H (1, 6) = 0.d0 D2H (2, 1) = 4.d0 D2H (2, 2) = 0.d0 D2H (2, 3) = 0.d0 D2H (2, 4) = -4.d0 D2H (2, 5) = 4.d0 D2H (2, 6) = -4.d0 D2H (3, 1) = 4.d0 D2H (3, 2) = 0.d0 D2H (3, 3) = 4.d0 D2H (3, 4) = 0.d0 D2H (3, 5) = 0.d0 D2H (3, 6) = -8.d0 END SUBROUTINE DERIV2_6_T SUBROUTINE DERIV_8_Q (S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIND LOCAL DERIVATIVES OF SHAPE FUNCTIONS FOR AN ! EIGHT NODE PARAMETRIC QUADRILATERAL ELEMENT ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: S, T REAL(DP), INTENT(OUT) :: DH (2, 8) REAL(DP) :: S_P, S_M, T_P, T_M ! S,T = LOCAL COORDINATES OF POINT, SEE SHAPE_8_Q ! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS, H ! DH(1,J) = DH(J)/DS, DH(2,J) = DH(J)/DT ! H = SHAPE FUNCTION ARRAY S_P = 1.d0 + S ; S_M = 1.d0 - S T_P = 1.d0 + T ; T_M = 1.d0 - T DH (1, 1) = - 0.25d0 * T_M * (S_M + S_M + T_M - 3.d0) DH (2, 1) = - 0.25d0 * S_M * (T_M + S_M + T_M - 3.d0) DH (1, 2) = 0.25d0 * T_M * (S_P + S_P + T_M - 3.d0) DH (2, 2) = - 0.25d0 * S_P * (T_M + S_P + T_M - 3.d0) DH (1, 3) = 0.25d0 * T_P * (S_P + S_P + T_P - 3.d0) DH (2, 3) = 0.25d0 * S_P * (T_P + S_P + T_P - 3.d0) DH (1, 4) = - 0.25d0 * T_P * (S_M + S_M + T_P - 3.d0) DH (2, 4) = 0.25d0 * S_M * (T_P + S_M + T_P - 3.d0) DH (1, 5) = - S * T_M DH (2, 5) = - 0.5d0 * (1.d0 - S * S) DH (1, 6) = 0.5d0 * (1.d0 - T * T) DH (2, 6) = - T * S_P DH (1, 7) = - S * T_P DH (2, 7) = 0.5d0 * (1.d0 - S * S) DH (1, 8) = - 0.5d0 * (1.d0 - T * T) DH (2, 8) = - T * S_M END SUBROUTINE DERIV_8_Q SUBROUTINE DCHECK (DELTA, NOD_PER_EL, N_SPACE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CHECKING OF THE LOCAL COORDINATE DERIVATIVES OF THE ! NOD_PER_EL SHAPE FUNCTIONS AT A LOCAL POINT FOR A C0 ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE REAL(DP), INTENT(IN) :: DELTA (N_SPACE, NOD_PER_EL) REAL(DP) :: D_SUM, TOL = 1.0D-7 INTEGER :: I_ERR, J, N_PRT = 6 ! DELTA = LOCAL DERIVATIVES OF SHAPE FUNCTIONS ! NOD_PER_EL = NUMBER OF SHAPE FUNCTIONS ! N_SPACE = DIMENSION OF LOCAL SPACE I_ERR = 0 DO J = 1, N_SPACE D_SUM = SUM (DELTA (J, 1:NOD_PER_EL)) IF (ABS (D_SUM) > TOL) THEN IF (I_ERR == 0) WRITE (N_PRT, *) 'DERIVATIVES INCORRECT' I_ERR = 1 ; WRITE (N_PRT, * ) 'ROW, D_SUM', J, D_SUM END IF END DO IF (I_ERR /= 0) THEN CALL RPRINT (DELTA, NOD_PER_EL, N_SPACE, 1) WRITE (N_PRT, * ) 'END OF WARNING FROM DCHECK' END IF END SUBROUTINE DCHECK SUBROUTINE DERIV (PT, DLH, NOD_PER_EL, N_SPACE, L_SHAPE, & N_G_DOF) !b N_G_DOF, ELEM_NODES) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C0 INTERPOLATION LOCAL DERIVATIVES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, L_SHAPE, N_G_DOF !b INTEGER, INTENT(IN) :: ELEM_NODES (NOD_PER_EL) REAL(DP), INTENT(IN) :: PT (N_SPACE) REAL(DP), INTENT(OUT) :: DLH (N_SPACE, NOD_PER_EL * N_G_DOF) ! DLH = LOCAL DERIVATIVES AT PT !b ! ELEM_NODES = TOPOLOGY LIST, IF VARIABLE ! L_SHAPE = 1-LINE, 2-TRI, 3-QUAD, 4-HEX, 5-TET, 6-WEDGE ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NUMBER OF SPATIAL DIMENSIONS ! PT = LOCAL COORD OF A POINT ! BRANCH ON SHAPE, THEN NUMBER OF NODES IF ( L_SHAPE == 1 ) THEN !--> 1-D ELEMENTS IF ( NOD_PER_EL == 2 ) CALL DERIV_2_L (PT (1), DLH) IF ( NOD_PER_EL == 3 ) CALL DERIV_3_L (PT(1),DLH) RETURN ELSEIF ( L_SHAPE == 2 ) THEN !--> TRIANGULAR 2-D ELEMENTS IF ( NOD_PER_EL == 3 ) CALL DERIV_3_T (PT (1), PT (2), DLH) ! IF ( NOD_PER_EL == 4 ) CALL DERIV_4_T (PT(1),PT(2),DLH) IF ( NOD_PER_EL == 6 ) CALL DERIV_6_T (PT (1), PT (2), DLH) ! IF ( NOD_PER_EL == 7 ) CALL DERIV_7_T (PT(1),PT(2),DLH) IF ( NOD_PER_EL == 10 ) CALL DERIV_10_T (PT(1),PT(2),DLH) IF ( NOD_PER_EL == 15 ) CALL DERIV_15_T (PT(1),PT(2),DLH) RETURN ELSEIF ( L_SHAPE == 3 ) THEN !--> QUADRILATERAL 2-D ELEMENTS IF ( NOD_PER_EL == 4 ) CALL DERIV_4_Q (PT (1), PT (2), DLH) IF ( NOD_PER_EL == 8 ) CALL DERIV_8_Q (PT (1), PT (2), DLH) IF ( NOD_PER_EL == 9 ) CALL DERIV_9_Q (PT (1), PT (2), DLH) IF ( NOD_PER_EL == 12 ) CALL DERIV_12_Q (PT (1), PT (2), DLH) IF ( NOD_PER_EL == 16 ) CALL DERIV_16_Q (PT (1), PT (2), DLH) IF ( NOD_PER_EL == 17 ) CALL DERIV_17_Q (PT (1), PT (2), DLH) ! IF ( NOD_PER_EL == 25 ) CALL DERIV_25_Q (PT (1), PT (2), DLH) RETURN ELSEIF ( L_SHAPE == 4) THEN !--> HEXAHEDRA 3-D ELEMENTS IF ( NOD_PER_EL == 8 ) CALL DERIV_8_H (PT (1), PT (2), PT (3), DLH) IF ( NOD_PER_EL == 20 ) CALL DERIV_20_H (PT (1), PT (2), PT (3), DLH) !b ! IF ( NOD_PER_EL == 20 ) CALL DERIV_8_20_H (PT(1),PT(2),PT(3),& !b ! DLH,ELEM_NODES) ! IF ( NOD_PER_EL == 27 ) CALL DERIV_27_H (PT(1),PT(2),PT(3),DLH) ! IF ( NOD_PER_EL == 32 ) CALL DERIV_32_H (PT(1),PT(2),PT(3),DLH) RETURN ELSEIF ( L_SHAPE == 5) THEN !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) IF ( NOD_PER_EL == 4 ) CALL DERIV_4_P (PT(1),PT(2),PT(3),DLH) ! IF ( NOD_PER_EL == 10 ) CALL DERIV_10_P (PT(1),PT(2),PT(3),DLH) ! IF ( NOD_PER_EL == 20 ) CALL DERIV_20_P (PT(1),PT(2),PT(3),DLH) RETURN ELSEIF ( L_SHAPE == 6) THEN !--> WEDGE 3-D ELEMENTS STOP 'NO WEDGE IN SHAPE' ! IF ( NOD_PER_EL == 6 ) CALL DERIV_6_W (PT(1),PT(2),PT(3),DLH) ! IF ( NOD_PER_EL == 15 ) CALL DERIV_15_W (PT(1),PT(2),PT(3),DLH) ! RETURN ELSEIF ( L_SHAPE == 7 ) THEN !--> USER SUPPLIED ELEMENT !b ! CALL DERIV_USER (PT(1),PT(2),PT(3),DLH,ELEM_NODES) STOP 'NO USER ELEMENT IN DERIV' ELSEIF ( L_SHAPE < 1 .OR. L_SHAPE > 7 ) THEN !--> UNSUPPORTED OPTION STOP 'UNSUPPORTED ELEMENT IN DERIV' END IF END SUBROUTINE DERIV SUBROUTINE DERIV_6_T_BAR (R, S, T, DBH) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! BARACENTRIC DERIVATIVES FOR A T6 ELEMENT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: DBH (3, 6) REAL(DP), PARAMETER :: ZERO = 0.D0, ONE = 1.D0, FOUR = 4.D0 ! DBH = BARACENTRIC DERIVATIVES OF H ! R,S,T = BARACENTRIC COORDINATES, R+S+T=1 ! NODE R S T NODE R S T ORDER: ! 1 1 0 0 4 1/2 1/2 0 3 ! 2 0 1 0 5 0 1/2 1/2 6 5 ! 3 0 0 1 6 1/2 0 1/2 1 4 2 ! BARACENTRIC DERIVATIVES DBH (1, 1) = FOUR * R - ONE DBH (2, 1) = ZERO DBH (3, 1) = ZERO DBH (1, 2) = ZERO DBH (2, 2) = FOUR * S - ONE DBH (3, 2) = ZERO DBH (1, 3) = ZERO DBH (2, 3) = ZERO DBH (3, 3) = FOUR * T - ONE DBH (1, 4) = ZERO DBH (2, 4) = FOUR * T DBH (3, 4) = FOUR * S DBH (1, 5) = FOUR * T DBH (2, 5) = ZERO DBH (3, 5) = FOUR * R DBH (1, 6) = FOUR * S DBH (2, 6) = FOUR * R DBH (3, 6) = ZERO END SUBROUTINE DERIV_6_T_BAR SUBROUTINE DERIV_USER (PT, DH, NOD_PER_EL, N_SPACE, N_G_DOF) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! USER DEFINED ELEMENT C0 INTERPOLATION FUNCTION DERIVATIVES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF REAL(DP), INTENT(IN) :: PT (N_SPACE) REAL(DP), INTENT(OUT) :: DH (N_SPACE, NOD_PER_EL * N_G_DOF) ! DH = LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS AT PT ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! PT = LOCAL COORD OF A POINT ! DH = 0.d0 PRINT *, 'DERIV_USER: USER ERROR, DERIVATIVES NOT SUPPLIED FOR' PRINT *, 'ELEMENT WITH ', NOD_PER_EL, ' NODES,' PRINT *, 'IN ', N_SPACE, ' DIMENSIONS, AND WITH' PRINT *, N_G_DOF, ' DEGREES OF FREEDOM PER NODE.' PRINT *, 'WARNING: CHECK VALUE OF shape KEYWORD' STOP 'No user defined element interpolations, DERIV_USER' END SUBROUTINE DERIV_USER SUBROUTINE DERIV2_USER (PT, D2H, NOD_PER_EL, N_SPACE, N_G_DOF, N_2_DER) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! USER DEFINED ELEMENT C0 INTERPOLATION FUNCTION SECOND ! DERIVATIVES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF, N_2_DER REAL(DP), INTENT(IN) :: PT (N_SPACE) REAL(DP), INTENT(OUT) :: D2H (N_2_DER, NOD_PER_EL * N_G_DOF) INTEGER, SAVE :: WARN = 0 ! D2H = LOCAL SECOND DERIVATIVES OF FUNCTIONS AT PT ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! N_2_DER = NUMBER OF SECOND DERIVATIVES ! PT = LOCAL COORD OF A POINT ! D2H = 0.d0 IF ( WARN == 0 ) THEN WARN = 1 PRINT *, 'DERIV2_USER: USER ERROR, SECOND DERIVATIVES NOT SUPPLIED' PRINT *, 'FOR ELEMENT TYPE WITH ', NOD_PER_EL, ' NODES, IN' PRINT *, N_SPACE, ' DIMENSIONS, AND WITH ', N_G_DOF, ' DEGREES OF' PRINT *, ' FREEDOM PER NODE. CHECK VALUE OF shape KEYWORD' PRINT *, 'SECOND DERIVATIVES SET TO ZERO' END IF !b STOP 'No user defined element interpolations, DERIV2_USER' END SUBROUTINE DERIV2_USER SUBROUTINE DERIV_C_N_USER (PT, LENGTH, DH, NOD_PER_EL, & N_SPACE, N_G_DOF) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! USER DEFINED ELEMENT C_N INTERPOLATION FUNCTION DERIVATIVES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF REAL(DP), INTENT(IN) :: PT (N_SPACE), LENGTH REAL(DP), INTENT(OUT) :: DH (N_SPACE, NOD_PER_EL * N_G_DOF) ! DH = LOCAL DERIVATIVES OF C1 INTERPOLATIONS AT PT ! LENGTH = PHYSICAL LENGTH ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! PT = LOCAL COORD OF A POINT ! DH = 0.d0 STOP 'No user defined element interpolations, DERIV_C_N_USER' END SUBROUTINE DERIV_C_N_USER SUBROUTINE SHAPE_USER (PT, H, NOD_PER_EL, N_SPACE, N_G_DOF) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! USER DEFINED ELEMENT C0 INTERPOLATION FUNCTIONS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF REAL(DP), INTENT(IN) :: PT (N_SPACE) REAL(DP), INTENT(OUT) :: H (NOD_PER_EL * N_G_DOF) ! H = ELEMENT C0 INTERPOLATION FUNCTIONS AT PT ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! PT = LOCAL COORD OF A POINT H = 0.d0 PRINT *, 'USER ERROR, INTERPOLATIONS NOT SUPPLIED FOR' PRINT *, 'ELEMENT WITH ', NOD_PER_EL, ' NODES,' PRINT *, 'IN ', N_SPACE, ' DIMENSIONS, AND WITH' PRINT *, N_G_DOF, ' DEGREES OF FREEDOM PER NODE.' PRINT *, 'WARNING: CHECK VALUE OF shape KEYWORD' STOP 'No user defined element interpolations, SHAPE_USER' END SUBROUTINE SHAPE_USER SUBROUTINE SHAPE_C_N_USER (PT, LENGTH, H, NOD_PER_EL, & N_SPACE, N_G_DOF) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! USER DEFINED ELEMENT C_N INTERPOLATION FUNCTIONS ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF REAL(DP), INTENT(IN) :: PT (N_SPACE), LENGTH REAL(DP), INTENT(OUT) :: H (NOD_PER_EL * N_G_DOF) ! H = LOCAL C1 INTERPOLATION FUNCTIONS AT PT ! LENGTH = PHYSICAL LENGTH ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! PT = LOCAL COORD OF A POINT H = 0.d0 STOP 'No user defined element interpolations, SHAPE_C_N_USER' END SUBROUTINE SHAPE_C_N_USER SUBROUTINE SCALAR_DERIVS (PT_L, DH_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C0 ELEMENT DERIVATIVE FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM) REAL(DP), INTENT(OUT) :: DH_L (LT_PARM, LT_N) ! DH_L = ELEMENT INTERPOLATION DERIVATIVES AT PT_L ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT ELEMENT DERIVATIVES' N_WARN = N_WARN + 1 ; CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( LT_N ) CASE (2) SELECT CASE ( N_G_DOF ) CASE (1) ; CALL DERIV_2_L (PT_L (1), DH_L(1,:)) CASE DEFAULT CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF) END SELECT CASE (3) ; CALL DERIV_3_L (PT_L (1), DH_L(1,:)) CASE (4) ; CALL DERIV_4_L (PT_L (1), DH_L(1,:)) CASE DEFAULT CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( LT_N ) !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; CALL DERIV_3_T (PT_L (1), PT_L (2), DH_L) ! CASE ( 4) ; CALL DERIV_4_T (PT_L (1), PT_L (2), DH_L) CASE ( 6) ; CALL DERIV_6_T (PT_L (1), PT_L (2), DH_L) ! CASE ( 7) ; CALL DERIV_7_T (PT_L (1), PT_L (2), DH_L) CASE (10) ; CALL DERIV_10_T (PT_L (1), PT_L (2), DH_L) CASE (15) ; CALL DERIV_15_T (PT_L (1), PT_L (2), DH_L) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL DERIV_4_Q (PT_L (1), PT_L (2), DH_L) CASE ( 8) ; CALL DERIV_8_Q (PT_L (1), PT_L (2), DH_L) CASE ( 9) ; CALL DERIV_9_Q (PT_L (1), PT_L (2), DH_L) CASE (12) ; CALL DERIV_12_Q (PT_L (1), PT_L (2), DH_L) CASE (16) ; CALL DERIV_16_Q (PT_L (1), PT_L (2), DH_L) CASE (17) ; CALL DERIV_17_Q (PT_L (1), PT_L (2), DH_L) ! CASE (25) ; CALL DERIV_25_Q (PT_L (1), PT_L (2), DH_L) CASE DEFAULT CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( LT_N ) !--> HEXAHEDRA 3-D ELEMENTS CASE ( 8) ; CALL DERIV_8_H (PT_L (1), PT_L (2), PT_L (3), DH_L) CASE (20) ; CALL DERIV_20_H (PT_L (1), PT_L (2), PT_L (3), DH_L) ! CASE (27) ; CALL DERIV_27_H (PT_L (1), PT_L (2), PT_L (3), DH_L) ! CASE (32) ; CALL DERIV_32_H (PT_L (1), PT_L (2), PT_L (3), DH_L) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) CASE ( 4) ; CALL DERIV_4_P (PT_L (1), PT_L (2), PT_L (3), DH_L) ! CASE (10) ; CALL DERIV_10_P (PT_L (1), PT_L (2), PT_L (3), DH_L) ! CASE (20) ; CALL DERIV_20_P (PT_L (1), PT_L (2), PT_L (3), DH_L) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL DERIV_6_W (PT_L (1), PT_L (2), PT_L (3), DH_L) ! CASE (15) ; CALL DERIV_15_W (PT_L (1), PT_L (2), PT_L (3), DH_L) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL DERIV_USER (PT_L, DH_L, LT_N, LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_N, LT_PARM, N_G_DOF STOP 'INVALID SPACE SIZE IN SCALAR_DERIVS' END SELECT ! number of spatial dimensions END SUBROUTINE SCALAR_DERIVS SUBROUTINE SCALAR_2ND_DERIVS (PT_L, D2H_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C0 ELEMENT 2ND DERIVATIVE FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM) REAL(DP), INTENT(OUT) :: D2H_L (N_2_DER, LT_N) ! D2H_L = ELEMENT INTERPOLATION SECOND DERIVATIVES AT PT_L ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT ELEMENT 2ND DERIVATIVES' N_WARN = N_WARN + 1 ; CASE (1) ! 1-D space !--> 1-D ELEMENTS ! | V,AA | = (P_AJ) | V,II | A,B IMPLY LOCAL COORD ! I,J IMPLY GLOBAL COORD SELECT CASE ( LT_N ) CASE (2) SELECT CASE ( N_G_DOF ) CASE (1) ; D2H_L = 0.d0 CASE DEFAULT CALL DERIV2_USER (PT_L, D2H_L, LT_N, LT_PARM, & N_G_DOF, N_2_DER) END SELECT CASE (3) ; D2H_L (1, :) = (/ 1.d0, -2.d0, 1.d0 /) CASE (4) ; CALL DERIV2_4_L (PT_L (1), D2H_L(1,:)) CASE DEFAULT CALL DERIV2_USER (PT_L, D2H_L, LT_N, LT_PARM, & N_G_DOF, N_2_DER) END SELECT ! number of element nodes CASE (2) ! 2-D space ! | V,AA | | V,II | A,B IMPLY LOCAL COORD ! | V,AB | = (P_AJ) | V,IJ | I,J IMPLY GLOBAL COORD ! | V,BB | | V,JJ | SELECT CASE ( LT_N ) ! NUMBER OF NODES !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; D2H_L = 0.d0 ! CASE ( 4) ; CALL DERIV2_4_T (PT_L (1), PT_L (2), D2H_L) CASE ( 6) ; CALL DERIV2_6_T (PT_L (1), PT_L (2), D2H_L) ! CASE ( 7) ; CALL DERIV2_7_T (PT_L (1), PT_L (2), D2H_L) ! CASE (10) ; CALL DERIV2_10_T (PT_L (1), PT_L (2), D2H_L) !b CASE (15) ; CALL DERIV2_15_T (PT_L (1), PT_L (2), D2H_L) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL DERIV2_4_Q (PT_L (1), PT_L (2), D2H_L) !b CASE ( 8) ; CALL DERIV2_8_Q (PT_L (1), PT_L (2), D2H_L) CASE ( 9) ; CALL DERIV2_9_Q (PT_L (1), PT_L (2), D2H_L) ! CASE (12) ; CALL DERIV2_12_Q (PT_L (1), PT_L (2), D2H_L) ! CASE (16) ; CALL DERIV2_16_Q (PT_L (1), PT_L (2), D2H_L) ! CASE (17) ; CALL DERIV2_17_Q (PT_L (1), PT_L (2), D2H_L) ! CASE (25) ; CALL DERIV2_25_Q (PT_L (1), PT_L (2), D2H_L) CASE DEFAULT CALL DERIV2_USER (PT_L,D2H_L,LT_N,LT_PARM,N_G_DOF,N_2_DER) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( LT_N ) !--> HEXAHEDRA 3-D ELEMENTS !b CASE ( 8) ; CALL DERIV2_8_H (PT_L (1), PT_L (2), PT_L (3), & !b D2H_L) !b CASE (20) ; CALL DERIV2_20_H (PT_L (1), PT_L (2), PT_L (3), & !b D2H_L) ! CASE (27) ; CALL DERIV2_27_H (PT_L (1), PT_L (2), PT_L (3), & ! D2H_L) ! CASE (32) ; CALL DERIV2_32_H (PT_L (1), PT_L (2), PT_L (3), & ! D2H_L) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) ! CASE ( 4) ; D2H_L = 0.d0 ! D2H_L) ! CASE (10) ; CALL DERIV2_10_P (PT_L (1), PT_L (2), PT_L (3), & ! D2H_L) ! CASE (20) ; CALL DERIV2_20_P (PT_L (1), PT_L (2), PT_L (3), & ! D2H_L) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL DERIV2_6_W (PT_L (1), PT_L (2), PT_L (3), & ! D2H_L) ! CASE (15) CALL DERIV2_15_W (PT_L (1), PT_L (2), PT_L (3), & ! D2H_L) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL DERIV2_USER (PT_L, D2H_L, LT_N, LT_PARM, N_G_DOF, N_2_DER) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_N, LT_PARM, N_G_DOF STOP 'INVALID SPACE SIZE IN SCALAR_2ND_DERIVS' END SELECT ! number of spatial dimensions END SUBROUTINE SCALAR_2ND_DERIVS SUBROUTINE ELEM_C_N_DERIVS (PT_L, LENGTH, DH_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C_N ELEMENT DERIVATIVE FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM), LENGTH REAL(DP), INTENT(OUT) :: DH_L (LT_PARM, LT_N) ! DH_L = ELEMENT INTERPOLATION DERIVATIVES AT PT_L ! LENGTH = PHYSICAL LENGTH OF ELEMENT ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT ELEMENT C_N DERIVATIVES' N_WARN = N_WARN + 1 ; CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( LT_N ) CASE (2) ! 2 node line element SELECT CASE ( N_G_DOF ) CASE (2) ; CALL DERIV_C1_L (PT_L (1), LENGTH, DH_L(1,:)) CASE (3) ; CALL DERIV_C2_L (PT_L (1), LENGTH, DH_L(1,:)) CASE (4) ; CALL DERIV_C3_L (PT_L (1), LENGTH, DH_L(1,:)) CASE DEFAULT CALL DERIV_C_N_USER (PT_L, LENGTH, DH_L, LT_N, & LT_PARM, N_G_DOF) END SELECT CASE DEFAULT CALL DERIV_C_N_USER (PT_L, LENGTH, DH_L, LT_N, & LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_N, LT_PARM, N_G_DOF STOP 'INVALID SPACE SIZE IN ELEM_C_N_DERIVS' END SELECT ! number of spatial dimensions END SUBROUTINE ELEM_C_N_DERIVS SUBROUTINE GEN_ELEM_DERIV (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C0 ELEMENT DERIVATIVE FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF REAL(DP), INTENT(IN) :: PT (N_SPACE) REAL(DP), INTENT(OUT) :: DH_L (N_SPACE, NOD_PER_EL * N_G_DOF) ! DH_L = ELEMENT INTERPOLATION DERIVATIVES AT PT ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! PT = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES (and then N_G_DOF) SELECT CASE ( N_SPACE ) CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( NOD_PER_EL ) CASE (2) ; CALL DERIV_2_L (PT (1), DH_L(1,:)) CASE (3) ; CALL DERIV_3_L (PT (1), DH_L(1,:)) CASE (4) ; CALL DERIV_4_L (PT (1), DH_L(1,:)) CASE DEFAULT CALL DERIV_USER (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( NOD_PER_EL ) !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; CALL DERIV_3_T (PT (1), PT (2), DH_L) ! CASE ( 4) ; CALL DERIV_4_T (PT (1), PT (2), DH_L) CASE ( 6) ; CALL DERIV_6_T (PT (1), PT (2), DH_L) ! CASE ( 7) ; CALL DERIV_7_T (PT (1), PT (2), DH_L) CASE (10) ; CALL DERIV_10_T (PT (1), PT (2), DH_L) CASE (15) ; CALL DERIV_15_T (PT (1), PT (2), DH_L) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL DERIV_4_Q (PT (1), PT (2), DH_L) CASE ( 8) ; CALL DERIV_8_Q (PT (1), PT (2), DH_L) CASE ( 9) ; CALL DERIV_9_Q (PT (1), PT (2), DH_L) CASE (12) ; CALL DERIV_12_Q (PT (1), PT (2), DH_L) CASE (16) ; CALL DERIV_16_Q (PT (1), PT (2), DH_L) CASE (17) ; CALL DERIV_17_Q (PT (1), PT (2), DH_L) ! CASE (25) ; CALL DERIV_25_Q (PT (1), PT (2), DH_L) CASE DEFAULT CALL DERIV_USER (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( NOD_PER_EL ) !--> HEXAHEDRA 3-D ELEMENTS CASE ( 8) ; CALL DERIV_8_H (PT (1), PT (2), PT (3), DH_L) CASE (20) ; CALL DERIV_20_H (PT (1), PT (2), PT (3), DH_L) ! CASE (27) ; CALL DERIV_27_H (PT (1), PT (2), PT (3), DH_L) ! CASE (32) ; CALL DERIV_32_H (PT (1), PT (2), PT (3), DH_L) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) CASE ( 4) ; CALL DERIV_4_P (PT (1), PT (2), PT (3), DH_L) ! CASE (10) ; CALL DERIV_10_P (PT (1), PT (2), PT (3), DH_L) ! CASE (20) ; CALL DERIV_20_P (PT (1), PT (2), PT (3), DH_L) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL DERIV_6_W (PT (1), PT (2), PT (3), DH_L) ! CASE (15) ; CALL DERIV_15_W (PT (1), PT (2), PT (3), DH_L) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL DERIV_USER (PT, DH_L, NOD_PER_EL, N_SPACE, N_G_DOF) END SELECT ! number of element nodes CASE DEFAULT STOP 'INVALID SPACE SIZE IN GEN_ELEM_DERIV' END SELECT ! number of spatial dimensions END SUBROUTINE GEN_ELEM_DERIV SUBROUTINE SCALAR_SHAPES (PT_L, H_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C0 ELEMENT INTERPOLATION FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM) REAL(DP), INTENT(OUT) :: H_L (LT_N) ! H_L = ELEMENT INTERPOLATION FUNCTIONS AT PT_L ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT ELEMENT SHAPES' N_WARN = N_WARN + 1 ; H = 1.d0 CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( LT_N ) CASE (2) SELECT CASE ( N_G_DOF ) CASE (1) ; CALL SHAPE_2_L (PT_L (1), H_L) CASE DEFAULT CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF) END SELECT CASE (3) ; CALL SHAPE_3_L (PT_L (1), H_L) CASE (4) ; CALL SHAPE_4_L (PT_L (1), H_L) CASE DEFAULT CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( LT_N ) !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; CALL SHAPE_3_T (PT_L (1), PT_L (2), H_L) ! CASE ( 4) ; CALL SHAPE_4_T (PT_L (1), PT_L (2), H_L) CASE ( 6) ; CALL SHAPE_6_T (PT_L (1), PT_L (2), H_L) ! CASE ( 7) ; CALL SHAPE_7_T (PT_L (1), PT_L (2), H_L) CASE (10) ; CALL SHAPE_10_T (PT_L (1), PT_L (2), H_L) CASE (15) ; CALL SHAPE_15_T (PT_L (1), PT_L (2), H_L) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL SHAPE_4_Q (PT_L (1), PT_L (2), H_L) CASE ( 8) ; CALL SHAPE_8_Q (PT_L (1), PT_L (2), H_L) CASE ( 9) ; CALL SHAPE_9_Q (PT_L (1), PT_L (2), H_L) CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), H_L) CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), H_L) CASE (17) ; CALL SHAPE_17_Q (PT_L (1), PT_L (2), H_L) ! CASE (25) ; CALL SHAPE_25_Q (PT_L (1), PT_L (2), H_L) CASE DEFAULT PRINT *, 'WARNING: CHECK VALUE OF shape KEYWORD' N_WARN = N_WARN + 1 CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( LT_N ) !--> HEXAHEDRA 3-D ELEMENTS CASE ( 8) ; CALL SHAPE_8_H (PT_L (1), PT_L (2), PT_L (3), H_L) CASE (20) ; CALL SHAPE_20_H (PT_L (1), PT_L (2), PT_L (3), H_L) ! CASE (27) ; CALL SHAPE_27_H (PT_L (1), PT_L (2), PT_L (3), H_L) ! CASE (32) ; CALL SHAPE_32_H (PT_L (1), PT_L (2), PT_L (3), H_L) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) CASE ( 4) ; CALL SHAPE_4_P (PT_L (1), PT_L (2), PT_L (3), H_L) ! CASE (10) ; CALL SHAPE_10_P (PT_L (1), PT_L (2), PT_L (3), H_L) ! CASE (20) ; CALL SHAPE_20_P (PT_L (1), PT_L (2), PT_L (3), H_L) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL SHAPE_6_W (PT_L (1), PT_L (2), PT_L (3), H_L) ! CASE (15) ; CALL SHAPE_15_W (PT_L (1), PT_L (2), PT_L (3), H_L) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL SHAPE_USER (PT_L, H_L, LT_N, LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_N, LT_PARM, N_G_DOF STOP 'INVALID SPACE SIZE IN SCALAR_SHAPES' END SELECT ! number of spatial dimensions END SUBROUTINE SCALAR_SHAPES SUBROUTINE ELEM_C_N_SHAPES (PT_L, LEN_1, H_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C_N ELEMENT INTERPOLATION FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM), LEN_1 REAL(DP), INTENT(OUT) :: H_L (LT_N) !b XXX must be bigger ! H_L = ELEMENT INTERPOLATION FUNCTIONS AT PT_L ! LEN_1 = PHYSICAL LENGTH IN PARAMETRIC DIRECTION 1 ! LT_N = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT ELEMENT C_N SHAPES' N_WARN = N_WARN + 1 ; H = 1.d0 CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( LT_N ) CASE (2) SELECT CASE ( N_G_DOF ) CASE (2) ; CALL SHAPE_C1_L (PT_L (1), LEN_1, H_L) CASE (3) ; CALL SHAPE_C2_L (PT_L (1), LEN_1, H_L) CASE (4) ; CALL SHAPE_C3_L (PT_L (1), LEN_1, H_L) CASE DEFAULT CALL SHAPE_C_N_USER (PT_L, LEN_1, H_L, LT_N, & LT_PARM, N_G_DOF) END SELECT CASE DEFAULT CALL SHAPE_C_N_USER (PT_L, LEN_1, H_L, LT_N, & LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( LT_N ) !--> RECTANGULAR 2-D ELEMENTS ! CASE (17) ! CALL SHAPE_16_R (PT_L (1), PT_L (2), LEN_1, LEN_2, H_L) CASE DEFAULT CALL SHAPE_C_N_USER (PT_L, LEN_1, H_L, LT_N, & LT_PARM, N_G_DOF) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_N, LT_PARM, N_G_DOF STOP 'INVALID SPACE SIZE IN ELEM_C_N_SHAPES' END SELECT ! number of spatial dimensions END SUBROUTINE ELEM_C_N_SHAPES SUBROUTINE GEOM_DERIVS (PT_L, DLG_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! C0 ELEMENT GEOMETRIC DERIVATIVE FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM) REAL(DP), INTENT(OUT) :: DLG_L (LT_PARM, LT_GEOM) ! DLG_L = LOCAL DERIVATIVES OF GEOMETRIC INTERPOLATIONS ! LT_GEOM = NUMBER OF GEOMETRIC NODES PER ELEMENT ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT GEOMETRIC DERIVATIVES' N_WARN = N_WARN + 1 ; CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( LT_GEOM ) CASE (2) ; CALL DERIV_2_L (PT_L (1), DLG_L) CASE (3) ; CALL DERIV_3_L (PT_L (1), DLG_L) CASE (4) ; CALL DERIV_4_L (PT_L (1), DLG_L) CASE DEFAULT CALL DERIV_USER (PT_L, DLG_L, LT_GEOM, LT_PARM, 1) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( LT_GEOM ) !--> LINE ELEMENT IN 2_D CASE ( 2) ; CALL DERIV_2_L (PT_L (1), DLG_L) !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; CALL DERIV_3_T (PT_L (1), PT_L (2), DLG_L) ! CASE ( 4) ; CALL DERIV_4_T (PT_L (1), PT_L (2), DLG_L) CASE ( 6) ; CALL DERIV_6_T (PT_L (1), PT_L (2), DLG_L) ! CASE ( 7) ; CALL DERIV_7_T (PT_L (1), PT_L (2), DLG_L) CASE (10) ; CALL DERIV_10_T (PT_L (1), PT_L (2), DLG_L) CASE (15) ; CALL DERIV_15_T (PT_L (1), PT_L (2), DLG_L) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL DERIV_4_Q (PT_L (1), PT_L (2), DLG_L) CASE ( 8) ; CALL DERIV_8_Q (PT_L (1), PT_L (2), DLG_L) CASE ( 9) ; CALL DERIV_9_Q (PT_L (1), PT_L (2), DLG_L) CASE (12) ; CALL DERIV_12_Q (PT_L (1), PT_L (2), DLG_L) CASE (16) ; CALL DERIV_16_Q (PT_L (1), PT_L (2), DLG_L) CASE (17) ; CALL DERIV_17_Q (PT_L (1), PT_L (2), DLG_L) ! CASE (25) ; CALL DERIV_25_Q (PT_L (1), PT_L (2), DLG_L) CASE DEFAULT CALL DERIV_USER (PT_L, DLG_L, LT_GEOM, LT_PARM, 1) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( LT_GEOM ) !--> LINE ELEMENT IN 3_D CASE ( 2) ; CALL DERIV_2_L (PT_L (1), DLG_L) !--> TRIANGULAR 2-D ELEMENTS IN 3-D CASE ( 3) ; CALL DERIV_3_T (PT_L (1), PT_L (2), DLG_L) CASE ( 6) ; CALL DERIV_6_T (PT_L (1), PT_L (2), DLG_L) !--> QUADRILATERAL 2-D ELEMENTS IN 3-D, OR TET CASE ( 4) IF ( LT_SHAP == 3 ) CALL DERIV_4_Q (PT_L (1), PT_L (2), DLG_L) IF ( LT_SHAP == 5 ) CALL DERIV_4_P (PT_L (1), PT_L (2), & PT_L (3), DLG_L) CASE ( 9) ; CALL DERIV_9_Q (PT_L (1), PT_L (2), DLG_L) CASE (12) ; CALL DERIV_12_Q (PT_L (1), PT_L (2), DLG_L) CASE (16) ; CALL DERIV_16_Q (PT_L (1), PT_L (2), DLG_L) !--> HEXAHEDRA 3-D ELEMENTS CASE ( 8) ; CALL DERIV_8_H (PT_L (1), PT_L (2), PT_L (3), DLG_L) CASE (20) ; CALL DERIV_20_H (PT_L (1), PT_L (2), PT_L (3), DLG_L) ! CASE (27) ; CALL DERIV_27_H (PT_L (1), PT_L (2), PT_L (3), DLG_L) ! CASE (32) ; CALL DERIV_32_H (PT_L (1), PT_L (2), PT_L (3), DLG_L) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) ! CASE (10) ; CALL DERIV_10_P (PT_L (1), PT_L (2), PT_L (3), DLG_L) ! CASE (20) ; CALL DERIV_20_P (PT_L (1), PT_L (2), PT_L (3), DLG_L) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL DERIV_6_W (PT_L (1), PT_L (2), PT_L (3), DLG_L) ! CASE (15) ; CALL DERIV_15_W (PT_L (1), PT_L (2), PT_L (3), DLG_L) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL DERIV_USER (PT_L, DLG_L, LT_GEOM, LT_PARM, 1) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_GEOM, LT_PARM STOP 'INVALID SPACE SIZE IN GEOM_DERIVS' END SELECT ! number of spatial dimensions END SUBROUTINE GEOM_DERIVS SUBROUTINE GEOM_SHAPES (PT_L, G_L) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! C0 ELEMENT GEOMETRIC INTERPOLATION FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data IMPLICIT NONE REAL(DP), INTENT(IN) :: PT_L (LT_PARM) REAL(DP), INTENT(OUT) :: G_L (LT_GEOM) ! G_L = ELEMENT INTERPOLATION FUNCTIONS AT PT_L ! LT_GEOM = NUMBER OF GEOMETRIC NODES PER ELEMENT ! LT_PARM = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES SELECT CASE ( LT_PARM ) CASE (0) ! point PRINT *,'WARNING, NO POINT GEOMETRIC SHAPES' N_WARN = N_WARN + 1 ; G_L = 1.d0 CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( LT_GEOM ) CASE (2) ; CALL SHAPE_2_L (PT_L (1), G_L) CASE (3) ; CALL SHAPE_3_L (PT_L (1), G_L) CASE (4) ; CALL SHAPE_4_L (PT_L (1), G_L) CASE DEFAULT CALL SHAPE_USER (PT_L, G_L, LT_GEOM, LT_PARM, 1) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( LT_GEOM ) !--> LINE ELEMENT IN 2-D CASE (2) ; CALL SHAPE_2_L (PT_L (1), G_L) !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; CALL SHAPE_3_T (PT_L (1), PT_L (2), G_L) ! CASE ( 4) ; CALL SHAPE_4_T (PT_L (1), PT_L (2), G_L) CASE ( 6) ; CALL SHAPE_6_T (PT_L (1), PT_L (2), G_L) ! CASE ( 7) ; CALL SHAPE_7_T (PT_L (1), PT_L (2), G_L) CASE (10) ; CALL SHAPE_10_T (PT_L (1), PT_L (2), G_L) CASE (15) ; CALL SHAPE_15_T (PT_L (1), PT_L (2), G_L) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL SHAPE_4_Q (PT_L (1), PT_L (2), G_L) CASE ( 8) ; CALL SHAPE_8_Q (PT_L (1), PT_L (2), G_L) CASE ( 9) ; CALL SHAPE_9_Q (PT_L (1), PT_L (2), G_L) CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), G_L) CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), G_L) CASE (17) ; CALL SHAPE_17_Q (PT_L (1), PT_L (2), G_L) ! CASE (25) ; CALL SHAPE_25_Q (PT_L (1), PT_L (2), G_L) CASE DEFAULT CALL SHAPE_USER (PT_L, G_L, LT_GEOM, LT_PARM, 1) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( LT_GEOM ) !--> LINE ELEMENT IN 3-D CASE (2) ; CALL SHAPE_2_L (PT_L (1), G_L) !--> TRIANGULAR 2-D ELEMENTS IN 3-D CASE ( 3) ; CALL SHAPE_3_T (PT_L (1), PT_L (2), G_L) CASE ( 6) ; CALL SHAPE_6_T (PT_L (1), PT_L (2), G_L) !--> QUADRILATERAL 2-D ELEMENTS IN 3-D, OR TET CASE ( 4) IF ( LT_SHAP == 3 ) CALL SHAPE_4_Q (PT_L (1), PT_L (2), G_L) IF ( LT_SHAP == 5 ) CALL SHAPE_4_P (PT_L (1), PT_L (2), & PT_L (3), G_L) CASE ( 9) ; CALL SHAPE_9_Q (PT_L (1), PT_L (2), G_L) CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), G_L) CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), G_L) !--> HEXAHEDRA 3-D ELEMENTS CASE ( 8) ; CALL SHAPE_8_H (PT_L (1), PT_L (2), PT_L (3), G_L) CASE (20) ; CALL SHAPE_20_H (PT_L (1), PT_L (2), PT_L (3), G_L) ! CASE (27) ; CALL SHAPE_27_H (PT_L (1), PT_L (2), PT_L (3), G_L) ! CASE (32) ; CALL SHAPE_32_H (PT_L (1), PT_L (2), PT_L (3), G_L) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) !b CASE ( 4) ; CALL SHAPE_4_P (PT_L (1), PT_L (2), PT_L (3), G_L) ! CASE (10) ; CALL SHAPE_10_P (PT_L (1), PT_L (2), PT_L (3), G_L) ! CASE (20) ; CALL SHAPE_20_P (PT_L (1), PT_L (2), PT_L (3), G_L) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL SHAPE_6_W (PT_L (1), PT_L (2), PT_L (3), G_L) ! CASE (15) ; CALL SHAPE_15_W (PT_L (1), PT_L (2), PT_L (3), G_L) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL SHAPE_USER (PT_L, G_L, LT_GEOM, LT_PARM, 1) END SELECT ! number of element nodes CASE DEFAULT PRINT *, LT_GEOM, LT_PARM STOP 'INVALID SPACE SIZE IN GEOM_SHAPES' END SELECT ! number of spatial dimensions END SUBROUTINE GEOM_SHAPES SUBROUTINE GEN_ELEM_SHAPE (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EVALUATE C0 ELEMENT INTERPOLATION FUNCTIONS IN 1D, 2D, 3D ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: NOD_PER_EL, N_SPACE, N_G_DOF REAL(DP), INTENT(IN) :: PT_L (N_SPACE) REAL(DP), INTENT(OUT) :: H (NOD_PER_EL * N_G_DOF) ! H = ELEMENT INTERPOLATION FUNCTIONS AT PT_L ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_G_DOF = NUMBER OF DEGREES OF FREEDOM PER NODE ! N_SPACE = NO OF SPATIAL DIMENSIONS ! PT_L = LOCAL COORD OF A POINT ! BRANCH ON SPACE, THEN NUMBER OF NODES (and then N_G_DOF) SELECT CASE ( N_SPACE ) CASE (1) ! 1-D space !--> 1-D ELEMENTS SELECT CASE ( NOD_PER_EL ) CASE (2) ; CALL SHAPE_2_L (PT_L (1), H) CASE (3) ; CALL SHAPE_3_L (PT_L (1), H) CASE (4) ; CALL SHAPE_4_L (PT_L (1), H) CASE DEFAULT CALL SHAPE_USER (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF) END SELECT ! number of element nodes CASE (2) ! 2-D space SELECT CASE ( NOD_PER_EL ) !--> LINE ELEMENTS IN 2-D CASE ( 2) ; CALL SHAPE_2_L (PT_L (1), H) !--> TRIANGULAR 2-D ELEMENTS CASE ( 3) ; CALL SHAPE_3_T (PT_L (1), PT_L (2), H) ! CASE ( 4) ; CALL SHAPE_4_T (PT_L (1), PT_L (2), H) CASE ( 6) ; CALL SHAPE_6_T (PT_L (1), PT_L (2), H) ! CASE ( 7) ; CALL SHAPE_7_T (PT_L (1), PT_L (2), H) CASE (10) ; CALL SHAPE_10_T (PT_L (1), PT_L (2), H) CASE (15) ; CALL SHAPE_15_T (PT_L (1), PT_L (2), H) !--> QUADRILATERAL 2-D ELEMENTS CASE ( 4) ; CALL SHAPE_4_Q (PT_L (1), PT_L (2), H) CASE ( 8) ; CALL SHAPE_8_Q (PT_L (1), PT_L (2), H) CASE ( 9) ; CALL SHAPE_9_Q (PT_L (1), PT_L (2), H) CASE (12) ; CALL SHAPE_12_Q (PT_L (1), PT_L (2), H) CASE (16) ; CALL SHAPE_16_Q (PT_L (1), PT_L (2), H) CASE (17) ; CALL SHAPE_17_Q (PT_L (1), PT_L (2), H) ! CASE (25) ; CALL SHAPE_25_Q (PT_L (1), PT_L (2), H) CASE DEFAULT CALL SHAPE_USER (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF) END SELECT ! number of element nodes CASE (3) ! 3-D space SELECT CASE ( NOD_PER_EL ) !--> LINE ELEMENTS IN 3-D CASE ( 2) ; CALL SHAPE_2_L (PT_L (1), H) !--> HEXAHEDRA 3-D ELEMENTS CASE ( 8) ; CALL SHAPE_8_H (PT_L (1), PT_L (2), PT_L (3), H) CASE (20) ; CALL SHAPE_20_H (PT_L (1), PT_L (2), PT_L (3), H) ! CASE (27) ; CALL SHAPE_27_H (PT_L (1), PT_L (2), PT_L (3), H) ! CASE (32) ; CALL SHAPE_32_H (PT_L (1), PT_L (2), PT_L (3), H) !--> TETRAHEDRA 3-D ELEMENTS (PYRAMIDS) CASE ( 4) ; CALL SHAPE_4_P (PT_L (1), PT_L (2), PT_L (3), H) ! CASE (10) ; CALL SHAPE_10_P (PT_L (1), PT_L (2), PT_L (3), H) ! CASE (20) ; CALL SHAPE_20_P (PT_L (1), PT_L (2), PT_L (3), H) !--> WEDGE 3-D ELEMENTS ! CASE ( 6) ; CALL SHAPE_6_W (PT_L (1), PT_L (2), PT_L (3), H) ! CASE (15) ; CALL SHAPE_15_W (PT_L (1), PT_L (2), PT_L (3), H) !--> USER SUPPLIED ELEMENT CASE DEFAULT CALL SHAPE_USER (PT_L, H, NOD_PER_EL, N_SPACE, N_G_DOF) END SELECT ! number of element nodes CASE DEFAULT STOP 'INVALID SPACE SIZE IN GEN_ELEM_SHAPE' END SELECT ! number of spatial dimensions END SUBROUTINE GEN_ELEM_SHAPE SUBROUTINE DERIV_10_T (R, S, DH) ! ********************************************************** ! LOCAL DERIVATIVES OF SHAPE FUNCTIONS FOR 10-NODED TRIANGLE ! ********************************************************** Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DH (2, 10) REAL(DP) :: T0, T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, Z ! S ! THE CUBIC TRIANGLE : ! 3 ! R,S = LOCAL COORDS OF RT | \ ! H = ELEM SHARE FUNCTIONS | \ ! 6 8 ! ELEM_NODES = ELEM INCIDENCES LIST | \ ! ELEMENT SKETCH TO RIGHT 9 10 5 ! | \ ! 1@(0,0) 2@(1,0) 2@(0,1) 1---4---7---2 ... R Z = 1.d0 - R - S T0 = 4.5d0 * (6.d0 * Z - 1.d0) T1 = 4.5d0 * S * (3.d0 * S - 1.d0) T2 = 4.5d0 * Z * (3.d0 * Z - 1.d0) T3 = 4.5d0 * R * (3.d0 * R - 1.d0) T4 = 4.5d0 * (6.d0 * S - 1.d0) T5 = 4.5d0 * (6.d0 * R - 1.d0) T6 = T0 * S T7 = T0 * R T8 = -0.5d0 * (27.d0 * Z * Z - 18.d0 * Z + 2.d0) T9 = 0.5d0 * (27.d0 * S * S - 18.d0 * S + 2.d0) T10 = 0.5d0 * (27.d0 * R * R - 18.d0 * R + 2.d0) DH (1, 1 ) = T8 DH (1, 2 ) = T10 DH (1, 3 ) = 0.d0 DH (1, 4 ) = T2 - T7 DH (1, 5 ) = S * T5 DH (1, 6 ) = -T1 DH (1, 7 ) = Z * T5 - T3 DH (1, 8 ) = T1 DH (1, 9 ) = -T6 DH (1, 10) = 27.d0 * S * (Z - R) DH (2, 1 ) = T8 DH (2, 2 ) = 0.d0 DH (2, 3 ) = T9 DH (2, 4 ) = -T7 DH (2, 5 ) = T3 DH (2, 6 ) = Z * T4 - T1 DH (2, 7 ) = -T3 DH (2, 8 ) = R * T4 DH (2, 9 ) = T2 - T6 DH (2, 10) = 27.d0 * R * (Z - S) END SUBROUTINE DERIV_10_T SUBROUTINE SHAPE_10_T (R, S, H) ! ****************************************************** ! SHAPE FUNCTIONS FOR 10-NODED TRIANGLE ! ****************************************************** Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (10) REAL(DP) :: T1, T2, T3, Z ! S ! THE CUBIC TRIANGLE : ! 3 ! R,S = LOCAL COORDS OF PT | \ ! H = ELEM SHAPE FUNCTIONS | \ ! 6 8 ! ELEM_NODES = ELEM INCIDENCES LIST | \ ! ELEMENT SKETCH TO RIGHT 9 10 5 ! | \ ! 1@(0,0) 2@(1,0) 2@(0,1) 1---4---7---2 ... R Z = 1.d0 - R - S T1 = S * (3.d0 * S - 1.d0) T2 = Z * (3.d0 * Z - 1.d0) T3 = R * (3.d0 * R - 1.d0) H (1 ) = 0.5d0 * T2 * (3.d0 * Z - 2.d0) H (2 ) = 0.5d0 * T3 * (3.d0 * R - 2.d0) H (3 ) = 0.5d0 * T1 * (3.d0 * S - 2.d0) H (4 ) = 4.5d0 * R * T2 H (5 ) = 4.5d0 * S * T3 H (6 ) = 4.5d0 * Z * T1 H (7 ) = 4.5d0 * Z * T3 H (8 ) = 4.5d0 * R * T1 H (9 ) = 4.5d0 * S * T2 H (10) = 27.d0 * S * Z * R END SUBROUTINE SHAPE_10_T SUBROUTINE DERIV_12_Q (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! DERIVATIVES OF A 12 NODE QUAD, SEE SHAPE_12_Q ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DH (2, 12) REAL(DP), PARAMETER :: F1 = 1.d0/32, F2 = 9.d0/32 REAL(DP) :: R_P, R_M, S_P, S_M, R2, S2 ! R,S = LOCAL COORDS OF PT 4--11---7---3 ! H = ELEM SHAPE FUNCTIONS I I ! 8 S 10 ! ELEMENT SKETCH TO RIGHT I .R I ! LOCAL COORD OF NODES: 12 6 ! 1-(-1,-1) 3-(+1,+1) I I ! SIDES OF ORDER 1, 2, OR 3 1---5---9---2 R_P = 1 + R ; R_M = 1 - R ; S_P = 1 + S ; S_M = 1 - S R2 = R*R ; S2 = S*S DH (1, 1) = F1*S_M*( 10 - 27*R2 + 18*R - 9*S2) DH (2, 1) = F1*R_M*( 10 - 27*S2 + 18*S - 9*R2) DH (1, 2) = F1*S_M*(-10 + 27*R2 + 18*R + 9*S2) DH (2, 2) = F1*R_P*( 10 - 27*S2 + 18*S - 9*R2) DH (1, 3) = F1*S_P*(-10 + 27*R2 + 18*R + 9*S2) DH (2, 3) = F1*R_P*(-10 + 27*S2 + 18*S + 9*R2) DH (1, 4) = F1*S_P*( 10 - 27*R2 + 18*R - 9*S2) DH (2, 4) = F1*R_M*(-10 + 27*S2 + 18*S + 9*R2) DH (1, 5) = F2*S_M*( 9*R2 - 2*R - 3) DH (2, 5) = -F2*(1 - R2)*(1 - 3*R) DH (1, 6) = F2*(1 - S2)*(1 - 3*S) DH (2, 6) = F2*R_P*( 9*S2 - 2*S - 3) DH (1, 7) = F2*S_P*(-9*R2 - 2*R + 3) DH (2, 7) = F2*(1 - R2)*(1 + 3*R) DH (1, 8) = -F2*(1 - S2)*(1 + 3*S) DH (2, 8) = F2*R_M*(-9*S2 - 2*S + 3) DH (1, 9) = F2*S_M*(-9*R2 - 2*R + 3) DH (2, 9) = -F2*(1 - R2)*(1 + 3*R) DH (1, 10) = F2*(1 - S2)*(1 + 3*S) DH (2, 10) = F2*R_P*(-9*S2 - 2*S + 3) DH (1, 11) = F2*S_P*( 9*R2 - 2*R - 3) DH (2, 11) = F2*(1 - R2)*(1 - 3*R) DH (1, 12) = -F2*(1 - S2)*(1 - 3*S) DH (2, 12) = F2*R_M*( 9*S2 - 2*S - 3) END SUBROUTINE DERIV_12_Q SUBROUTINE DERIV_15_T (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES FOR 15-NODED TRIANGLE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), PARAMETER :: PT1 = 42.666666666666667d0, & PT2 = 10.666666666666667d0 REAL(DP), INTENT(OUT) :: DH (2, 15) REAL(DP) :: CC, T1, T2, T3, T4, T5, T6, T7, T8, T9 CC = 1.d0 - R - S T1 = R - 0.25d0 ; T2 = R - 0.5d0 ; T3 = R - 0.75d0 T4 = S - 0.25d0 ; T5 = S - 0.5d0 ; T6 = S - 0.75d0 T7 = CC - 0.25d0 ; T8 = CC - 0.5d0 ; T9 = CC - 0.75d0 DH (1, 2) = PT2 * (T2 * T3 * (T1 + R) + R * T1 * (T3 + T2) ) DH (1, 5) = PT1 * S * (T2 * (T1 + R) + R * T1) DH (1, 8) = 64.d0 * S * T4 * (T1 + R) DH (1, 11) = PT1 * S * T4 * T5 DH (1, 3) = 0.d0 DH (1, 6) = - PT1 * S * T4 * T5 DH (1, 9) = - 64.d0 * S * T4 * (T7 + CC) DH (1, 12) = - PT1 * S * (T8 * (T7 + CC) + CC * T7) DH (1, 1) = - PT2 * (T8 * T9 * (T7 + CC) + CC * T7 * (T8 + T9) ) DH (1, 4) = PT1 * (CC * T7 * T8 - R * (T8 * (T7 + CC) + CC * T7) ) DH (1, 7) = 64.d0 * (CC * T7 * (T1 + R) - R * T1 * (T7 + CC) ) DH (1, 10) = PT1 * (CC * (T2 * (T1 + R) + R * T1) - R * T1 * T2) DH (1, 14) = 128.d0 * S * (CC * (T1 + R) - R * T1) DH (1, 15) = 128.d0 * S * T4 * (CC - R) DH (1, 13) = 128.d0 * S * (CC * T7 - R * (T7 + CC) ) DH (2, 2) = 0.d0 DH (2, 5) = PT1 * R * T1 * T2 DH (2, 8) = 64.d0 * R * T1 * (T4 + S) DH (2, 11) = PT1 * R * (T5 * (T4 + S) + S * T4) DH (2, 3) = PT2 * (T5 * T6 * (T4 + S) + S * T4 * (T6 + T5) ) DH (2, 6) = PT1 * ((CC * (T5 * (T4 + S) + S * T4)) - S * T4 * T5) DH (2, 9) = 64.d0 * (CC * T7 * (T4 + S) - S * T4 * (T7 + CC) ) DH (2, 12) = PT1 * (CC * T7 * T8 - S * (T8 * (T7 + CC) + CC * T7)) DH (2, 1) = - PT2 * (T8 * T9 * (T7 + CC) + CC * T7 * (T8 + T9) ) DH (2, 4) = - PT1 * R * (T8 * (T7 + CC) + CC * T7) DH (2, 7) = - 64.d0 * R * T1 * (T7 + CC) DH (2, 10) = - PT1 * R * T1 * T2 DH (2, 14) = 128.d0 * R * T1 * (CC - S) DH (2, 15) = 128.d0 * R * (CC * (T4 + S) - S * T4) DH (2, 13) = 128.d0 * R * (CC * T7 - S * (CC + T7) ) END SUBROUTINE DERIV_15_T SUBROUTINE DERIV_16_Q (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! DERIVATIVES OF A LAGRANGIAN 16 NODE QUAD ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! SEE SHAPE_16_Q Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DH (2, 16) REAL(DP) :: HR (4), HS (4), DHR (4), DHS (4) ! R,S = LOCAL COORDS OF PT 4--11---7---3 ! H = ELEM SHAPE FUNCTIONS I I ! 8 16 S 15 10 ! ELEMENT SKETCH TO RIGHT I .R I ! LOCAL COORD OF NODES: 12 13 14 6 ! 1-(-1,-1) 3-(+1,+1) I I ! 1---5---9---2 ! Evaluate the 1D interpolations CALL SHAPE_4_L (R, HR) CALL SHAPE_4_L (S, HS) CALL DERIV_4_L (R, DHR) CALL DERIV_4_L (S, DHS) DH (1, 1) = DHR (1)*HS (1) DH (2, 1) = HR (1)*DHS (1) DH (1, 2) = DHR (4)*HS (1) DH (2, 2) = HR (4)*DHS (1) DH (1, 3) = DHR (4)*HS (4) DH (2, 3) = HR (4)*DHS (4) DH (1, 4) = DHR (1)*HS (4) DH (2, 4) = HR (1)*DHS (4) DH (1, 5) = DHR (2)*HS (1) DH (2, 5) = HR (2)*DHS (1) DH (1, 6) = DHR (4)*HS (2) DH (2, 6) = HR (4)*DHS (2) DH (1, 7) = DHR (3)*HS (4) DH (2, 7) = HR (3)*DHS (4) DH (1, 8) = DHR (1)*HS (3) DH (2, 8) = HR (1)*DHS (3) DH (1, 9) = DHR (3)*HS (1) DH (2, 9) = HR (3)*DHS (1) DH (1, 10) = DHR (4)*HS (3) DH (2, 10) = HR (4)*DHS (3) DH (1, 11) = DHR (2)*HS (4) DH (2, 11) = HR (2)*DHS (4) DH (1, 12) = DHR (1)*HS (2) DH (2, 12) = HR (1)*DHS (2) DH (1, 13) = DHR (2)*HS (2) DH (2, 13) = HR (2)*DHS (2) DH (1, 14) = DHR (3)*HS (2) DH (2, 14) = HR (3)*DHS (2) DH (1, 15) = DHR (3)*HS (3) DH (2, 15) = HR (3)*DHS (3) DH (1, 16) = DHR (2)*HS (3) DH (2, 16) = HR (2)*DHS (3) END SUBROUTINE DERIV_16_Q SUBROUTINE DERIV_16_QS (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! LOCAL DERIVATIVES FOR A SERENDIPITY 16 NODE QUAD ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! SEE SHAPE_16_QS FOR NODE LOCATIONS Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), PARAMETER :: PT667 = 0.66666666666667d0, & PT1D12 = 0.083333333333333d0, HALF = 0.5d0 REAL(DP), INTENT(OUT) :: DH (2, 16) REAL(DP) :: RR, SS, RRR, SSS, RP, RM, S_P, SM, R1 RR = R * R ; SS = S * S RRR = R * R * R ; SSS = S * S * S RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S R1 = R - 1.d0 DH (1, 1) = PT1D12*SM*(16.d0*RRR - 12.d0*RR - 2.d0*R + 4.d0*SSS - S + 4.d0) DH (1, 2) = PT1D12*SM*(16.d0*RRR + 12.d0*RR - 2.d0*R - 4.d0*SSS + S - 4.d0) DH (1, 3) = PT1D12*S_P*(16.d0*RRR + 12.d0*RR - 2.d0*R + 4.d0*SSS - S - 4.d0) DH (1, 4) = PT1D12*S_P*(16.d0*RRR - 12.d0*RR - 2.d0*R - 4.d0*SSS + S + 4.d0) DH (1, 5) = PT667 * (4.d0 * R - 1.d0 + 3.d0 * RR - 8.d0 * RRR) * SM DH (1, 6) = PT667 * S * SM * S_P * ( -1.d0 + 2.d0 * S) DH (1, 7) = PT667 * (1.d0 + 4.d0 * R - 3.d0 * RR - 8.d0 * RRR) * S_P DH (1, 8) = -PT667 * S * SM * S_P * (1.d0 + 2.d0 * S) DH (1, 9) = R * ( -5.d0 + 8.d0 * RR) * SM DH (1, 10) = HALF * SM * S_P * (1.d0 - 2.d0 * S) * (1.d0 + 2.d0 * S) DH (1, 11) = R * ( -5.d0 + 8.d0 * RR) * S_P DH (1, 12) = HALF * SM * S_P * ( -1.d0 + 2.d0 * S) * (1.d0 + 2.d0 * S) DH (1, 13) = PT667 * (1.d0 + 4.d0 * R - 3.d0 * RR - 8.d0 * RRR) * SM DH (1, 14) = PT667 * S * SM * S_P * (1.d0 + 2.d0 * S) DH (1, 15) = PT667 * (4.d0 * R - 1.d0 + 3.d0 * RR - 8.d0 * RRR) * S_P DH (1, 16) = PT667 * S * SM * S_P * (1.d0 - 2.d0 * S) DH (2, 1) = PT1D12 * RM * (16.d0 * SSS - 12.d0 * SS - 2.d0 * S & + 4.d0 * RRR - R + 4.d0) DH (2, 2) = PT1D12 * RP * (16.d0 * SSS - 12.d0 * SS - 2.d0 * S & - 4.d0 * RRR + R + 4.d0) DH (2, 3) = PT1D12 * RP * (16.d0 * SSS + 12.d0 * SS - 2.d0 * S & + 4.d0 * RRR - R - 4.d0) DH (2, 4) = PT1D12 * R1 * ( -16.d0 * SSS - 12.d0 * SS & + 2.d0 * S + 4.d0 * RRR - R + 4.d0) DH (2, 5) = PT667 * R * R1 * RP * (2.d0 * R - 1.d0) DH (2, 6) = PT667 * ( -1.d0 + 4.d0 * S + 3.d0 * SS - 8.d0 * SSS) * RP DH (2, 7) = PT667 * R * RM * RP * (1.d0 + 2.d0 * R) DH (2, 8) = PT667 * ( -1.d0 - 4.d0 * S + 3.d0 * SS + 8.d0 * SSS) * R1 DH (2, 9) = HALF * RM * RP * (2.d0 * R - 1.d0) * (1.d0 + 2.d0 * R) DH (2, 10) = S * ( -5.d0 + 8.d0 * SS) * RP DH (2, 11) = HALF * R1 * RP * (2.d0 * R - 1.d0) * (1.d0 + 2.d0 * R) DH (2, 12) = S * (5.d0 - 8.d0 * SS) * R1 DH (2, 13) = PT667 * R * R1 * RP * (1.d0 + 2.d0 * R) DH (2, 14) = PT667 * (1.d0 + 4.d0 * S - 3.d0 * SS - 8.d0 * SSS) * RP DH (2, 15) = PT667 * R * RM * RP * (2.d0 * R - 1.d0) DH (2, 16) = PT667 * (1.d0 - 4.d0 * S - 3.d0 * SS + 8.d0 * SSS) * R1 END SUBROUTINE DERIV_16_QS SUBROUTINE DERIV_16_R (R, S, A, B, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! WARNING, this code is not fully checked. No known bugs. ! FIRST DERIVATIVES OF A C1 RECTANGULAR ELEMENT IN ! UNIT COORDINATES USING TENSOR PRODUCTS OF 1D BASIS ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, A, B REAL(DP), INTENT(OUT) :: DH (2, 16) REAL(DP) :: HR (4), DHR (4), HS (4), DHS (4) ! DOF ARE W W,X W,Y W,XY AT EACH NODE (N_G_DOF=4) ! X // R, Y // S. S ! A = PHYSICAL LENGTH IN X 4 -------- 3 ! B = PHYSICAL LENGTH IN Y I I ! R,S = LOCAL UNIT COORDS I I ! 1@(0,0), 3@(1,1) 1 -------- 2 ->R ! Evaluate the 1D interpolations CALL SHAPE_C1_L (R, A, HR) CALL SHAPE_C1_L (S, B, HS) CALL DERIV_C1_L (R, A, DHR) CALL DERIV_C1_L (S, B, DHS) ! Form tensor products for R direction DH (1, 1) = DHR (1) * HS (1) DH (1, 2) = DHR (2) * HS (1) DH (1, 3) = DHR (1) * HS (2) DH (1, 4) = DHR (2) * HS (2) DH (1, 5) = DHR (3) * HS (1) DH (1, 6) = DHR (4) * HS (1) DH (1, 7) = DHR (3) * HS (2) DH (1, 8) = DHR (4) * HS (2) DH (1, 9) = DHR (3) * HS (3) DH (1, 10) = DHR (4) * HS (3) DH (1, 11) = DHR (3) * HS (4) DH (1, 12) = DHR (4) * HS (4) DH (1, 13) = DHR (1) * HS (3) DH (1, 14) = DHR (2) * HS (3) DH (1, 15) = DHR (1) * HS (4) DH (1, 16) = DHR (2) * HS (4) ! Form tensor products for S direction DH (2, 1) = HR (1) * DHS (1) DH (2, 2) = HR (2) * DHS (1) DH (2, 3) = HR (1) * DHS (2) DH (2, 4) = HR (2) * DHS (2) DH (2, 5) = HR (3) * DHS (1) DH (2, 6) = HR (4) * DHS (1) DH (2, 7) = HR (3) * DHS (2) DH (2, 8) = HR (4) * DHS (2) DH (2, 9) = HR (3) * DHS (3) DH (2, 10) = HR (4) * DHS (3) DH (2, 11) = HR (3) * DHS (4) DH (2, 12) = HR (4) * DHS (4) DH (2, 13) = HR (1) * DHS (3) DH (2, 14) = HR (2) * DHS (3) DH (2, 15) = HR (1) * DHS (4) DH (2, 16) = HR (2) * DHS (4) END SUBROUTINE DERIV_16_R SUBROUTINE DERIV_17_Q (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES OF SERENDIPITY QUAD WITH 17 NODES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SEE SHAPE_17_Q FOR NODE LOCATIONS Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: PT667 = 0.66666666666667d0, & PT1D12 = 0.083333333333333d0, HALF = 0.5d0 REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DH (2, 17) REAL(DP) :: RR, SS, RRR, SSS, RP, RM, S_P, SM, R1 RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S RR = R * R ; SS = S * S RRR = R * R * R ; SSS = S * S * S R1 = R - 1. DH (1, 1) = PT1D12 * SM * (16.d0 * RRR - 12.d0 * RR & - 6.d0 * R * S - 8.d0 * R + 4.d0 * SSS - S + 4.d0) DH (1, 2) = PT1D12 * SM * (16.d0 * RRR + 12.d0 * RR & - 6.d0 * R * S - 8.d0 * R - 4.d0 * SSS + S - 4.d0) DH (1, 3) = PT1D12 * S_P * (16.d0 * RRR + 12.d0 * RR & + 6.d0 * R * S - 8.d0 * R + 4.d0 * SSS - S - 4.d0) DH (1, 4) = PT1D12 * S_P * (16.d0 * RRR - 12.d0 * RR & + 6.d0 * R * S - 8.d0 * R - 4.d0 * SSS + S + 4.d0) DH (1, 5) = -PT667 * (1.d0 - 4.d0 * R - 3.d0 * RR + 8.d0 * RRR) * SM DH (1, 6) = PT667 * S * SM * S_P * ( - 1.d0 + 2.d0 * S) DH (1, 7) = -PT667 * (-1.d0 - 4.d0 * R + 3.d0 * RR + 8.d0 * RRR) * S_P DH (1, 8) = - PT667 * S * SM * S_P * (1.d0 + 2.d0 * S) DH (1, 9) = R * SM * (8.d0 * RR + S - 4.d0) DH (1, 10) = HALF * SM * S_P * (2.d0 * R - 4.d0 * SS + 1.d0) DH (1, 11) = R * S_P * (8.d0 * RR - S - 4.d0) DH (1, 12) = HALF * SM * S_P * (2.d0 * R - 1.d0 + 4.d0 * SS) DH (1, 13) = PT667 * (1.d0 + 4.d0 * R - 3.d0 * RR - 8.d0 * RRR) * SM DH (1, 14) = PT667 * S * SM * S_P * (1.d0 + 2.d0 * S) DH (1, 15) = -PT667 * (1.d0 - 4.d0 * R - 3.d0 * RR + 8.d0 * RRR) * S_P DH (1, 16) = PT667 * S * SM * S_P * (1.d0 - 2.d0 * S) DH (1, 17) = - 2.d0 * R * SM * S_P DH (2, 1) = PT1D12 * RM * (16.d0 * SSS - 12.d0 * SS & - 6.d0 * R * S - 8.d0 * S + 4.d0 * RRR - R + 4.d0) DH (2, 2) = - PT1D12 * RP * ( -16.d0 * SSS + 12.d0 * SS & - 6.d0 * R * S + 8.d0 * S + 4.d0 * RRR - R - 4.d0) DH (2, 3) = PT1D12 * RP * (16.d0 * SSS + 12.d0 * SS & + 6.d0 * R * S - 8.d0 * S + 4.d0 * RRR - R - 4.d0) DH (2, 4) = PT1D12 * R1 * ( -16.d0 * SSS - 12.d0 * SS & + 6.d0 * R * S + 8.d0 * S + 4.d0 * RRR - R + 4.d0) DH (2, 5) = PT667 * R * R1 * RP * (2.d0 * R - 1.d0) DH (2, 6) = - PT667 * (1.d0 - 4.d0 * S - 3.d0 * SS + 8.d0 * SSS) * RP DH (2, 7) = - PT667 * R * R1 * RP * (1.d0 + 2.d0 * R) DH (2, 8) = PT667 * ( - 1.d0 - 4.d0 * S + 3.d0 * SS + 8.d0 * SSS) * R1 DH (2, 9) = - HALF * R1 * RP * (2.d0 * S - 1.d0 + 4.d0 * RR) DH (2, 10) = - S * RP * ( - 8.d0 * SS + R + 4.d0) DH (2, 11) = HALF * R1 * RP * ( - 2.d0 * S + 4.d0 * RR - 1.d0) DH (2, 12) = - S * R1 * (8.d0 * SS + R - 4.d0) DH (2, 13) = PT667 * R * R1 * RP * (1.d0 + 2.d0 * R) DH (2, 14) = -PT667 * (-1.d0 - 4.d0 * S + 3.d0 * SS + 8.d0 * SSS) * RP DH (2, 15) = -PT667 * R * R1 * RP * (2.d0 * R - 1.d0) DH (2, 16) = PT667 * (1.d0 - 4.d0 * S - 3.d0 * SS + 8.d0 * SSS) * R1 DH (2, 17) = 2.d0 * S * R1 * RP END SUBROUTINE DERIV_17_Q SUBROUTINE DERIV_8_20_H (R, S, T, DH, ELEM_NODES) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS FOR AN ! 8 TO 20 NODE HEXAHEDRON, SEE SHAPE_8_20_H FOR TOPOLOGY FIGURE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: DH (3, 20) INTEGER :: ELEM_NODES (20), I1 (20), I2 (20), K, K1, K2 REAL(DP) :: FR, FS, FT, RP, RM, S_P, SM, TP, TM, & RZ, SZ, TZ, DH1, DH2, DH3, FP, FM DATA I1 / 8 * 0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4 / DATA I2 / 8 * 0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8 / DATA FP, FM / 0.5d0, -0.5d0 / ! FOR PARAMETER DEFINITIONS SEE SUBROUTINE SHAPE_8_20_H ! R,S,T = LOCAL COORDINATES OF THE POINT -1 <= (R,S,T) <= +1 ! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS, 0 IF ELEM_NODES(I) = 0 ! ELEM_NODES = ARRAY OF ELEMENT INCIDENCES, ! IF ELEM_NODES(I)=0 THEN LOCAL NODE I IS NOT CONSIDERED IN ANALYSIS ! I1, I2 = CORNER NODES OF TWELVE EDGES RP = 0.5d0 * (1.d0 + R) ; S_P = 0.5d0 * (1.d0 + S) ; TP = 0.5d0 * (1.d0 + T) RM = 0.5d0 * (1.d0 - R) ; SM = 0.5d0 * (1.d0 - S) ; TM = 0.5d0 * (1.d0 - T) RZ = 1.d0 - R * R ; SZ = 1.d0 - S * S ; TZ = 1.d0 - T * T FR = - 2.d0 * R ; FS = - 2.d0 * S ; FT = - 2.d0 * T ! DERIVATIVES OF TRI-LINEAR CORNERS DH (1, 1) = TP * S_P * FP DH (2, 1) = TP * FP * RP DH (3, 1) = FP * S_P * RP DH (1, 2) = TP * S_P * FM DH (2, 2) = TP * FP * RM DH (3, 2) = FP * S_P * RM DH (1, 3) = TP * SM * FM DH (2, 3) = TP * FM * RM DH (3, 3) = FP * SM * RM DH (1, 4) = TP * SM * FP DH (2, 4) = TP * FM * RP DH (3, 4) = FP * SM * RP DH (1, 5) = TM * S_P * FP DH (2, 5) = TM * FP * RP DH (3, 5) = FM * S_P * RP DH (1, 6) = TM * S_P * FM DH (2, 6) = TM * FP * RM DH (3, 6) = FM * S_P * RM DH (1, 7) = TM * SM * FM DH (2, 7) = TM * FM * RM DH (3, 7) = FM * SM * RM DH (1, 8) = TM * SM * FP DH (2, 8) = TM * FM * RP DH (3, 8) = FM * SM * RP ! DERIVATIVES OF EDGE BUBBLES DH (1, 9) = TP * S_P * FR * 0.5d0 DH (2, 9) = TP * FP * RZ * 0.5d0 DH (3, 9) = FP * S_P * RZ * 0.5d0 DH (1, 10) = TP * SZ * FM * 0.5d0 DH (2, 10) = TP * FS * RM * 0.5d0 DH (3, 10) = FP * SZ * RM * 0.5d0 DH (1, 11) = TP * SM * FR * 0.5d0 DH (2, 11) = TP * FM * RZ * 0.5d0 DH (3, 11) = FP * SM * RZ * 0.5d0 DH (1, 12) = TP * SZ * FP * 0.5d0 DH (2, 12) = TP * FS * RP * 0.5d0 DH (3, 12) = FP * SZ * RP * 0.5d0 DH (1, 13) = TM * S_P * FR * 0.5d0 DH (2, 13) = TM * FP * RZ * 0.5d0 DH (3, 13) = FM * S_P * RZ * 0.5d0 DH (1, 14) = TM * SZ * FM * 0.5d0 DH (2, 14) = TM * FS * RM * 0.5d0 DH (3, 14) = FM * SZ * RM * 0.5d0 DH (1, 15) = TM * SM * FR * 0.5d0 DH (2, 15) = TM * FM * RZ * 0.5d0 DH (3, 15) = FM * SM * RZ * 0.5d0 DH (1, 16) = TM * SZ * FP * 0.5d0 DH (2, 16) = TM * FS * RP * 0.5d0 DH (3, 16) = FM * SZ * RP * 0.5d0 DH (1, 17) = TZ * S_P * FP * 0.5d0 DH (2, 17) = TZ * FP * RP * 0.5d0 DH (3, 17) = FT * S_P * RP * 0.5d0 DH (1, 18) = TZ * S_P * FM * 0.5d0 DH (2, 18) = TZ * FP * RM * 0.5d0 DH (3, 18) = FT * S_P * RM * 0.5d0 DH (1, 19) = TZ * SM * FM * 0.5d0 DH (2, 19) = TZ * FM * RM * 0.5d0 DH (3, 19) = FT * SM * RM * 0.5d0 DH (1, 20) = TZ * SM * FP * 0.5d0 DH (2, 20) = TZ * FM * RP * 0.5d0 DH (3, 20) = FT * SM * RP * 0.5d0 ! LOOP OVER TWELVE EDGES DO K = 9, 20 IF (ELEM_NODES (K) == 0) THEN ! ZERO EDGE BUBBLE DERIVATIVES DH (1, K) = 0.d0 DH (2, K) = 0.d0 DH (3, K) = 0.d0 ELSE ! ENRICH DERIVATIVES AT TWO ENDS OF EDGE DH1 = DH (1, K) DH2 = DH (2, K) DH3 = DH (3, K) K1 = I1 (K) K2 = I2 (K) DH (1, K1) = DH (1, K1) - DH1 DH (2, K1) = DH (2, K1) - DH2 DH (3, K1) = DH (3, K1) - DH3 DH (1, K2) = DH (1, K2) - DH1 DH (2, K2) = DH (2, K2) - DH2 DH (3, K2) = DH (3, K2) - DH3 DH (1, K) = DH1 + DH1 DH (2, K) = DH2 + DH2 DH (3, K) = DH3 + DH3 END IF END DO END SUBROUTINE DERIV_8_20_H SUBROUTINE DERIV_20_H (R, S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! LOCAL DERIVATIVES OF INTERPOLATION FUNCTIONS FOR AN ! 20 NODE HEXAHEDRON, SEE SHAPE_20_H FOR TOPOLOGY FIGURE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: DH (3, 20) INTEGER :: I1 (20), I2 (20), K, K1, K2 REAL(DP) :: FR, FS, FT, RP, RM, S_P, SM, TP, TM, & RZ, SZ, TZ, DH1, DH2, DH3, FP, FM DATA I1 / 8 * 0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4 / DATA I2 / 8 * 0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8 / DATA FP, FM / 0.5d0, -0.5d0 / ! FOR PARAMETER DEFINITIONS SEE SUBROUTINE SHAPE_20_H ! R,S,T = LOCAL COORDINATES OF THE POINT -1 <= (R,S,T) <= +1 ! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS ! I1, I2 = CORNER NODES OF TWELVE EDGES RP = 0.5d0 * (1.d0 + R) ; S_P = 0.5d0 * (1.d0 + S) ; TP = 0.5d0 * (1.d0 + T) RM = 0.5d0 * (1.d0 - R) ; SM = 0.5d0 * (1.d0 - S) ; TM = 0.5d0 * (1.d0 - T) RZ = 1.d0 - R * R ; SZ = 1.d0 - S * S ; TZ = 1.d0 - T * T FR = - 2.d0 * R ; FS = - 2.d0 * S ; FT = - 2.d0 * T ! DERIVATIVES OF TRI-LINEAR CORNERS DH (1, 1) = TP * S_P * FP DH (2, 1) = TP * FP * RP DH (3, 1) = FP * S_P * RP DH (1, 2) = TP * S_P * FM DH (2, 2) = TP * FP * RM DH (3, 2) = FP * S_P * RM DH (1, 3) = TP * SM * FM DH (2, 3) = TP * FM * RM DH (3, 3) = FP * SM * RM DH (1, 4) = TP * SM * FP DH (2, 4) = TP * FM * RP DH (3, 4) = FP * SM * RP DH (1, 5) = TM * S_P * FP DH (2, 5) = TM * FP * RP DH (3, 5) = FM * S_P * RP DH (1, 6) = TM * S_P * FM DH (2, 6) = TM * FP * RM DH (3, 6) = FM * S_P * RM DH (1, 7) = TM * SM * FM DH (2, 7) = TM * FM * RM DH (3, 7) = FM * SM * RM DH (1, 8) = TM * SM * FP DH (2, 8) = TM * FM * RP DH (3, 8) = FM * SM * RP ! DERIVATIVES OF EDGE BUBBLES DH (1, 9) = TP * S_P * FR * 0.5d0 DH (2, 9) = TP * FP * RZ * 0.5d0 DH (3, 9) = FP * S_P * RZ * 0.5d0 DH (1, 10) = TP * SZ * FM * 0.5d0 DH (2, 10) = TP * FS * RM * 0.5d0 DH (3, 10) = FP * SZ * RM * 0.5d0 DH (1, 11) = TP * SM * FR * 0.5d0 DH (2, 11) = TP * FM * RZ * 0.5d0 DH (3, 11) = FP * SM * RZ * 0.5d0 DH (1, 12) = TP * SZ * FP * 0.5d0 DH (2, 12) = TP * FS * RP * 0.5d0 DH (3, 12) = FP * SZ * RP * 0.5d0 DH (1, 13) = TM * S_P * FR * 0.5d0 DH (2, 13) = TM * FP * RZ * 0.5d0 DH (3, 13) = FM * S_P * RZ * 0.5d0 DH (1, 14) = TM * SZ * FM * 0.5d0 DH (2, 14) = TM * FS * RM * 0.5d0 DH (3, 14) = FM * SZ * RM * 0.5d0 DH (1, 15) = TM * SM * FR * 0.5d0 DH (2, 15) = TM * FM * RZ * 0.5d0 DH (3, 15) = FM * SM * RZ * 0.5d0 DH (1, 16) = TM * SZ * FP * 0.5d0 DH (2, 16) = TM * FS * RP * 0.5d0 DH (3, 16) = FM * SZ * RP * 0.5d0 DH (1, 17) = TZ * S_P * FP * 0.5d0 DH (2, 17) = TZ * FP * RP * 0.5d0 DH (3, 17) = FT * S_P * RP * 0.5d0 DH (1, 18) = TZ * S_P * FM * 0.5d0 DH (2, 18) = TZ * FP * RM * 0.5d0 DH (3, 18) = FT * S_P * RM * 0.5d0 DH (1, 19) = TZ * SM * FM * 0.5d0 DH (2, 19) = TZ * FM * RM * 0.5d0 DH (3, 19) = FT * SM * RM * 0.5d0 DH (1, 20) = TZ * SM * FP * 0.5d0 DH (2, 20) = TZ * FM * RP * 0.5d0 DH (3, 20) = FT * SM * RP * 0.5d0 ! LOOP OVER TWELVE EDGES DO K = 9, 20 ! ENRICH DERIVATIVES AT TWO ENDS OF EDGE DH1 = DH (1, K) DH2 = DH (2, K) DH3 = DH (3, K) K1 = I1 (K) K2 = I2 (K) DH (1, K1) = DH (1, K1) - DH1 DH (2, K1) = DH (2, K1) - DH2 DH (3, K1) = DH (3, K1) - DH3 DH (1, K2) = DH (1, K2) - DH1 DH (2, K2) = DH (2, K2) - DH2 DH (3, K2) = DH (3, K2) - DH3 DH (1, K) = DH1 + DH1 DH (2, K) = DH2 + DH2 DH (3, K) = DH3 + DH3 END DO END SUBROUTINE DERIV_20_H SUBROUTINE DERIV_2_L (R, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! DERIVATIVES OF A 2 NODE LINE ELEMENT ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R REAL(DP), INTENT(OUT) :: DH (2) ! R IS UNIT COORD. R=-1 1------------2 R=1 (unused) DH (1) = - 0.5d0 DH (2) = 0.5d0 END SUBROUTINE DERIV_2_L SUBROUTINE DERIV_3_T (S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES OF A THREE NODE UNIT TRIANGLE ! SEE SUBROUTINE SHAPE_3_T ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: S, T REAL(DP), INTENT(OUT) :: DH (2, 3) ! S,T = LOCAL COORDINATES OF THE POINT ! DH(1,K) = DH(K)/DS ! DH(2,K) = DH(K)/DT ! NODAL COORDS ARE : 1-(0,0) 2-(1,0) 3-(0,1) DH (1, 1) = - 1.d0 DH (1, 2) = 1.d0 DH (1, 3) = 0.d0 DH (2, 1) = - 1.d0 DH (2, 2) = 0.d0 DH (2, 3) = 1.d0 END SUBROUTINE DERIV_3_T !SUBROUTINE DERIV_4_12_Q (R, S, DH, ELEM_NODES) !! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* !! DERIVATIVES OF 4 TO 12 NODE QUADRILATERAL !! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* !Use Precision_Module ! IMPLICIT NONE ! REAL(DP), INTENT(IN) :: R, S ! INTEGER, INTENT(IN) :: ELEM_NODES (12) ! REAL(DP), INTENT(OUT) :: DH (2, 12) ! DH = 0.d0 ! STOP 'CONSTRUCT DERIV_4_12_Q FROM SHAPE_4_12_Q' !END SUBROUTINE DERIV_4_12_Q SUBROUTINE DERIV_4_L (X, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! CALCULATE DERIVATIVE OF A 4 NODE LINE ELEMENT ! IN NATURAL COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: DH (4) ! DH = LOCAL DERIVATIVE OF SHAPE_4_L ! X = LOCAL COORDIATE OF POINT, -1 <= X <= +1 ! LOCAL NODE COORD. ARE 1----2---3----4 DH (1) = (1.d0 + 18.d0 * X - 27.d0 * X * X) / 16.d0 DH (2) = 9.d0 / 16.d0 * ( - 3.d0 - 2.d0 * X + 9.d0 * X * X) DH (3) = 9.d0 / 16.d0 * (3.d0 - 2.d0 * X - 9.d0 * X * X) DH (4) = ( - 1.d0 + 18.d0 * X + 27.d0 * X * X) / 16.d0 END SUBROUTINE DERIV_4_L SUBROUTINE DERIV2_4_L (X, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! CALCULATE SECOND DERIVATIVE OF A 4 NODE LINE ELEMENT ! IN NATURAL COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: D2H (4) ! D2H = SECOND LOCAL DERIVATIVE OF SHAPE_4_L ! X = LOCAL COORDINATE OF POINT, -1 <= X <= +1 ! LOCAL NODE COORD. ARE 1----2---3----4 D2H (1) = (18.d0 - 54.d0 * X) / 16.d0 D2H (2) = 9.d0 / 16.d0 * ( - 2.d0 + 18.d0 * X) D2H (3) = 9.d0 / 16.d0 * ( - 2.d0 - 18.d0 * X) D2H (4) = (18.d0 + 54.d0 * X) / 16.d0 END SUBROUTINE DERIV2_4_L SUBROUTINE DERIV_4_P (R, S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVE OF A 4 NODE TETRAHEDRON (PYRAMID) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: DH (3, 4) DH (1, 1) = - 1.d0 DH (1, 2) = 1.d0 DH (1, 3) = 0.d0 DH (1, 4) = 0.d0 DH (2, 1) = - 1.d0 DH (2, 2) = 0.d0 DH (2, 3) = 1.d0 DH (2, 4) = 0.d0 DH (3, 1) = - 1.d0 DH (3, 2) = 0.d0 DH (3, 3) = 0.d0 DH (3, 4) = 1.d0 END SUBROUTINE DERIV_4_P SUBROUTINE DERIV_4_Q (R, S, DELTA) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES OF THE SHAPE FUNCTIONS FOR AN ! PARAMETRIC QUADRILATERAL WITH FOUR NODES ! SEE SHAPE_4_Q ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DELTA (2, 4) REAL(DP) :: RP, RM, S_P, SM ! DELTA(1,I) = DH/DR ! DELTA(2,I) = DH/DS ! H = LOCAL INTERPOLATION FUNCTIONS ! (R,S) = A POINT IN THE LOCAL COORDINATES ! HERE D(H(I))/DR = 0.25d0*R(I)*(1+S*S(I)), ETC. RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S DELTA (1, 1) = - 0.25d0 * SM DELTA (1, 2) = 0.25d0 * SM DELTA (1, 3) = 0.25d0 * S_P DELTA (1, 4) = - 0.25d0 * S_P DELTA (2, 1) = - 0.25d0 * RM DELTA (2, 2) = - 0.25d0 * RP DELTA (2, 3) = 0.25d0 * RP DELTA (2, 4) = 0.25d0 * RM END SUBROUTINE DERIV_4_Q SUBROUTINE DERIV2_4_Q (R, S, D2L_H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL SECOND DERIVATIVES OF THE SHAPE FUNCTIONS FOR ! A PARAMETRIC QUADRILATERAL WITH FOUR NODES ! H,RR H,RS H,SS ! SEE SHAPE_4_Q, DERIV_4_Q ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: D2L_H (3, 4) ! D2L_H(1,I) = D2 H/DR DR FOR NODE I ! D2L_H(2,I) = D2 H/DR DS ! D2L_H(3,I) = D2 H/DS DS ! H = LOCAL INTERPOLATION FUNCTIONS ! (R,S) = A POINT IN THE LOCAL COORDINATES ! HERE D(H(I))/DR = 0.25d0*R(I)*(1+S*S(I)), ETC. D2L_H (1, 1:4) = 0.d0 D2L_H (2, 1:4) = (/ -0.25d0, 0.25d0, 0.25d0, -0.25d0 /) D2L_H (3, 1:4) = 0.d0 END SUBROUTINE DERIV2_4_Q SUBROUTINE DERIV_4_R (R, S, A, B, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! FIRST DERIVATIVES OF A ! C0 RECTANGULAR ELEMENT IN UNIT COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, A, B REAL(DP), INTENT(OUT) :: DH (2, 4) ! DOF ARE W AT EACH NODE (N_G_DOF=1) ! X // R, Y // S. S ! A = PHYSICAL LENGTH IN X 4 -------- 3 ! B = PHYSICAL LENGTH IN Y I I ! R,S = LOCAL UNIT COORDS I I ! 1@(0,0), 3@(1,1) 1 -------- 2 ->R DH (1, 1) = (S - 1.d0) / A DH (1, 2) = (1.d0 - S) / A DH (1, 3) = S / A DH (1, 4) = - S / A DH (2, 1) = (R - 1.d0) / B DH (2, 2) = - R / B DH (2, 3) = R / B DH (2, 4) = (1.d0 - R) / B END SUBROUTINE DERIV_4_R SUBROUTINE DERIV_4_T (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES OF A ! UNIT TRIANGLE WITH CENTROID (4TH) NODE ! This is untested ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DH (2, 4) ! R,S = UNIT COORDINATES, SEE SHAPE_4_T ! DH = DERIVATIVES OF INTERPOLATION FUNCTIONS ! NODE 1@(0,0), 2@(1,0), 3@(0,1), 4@(1/3,1/3) DH (1, 4) = 27.d0 * S * (1.d0 - R - S) - 27.d0 * R * S DH (2, 4) = 27.d0 * R * (1.d0 - R - S) - 27.d0 * R * S DH (1, 1) = - 1.d0 - DH (1, 4) / 3.d0 DH (2, 1) = - 1.d0 - DH (2, 4) / 3.d0 DH (1, 2) = 1.d0 - DH (1, 4) / 3.d0 DH (2, 2) = - DH (2, 4) / 3.d0 DH (1, 3) = - DH (1, 4) / 3.d0 DH (2, 3) = 1.d0 - DH (2, 4) / 3.d0 END SUBROUTINE DERIV_4_T SUBROUTINE DERIV_7_T (S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES FOR A SIX NODE UNIT TRIANGLE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: S, T REAL(DP), INTENT(OUT) :: DH (2, 7) ! S,T = LOCAL COORDINATES, SEE SHAPE_7_T ! DH = LOCAL DERIVATIVES OF SHAPE FUNCTIONS ! DH(1,K) = DH(K)/DS, DH(2,K)=DH(K)/DT DH = 0.d0 STOP 'CONSTRUCT DERIV_7_T FROM SHAPE_7_T' END SUBROUTINE DERIV_7_T SUBROUTINE DERIV_8_H (R, S, T, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LOCAL DERIVATIVES FOR EIGHT NODE HEXAHEDRON ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: DH (3, 8) REAL(DP) :: RP, RM, S_P, SM, TP, TM ! R,S,T = LOCAL COORDINATES OF THE POINT ! DH(1,K)=DH/DR, DH(2,K)=DH/DS, DH(3,K)=DH/DT ! H = ELEMENT SHAPE FUNCTIONS, SEE SHAPE_8_H RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S TP = 1.d0 + T ; TM = 1.d0 - T DH (1, 1) = 0.125d0 * S_P * TP DH (2, 1) = 0.125d0 * RP * TP DH (3, 1) = 0.125d0 * RP * S_P DH (1, 2) = 0.125d0 * SM * TP DH (2, 2) = - 0.125d0 * RP * TP DH (3, 2) = 0.125d0 * RP * SM DH (1, 3) = 0.125d0 * SM * TM DH (2, 3) = - 0.125d0 * RP * TM DH (3, 3) = - 0.125d0 * RP * SM DH (1, 4) = 0.125d0 * S_P * TM DH (2, 4) = 0.125d0 * RP * TM DH (3, 4) = - 0.125d0 * RP * S_P DH (1, 5) = - 0.125d0 * S_P * TP DH (2, 5) = 0.125d0 * RM * TP DH (3, 5) = 0.125d0 * RM * S_P DH (1, 6) = - 0.125d0 * SM * TP DH (2, 6) = - 0.125d0 * RM * TP DH (3, 6) = 0.125d0 * RM * SM DH (1, 7) = - 0.125d0 * SM * TM DH (2, 7) = - 0.125d0 * RM * TM DH (3, 7) = - 0.125d0 * RM * SM DH (1, 8) = - 0.125d0 * S_P * TM DH (2, 8) = 0.125d0 * RM * TM DH (3, 8) = - 0.125d0 * RM * S_P END SUBROUTINE DERIV_8_H SUBROUTINE DERIV_9_Q (R, S, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! LOCAL DERIVATIVES FOR 9-NODED QUAD ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! SEE SHAPE_9_Q FOR TOPOLOGY Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: DH (2, 9) REAL(DP) :: RP, RM, S_P, SM, R2P1, R2M1, S2P1, S2M1 RM = R - 1.D0 ; SM = S - 1.D0 RP = R + 1.D0 ; S_P = S + 1.D0 S2P1 = S + S + 1.D0 ; S2M1 = S + S - 1.D0 R2P1 = R + R + 1.D0 ; R2M1 = R + R - 1.D0 DH (1, 1) = 0.25D0 * S * SM * R2M1 DH (1, 2) = 0.25D0 * S * SM * R2P1 DH (1, 3) = 0.25D0 * S * S_P * R2P1 DH (1, 4) = 0.25D0 * S * S_P * R2M1 DH (1, 5) = - S * SM * R DH (1, 6) = - 0.5D0 * S_P * SM * R2P1 DH (1, 7) = - S * S_P * R DH (1, 8) = - 0.5D0 * S_P * SM * R2M1 DH (1, 9) = 2.D0 * S_P * SM * R DH (2, 1) = 0.25D0 * S2M1 * R * RM DH (2, 2) = 0.25D0 * S2M1 * R * RP DH (2, 3) = 0.25D0 * S2P1 * R * RP DH (2, 4) = 0.25D0 * S2P1 * R * RM DH (2, 5) = - 0.5D0 * S2M1 * RP * RM DH (2, 6) = - S * R * RP DH (2, 7) = - 0.5D0 * S2P1 * RP * RM DH (2, 8) = - S * R * RM DH (2, 9) = 2.D0 * S * RP * RM END SUBROUTINE DERIV_9_Q SUBROUTINE DERIV2_9_Q (R, S, D2L_H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! LOCAL SECOND DERIVATIVES FOR 9-NODED QUAD ! H,RR H,RS H,SS ! SEE SHAPE_9_Q FOR TOPOLOGY ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: D2L_H (3, 9) REAL(DP) :: RP, RM, S_P, SM, R2P1, R2M1, S2P1, S2M1 ! D2L_H(1,I) = D2 H/DR DR FOR NODE I ! D2L_H(2,I) = D2 H/DR DS ! D2L_H(3,I) = D2 H/DS DS ! H = LOCAL INTERPOLATION FUNCTIONS ! (R,S) = A POINT IN THE LOCAL COORDINATES RM = R - 1.D0 ; SM = S - 1.D0 RP = R + 1.D0 ; S_P = S + 1.D0 S2P1 = S + S + 1.D0 ; S2M1 = S + S - 1.D0 R2P1 = R + R + 1.D0 ; R2M1 = R + R - 1.D0 D2L_H (1, 1) = 0.5D0 * S * SM D2L_H (1, 2) = 0.5D0 * S * SM D2L_H (1, 3) = 0.5D0 * S * S_P D2L_H (1, 4) = 0.5D0 * S * S_P D2L_H (1, 5) = - S * SM D2L_H (1, 6) = - S_P * SM D2L_H (1, 7) = - S * S_P D2L_H (1, 8) = - S_P * SM D2L_H (1, 9) = 2.D0 * S_P * SM D2L_H (2, 1) = 0.25D0 * S2M1 * R2M1 D2L_H (2, 2) = 0.25D0 * S2M1 * R2P1 D2L_H (2, 3) = 0.25D0 * S2P1 * R2P1 D2L_H (2, 4) = 0.25D0 * S2P1 * R2M1 D2L_H (2, 5) = - S2M1 * R D2L_H (2, 6) = - S * R2P1 D2L_H (2, 7) = - S2P1 * R D2L_H (2, 8) = - S * R2M1 D2L_H (2, 9) = 4.D0 * S * R D2L_H (3, 1) = 0.5D0 * R * RM D2L_H (3, 2) = 0.5D0 * R * RP D2L_H (3, 3) = 0.5D0 * R * RP D2L_H (3, 4) = 0.5D0 * R * RM D2L_H (3, 5) = - RP * RM D2L_H (3, 6) = - R * RP D2L_H (3, 7) = - RP * RM D2L_H (3, 8) = - R * RM D2L_H (3, 9) = 2.D0 * RP * RM END SUBROUTINE DERIV2_9_Q SUBROUTINE DERIV_C0_L (R, A, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIRST DERIVATIVE OF C0 LINE ELEMENT IN UNIT COORDINATES ! (SEE SHAPE_C0_L) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: DH (2) ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! DH = FIRST PHYSICAL DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE NODAL VALUES ONLY DH (1) = - 1.d0 / A DH (2) = 1.d0 / A END SUBROUTINE DERIV_C0_L SUBROUTINE DERIV_C1_L (R, A, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIRST DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES ! (A C1 ELEMENT, SEE SHAPE_C1_L) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: DH (4) ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! DH = FIRST PHYSICAL DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE DH (1) = 6.d0 * (R * R - R) / A DH (2) = 1.d0 - 4.d0 * R + 3.d0 * R * R DH (3) = 6.d0 * (R - R * R) / A DH (4) = 3.d0 * R * R - 2.d0 * R END SUBROUTINE DERIV_C1_L SUBROUTINE DERIV_C1_L2_NC (R, A, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIRST DERIVATIVES OF CUBIC HERMITE IN NATURAL COORDINATES ! (A C1 ELEMENT, SEE SHAPE_C1_L2_NC) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: DH (4) ! A = PHYSICAL LENGTH OF ELEMENT 1 ----+---- 2 --> R ! R = LOCAL COORDINATE OF POINT R=-1 0 R=1 ! DH = FIRST PHYSICAL DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 2/A * D()/DR ! DOF ARE FUNCTION AND SLOPE, D/DX, AT EACH NODE DH(1) = (-3 + 3*R*R) * 0.5d0 / A DH(2) = (-1 - 2*R + 3*R*R) * 0.25d0 DH(3) = ( 3 - 3*R*R) * 0.5d0 / A DH(4) = (-1 + 2*R + 3*R*R) * 0.25d0 END SUBROUTINE DERIV_C1_L2_NC SUBROUTINE DERIV2_C1_L (R, A, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! 2ND DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES ! (A C1 ELEMENT, SEE SHAPE_C1L & DERIV_C1_L) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: D2H (4) ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! D2H = SECOND DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE VALUE, SLOPE AT EACH END (WRT X) D2H (1) = 6.d0 * (R + R - 1.d0) / A**2 D2H (2) = ( - 4.d0 + 6.d0 * R) / A D2H (3) = 6.d0 * (1.d0 - R - R) / A**2 D2H (4) = (6.d0 * R - 2.d0) / A END SUBROUTINE DERIV2_C1_L SUBROUTINE DERIV2_C1_L2_NC (R, A, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! 2ND DERIVATIVES OF CUBIC HERMITE IN UNIT COORDINATES ! (A C1 ELEMENT, SEE SHAPE_C1_L2_NC & DERIV_C1_L2_NC) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: D2H (4) ! A = PHYSICAL LENGTH OF ELEMENT 1 ----+---- 2 --> R ! R = LOCAL COORDINATE OF POINT R=-1 0 R=1 ! D2H = SECOND DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE VALUE, SLOPE AT EACH END (WRT X) D2H(1) = ( 6*R) / (A*A) D2H(2) = (-2 + 6*R) * 0.5d0 / A D2H(3) = (-6*R) / (A*A) D2H(4) = ( 2 + 6*R) * 0.5d0 / A END SUBROUTINE DERIV2_C1_L2_NC SUBROUTINE DERIV_C2_L (R, A, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! FIRST DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT ! COORDINATES ( A C2 ELEMENT, SEE SURR SHAPE_C2_L ) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: DH (6) REAL(DP) :: R3 ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! DH = FIRST PHYSICAL DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE VALUE, SLOPE, CURVATURE AT EACH NODE (WRT X) R3 = R * R * R DH (1) = 30.d0 * (2.d0 * R3 - R * R - R3 * R) / A DH (2) = 1.d0 - 18.d0 * R * R + 32.d0 * R3 - 15.d0 * R3 * R DH (3) = 0.5d0*(2.d0 * R - 9.d0 * R * R + 12.d0 * R3 - 5.d0*R3*R)*A DH (4) = 30.d0 * R * R * (1.d0 - 2.d0 * R + R * R) / A DH (5) = 28.d0 * R3 - 15.d0 * R3 * R - 12.d0 * R * R DH (6) = 0.5d0 * R * R * (3.d0 - 8.d0 * R + 5.d0 * R * R) * A END SUBROUTINE DERIV_C2_L SUBROUTINE DERIV2_C2_L (R, A, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! 2ND DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT ! COORDINATES ( A C2 ELEMENT, SEE SUBR SHAPE_C2_L & DERIV_C2_L ) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: D2H (6) REAL(DP) :: R3 ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! D2H = SECOND DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE VALUE, SLOPE, 2ND DERIV AT EACH END (WRT X) R3 = R * R * R D2H (1) = 30.d0 * (6.d0 * R**2 - 2.d0 * R - 4.d0 * R3) / A**2 D2H (2) = ( - 36.d0 * R + 96.d0 * R * R - 60.d0 * R3) / A D2H (3) = 0.5d0 * (2.d0 - 18.d0 * R + 36.d0 * R * R - 20.d0 * R3) D2H (4) = 30.d0 * (2.d0 * R - 6.d0 * R * R + 4.d0 * R3) / A**2 D2H (5) = (84.d0 * R * R - 60.d0 * R3 - 24.d0 * R) / A D2H (6) = 0.5d0 * (6.d0 * R - 24.d0 * R * R + 20.d0 * R3) END SUBROUTINE DERIV2_C2_L SUBROUTINE DERIV3_C2_L (R, A, D3H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! 3RD DERIVATIVES OF FIFTH ORDER HERMITE ELEMENT IN UNIT ! COORDINATES (A C2 ELEMENT, SEE SUBR SHAPE_C2_L & DERIV_C2_L) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: D3H (6) ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! D3H = THIRD DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE VALUE, SLOPE, 2ND DERIV AT EACH END (WRT X) !D3H (1) = 30.d0 * (6.d0 * R**2 - 2.d0 * R - 4.d0 * R3) / A**2 D3H (1) = 30.d0 * (12.d0 * R - 2.d0 - 12.d0 * R * R) / A**3 !b !D3H (2) = ( - 36.d0 * R + 96.d0 * R * R - 60.d0 * R3) / A D3H (2) = ( - 36.d0 + 192.d0 * R - 180.d0 * R * R) / A**2 !b !D3H (3) = 0.5d0 * (2.d0 - 18.d0 * R + 36.d0 * R * R - 20.d0 * R3) D3H (3) = 0.5d0 * (-18.d0 + 72.d0 * R - 60.d0 * R * R) / A !b !D3H (4) = 30.d0 * (2.d0 * R - 6.d0 * R * R + 4.d0 * R3) / A**2 D3H (4) = 30.d0 * (2.d0 - 12.d0 * R + 12.d0 * R * R) / A**3 !b !D3H (5) = (84.d0 * R * R - 60.d0 * R3 - 24.d0 * R) / A D3H (5) = (168.d0 * R - 180.d0 * R * R - 24.d0) / A**2 !b !D3H (6) = 0.5d0 * (6.d0 * R - 24.d0 * R * R + 20.d0 * R3) D3H (6) = 0.5d0 * (6.d0 - 48.d0 * R + 60.d0 * R * R) / A !b END SUBROUTINE DERIV3_C2_L SUBROUTINE DERIV_C3_L (R, A, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIRST DERIVATIVES FOR 7TH ORDER HERMITE LINE ELEMENT ! ( A C3 ELEMENT, SEE SHAPE_C3_L ) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: DH (8) REAL(DP) :: R2, R3, R4 ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! DH = FIRST PHYSICAL DERIVATIVES OF H ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X) R2 = R * R R3 = R2 * R R4 = R2 * R2 DH (1) = 140.d0 * R3 * (R3 - 3.d0 * R2 + 3.d0 * R - 1.d0) / A DH (2) = 1.d0 - R3 * (80.d0 - 225.d0 * R + 216.d0 * R2 - 70.d0 * R3) DH (3) = R * (1.d0 - 20.d0*R2 + 50.d0*R3 - 45.d0*R4 + 14.d0*R2*R3) * A DH (4) = R2 * (3.d0 - 16.d0 * R + 30.d0 * R2 - 24.d0 * R3 & + 7.d0 * R4) * A * A / 6. DH (5) = 140.d0 * R3 * (1.d0 - 3.d0 * R + 3.d0 * R2 - R3) / A DH (6) = R3 * (70.d0 * R3 - 204.d0 * R2 + 195.d0 * R - 60.) DH (7) = R3 * (10.d0 - 35.d0 * R + 39.d0 * R2 - 14.d0 * R3) * A DH (8) = R3 * (7.d0 * R3 - 18.d0 * R2 + 15.d0 * R - 4.d0) * A * A / 6. END SUBROUTINE DERIV_C3_L SUBROUTINE DERIV2_C3_L (R, A, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! 2ND DERIV FOR SEVENTH ORDER HERMITE LINE ELEMENT ! ( A C3 ELEMENT, SEE SHAPE_C3_L & DERIV_C3_L ) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: D2H (8) REAL(DP) :: R2, R3, R4 ! A = PHYSICAL LENGTH OF ELEMENT 1 -------- 2 --> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! D()/DX = D()/DR DR/DX = 1/A * D()/DR ! D2H = SECOND DERIVATIVES OF SHAPE FUNCTIONS ARRAY ! DOF ARE VALUE, SLOPE, 2ND, 3RD DERIV AT EACH END (WRT X) R2 = R*R R3 = R2*R R4 = R2*R2 D2H (1) = 140.d0*(6.d0*R2*R3 - 15.d0*R4 + 12.d0*R3 - 3.d0*R2) / A**2 D2H (2) = - (240.d0*R2 - 900.d0*R3 + 1080.d0*R4 - 420.d0*R3*R2) / A D2H (3) = (1.d0 - 60.d0*R2 + 200.d0*R3 - 225.d0*R4 + 84.d0*R2*R3) D2H (4) = (R - 8.d0*R2 + 20.d0*R3 - 20.d0*R4 + 7.d0*R2*R3)*A D2H (5) = 140.d0*(3.d0*R2 - 12.d0*R3 + 15.d0*R4 - 6.d0*R2*R3) / A**2 D2H (6) = (420.d0*R3*R2 - 1020.d0*R4 + 780.d0*R3 - 180.d0*R2) / A D2H (7) = (30.d0*R2 - 140.d0*R3 + 195.d0*R4 - 84.d0*R3*R2) D2H (8) = (7.d0*R3*R2 - 15.d0*R4 + 10.d0*R3 - 2.d0*R2)*A END SUBROUTINE DERIV2_C3_L SUBROUTINE DERIV_CUBIC (B, A, DH) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! FIRST DERIVATIVES OF SHAPE FUNCTIONS FOR 1-D ! CUBIC HERMITE ELEMENT (A C1 ELEMENT) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: A, B REAL(DP), INTENT(OUT) :: DH (4) ! A = LENGTH OF ELEMENT (SEE SUBR SHAPE_CUBIC) ! B = COORDINATE OF POINT ! DH = FIRST DERIVATIVES OF H DH (1) = 6.d0 * (B * B - B) / A DH (2) = 1.d0 - 4.d0 * B + 3.d0 * B * B DH (3) = 6.d0 * (B - B * B) / A DH (4) = 3.d0 * B * B - 2.d0 * B END SUBROUTINE DERIV_CUBIC SUBROUTINE DERIV2_CUBIC (B, A, D2H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SECOND DERIVATIVES OF SHAPE FUNCTIONS FOR 1-D ! CUBIC HERMITE ELEMENT (A C1 ELEMENT) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: A, B REAL(DP), INTENT(OUT) :: D2H (4) ! A = LENGTH OF ELEMENT (SEE SUBR SHAPE_CUBIC) ! B = COORDINATE OF POINT ! D2H = SECOND DERIVATIVES OF H D2H (1) = (12.d0 * B - 6.d0) / A / A D2H (2) = (6.d0 * B - 4.d0) / A D2H (3) = (6.d0 - 12.d0 * B) / A / A D2H (4) = (6.d0 * B - 2.d0) / A END SUBROUTINE DERIV2_CUBIC SUBROUTINE DERIV_SHAPE_L_Q_H (NODEDG, LOCATE, NEDGE, N_SPACE, RST,& VALUE, DERIV) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * ! SHAPE FUNCTIONS AND DERIVATIVES FOR GENERAL SERENDIPITY ! LINE, QUAD, OR OR HEXAHEDRON WITH AN ! ARBITRARY NUMBER OF NODES ON EACH EDGE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* * Use Precision_Module IMPLICIT NONE INTEGER, PARAMETER :: MAXDEG = 20 REAL(DP), INTENT(IN) :: RST (3) INTEGER, INTENT(IN) :: NODEDG (12), LOCATE, NEDGE, N_SPACE REAL(DP), INTENT(OUT) :: DERIV (3), VALUE ! NODEDG = NUMBER ON NODES ON 1,4, OR 12 EDGES ! LOCAL = LOCAL COORDINATE PARALLEL TO EACH EDGE ! NEDGE = EDGE NUMBER OR CORNER NUMBER OF THE NODE COMPUTED ! N_SPACE = NUMBER OF SPATIAL DIMENSIONS ! RST = LOCAL COORDINATES FOR EVALUATION ! VALUE = SHAPE FUNCTION VALUE (RETURNED) ! DERIV = SHAPE FUNCTION DERIVATIVES (RETURNED) REAL(DP) :: BLKCRD (3, 8), POLI2 (3), CDRFN (3), FARSID (3), & CRDEDG (3, MAXDEG + 1) INTEGER :: NEATC (3, 8), NODEOP (2, 12), NODATC (3), LOCAL (12) INTEGER :: ICORD, NSIDE, INODE, NOPV1, NOPV2, ISRFN, & INOD1, INOD2, ICOR1, ICOR2 REAL(DP) :: POLI1, POLI3, CPNUL, PLAN2, DPLA2, DETP2, DETP3, & DPOL1, DPOL2, DPOL3 DATA BLKCRD / -1.d0, -1.d0, -1.d0, 1.d0, -1.d0, -1.d0, & 1.d0, 1.d0, -1.d0, -1.d0, 1.d0, -1.d0, & -1.d0, -1.d0, 1.d0, 1.d0, -1.d0, 1.d0, & 1.d0, 1.d0, 1.d0, -1.d0, 1.d0, 1.d0 / DATA NEATC /1, 4, 9, 1, 2, 10, 3, 2, 11, 3, 4, 12, & 5, 8, 9, 5, 6, 10, 7, 6, 11, 7, 8, 12 / DATA NODEOP / 7, 8, 8, 5, 5, 6, 6, 7, & 3, 4, 4, 1, 1, 2, 2, 3, & 3, 7, 4, 8, 1, 5, 2, 6 / DATA LOCAL / -1, -2, 1, 2, -1, -2, & 1, 2, -3, -3, -3, -3 / ! BLKCRD = BLOCK CORNER LOCAL COORDINATES ! CRDEDG = LOCAL COORDINATES OF SIDE NODES JOINING CORNER ! FARSID = FAR SIDE LOCAL COORDINATE ! LOCATE = POSITION NUMBER ON EDGE, 0 IF CORNER ! MAXDEG = MAXIMUM PLOYNOMIAL DEGREE ON ANY SIDE ! NEATC = THE 1, 2, OR 3 EDGES AT A CORNER ! NODATC = NUMBER OF SIDE NODES JOINING A CORNER ! NODEOP = 2 DIAGONALLY OPPOSITE NODES FOR EACH EDGE ! ! VALUE = A(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT ) ! DERIV = DA(R,S,T)*( P1(R) + P2(S) + P3(T) + CONSTANT ) ! + A(R,S,T)*( DP1(R) + DP2(S) + DP3(T) ) ! ! REF: G. ZAVARISE, ET AL, "AN ALGORITHM FOR GENERATION OF SHAPE ! FUNCTIONS IN SERENDIPITY ELEMENTS, I.J.N.M.E." ! ! T: S C8 *---E7----* C7 T: S 8---15----7 ! : / /. /: : / /. /: ! :/ / . / : :/ / 22 / : ! *---R / E12 / E11 *---R / . / 20 ! E8 . E6 : 16 21 / : ! / . / : / . / : ! / C4*.../.E3..* C3 / 4.13/.12..3 ! / . / / / . / / ! C5 *--E5-----* C6 / 5---------6 / ! : . : / : . : 11 ! : E4 : E2 : 14 19 / ! E9 . E10 / 17 . : 10 ! : . : / : . 18 / ! :. :/ :. :/ ! C1 *---E1----* C2 1----9----2 ! CORNER NODE & EDGE NUMBERS. 22 NODES: CORNERS, THEN BY EDGES. ! CCW IF |T|=1, ELSE IN POSITIVE T. ! === 3-D FORM === ! ! C4 *---E3----* C3 4----8----3 ! : . : :S : : ! : : : : : ! E4 E2 : 9 : ! : : *---R : : ! : : : : ! C1 *---E1----* C2 1--5-6-7--2 ! CORNER NODE & EDGE NUMBERS. 9 NODES: CORNERS, THEN BY EDGE ORDER. ! === 2-D FORM === ! ! C1 *---E1----* C2 1--3-4-5--2 ! CORNER NODE & EDGE NUMBERS. 5 NODES NUMBERED BY EDGE ORDER. ! === 1-D FORM === POLI1 = 1. IF (LOCATE == 0) THEN ! ! SHAPE FUNCTION FOR CORNER NODES ! DO ICORD = 1, N_SPACE POLI1 = POLI1*(RST (ICORD) + BLKCRD (ICORD, NEDGE)) & /(2*BLKCRD (ICORD, NEDGE)) END DO CPNUL = 1.d0 POLI2 (1) = 0.d0 POLI2 (2) = 0.d0 POLI2 (3) = 0.d0 DO ICORD = 1, N_SPACE NSIDE = NEATC (ICORD, NEDGE) NODATC (ICORD) = NODEDG (NSIDE) - 2 IF (NODATC (ICORD) > 0) THEN IF (NODATC (ICORD) > MAXDEG) STOP 'ERROR: MAXDEG, DERIV_SHAFN' CPNUL = CPNUL - 1. POLI2 (ICORD) = 1. FARSID (ICORD) = 2.d0/(NODEDG (NSIDE) - 1) DO INODE = 1, NODATC (ICORD) CRDEDG (ICORD, INODE) = - 1.d0 + FARSID (ICORD)*INODE POLI2 (ICORD) = POLI2 (ICORD)*(RST (ICORD) & - CRDEDG (ICORD, INODE))/(BLKCRD (ICORD, NEDGE) & - CRDEDG (ICORD, INODE)) END DO END IF END DO VALUE = POLI1*(POLI2 (1) + POLI2 (2) + POLI2 (3) + CPNUL) ELSE ! ! SHAPE FUNCTION FOR EDGE NODES ! NOPV1 = NODEOP (1, NEDGE) NOPV2 = NODEOP (2, NEDGE) ISRFN = ABS (LOCAL (NEDGE)) FARSID (1) = 2.d0/(NODEDG (NEDGE) - 1) CDRFN (1) = - BLKCRD (1, NOPV1) CDRFN (2) = - BLKCRD (2, NOPV1) CDRFN (3) = - BLKCRD (3, NOPV1) CDRFN (ISRFN) = (1.d0 - FARSID (1)*LOCATE)*LOCAL (NEDGE)/ISRFN DO ICORD = 1, N_SPACE POLI1 = POLI1*(RST (ICORD) - BLKCRD (ICORD, NOPV1)) & /(CDRFN (ICORD) - BLKCRD (ICORD, NOPV1)) END DO PLAN2 = (RST (ISRFN) - BLKCRD (ISRFN, NOPV2)) & /(CDRFN (ISRFN) - BLKCRD (ISRFN, NOPV2)) POLI3 = 1. NODATC (1) = NODEDG (NEDGE) - 2 IF (NODATC (1) > 0) THEN IF (NODATC (1) > MAXDEG) STOP 'MAXDEG, DERIV_SHAFN' DO INODE = 1, NODATC (1) CRDEDG (1, INODE) = - 1.d0 + FARSID (1)*INODE IF (ABS (CRDEDG (1, INODE) - CDRFN (ISRFN)) > 0.0001) THEN POLI3 = POLI3*(RST (ISRFN) - CRDEDG (1, INODE)) & /(CDRFN (ISRFN) - CRDEDG (1, INODE)) END IF END DO END IF VALUE = POLI1*PLAN2*POLI3 END IF ! ! DERIVATIVES OF SHAPE FUNCTIONS ! DO ICOR1 = 1, N_SPACE IF (LOCATE == 0) THEN ! ! DERIVATIVES FOR CORNER NODES ! DPOL1 = POLI2 (1) + POLI2 (2) + POLI2 (3) + CPNUL DO ICOR2 = 1, N_SPACE IF (ICOR2 /= ICOR1) THEN DPOL1 = DPOL1*(RST (ICOR2) + BLKCRD (ICOR2, NEDGE)) & /(2*BLKCRD (ICOR2, NEDGE)) ELSE DPOL1 = DPOL1/(2*BLKCRD (ICOR2, NEDGE)) END IF END DO DPOL2 = 0.d0 DO INOD1 = 1, NODATC (ICOR1) DETP2 = 1. DO INOD2 = 1, NODATC (ICOR1) IF (INOD2 /= INOD1) THEN DETP2 = DETP2*(RST (ICOR1) - CRDEDG (ICOR1, INOD2)) & /(BLKCRD (ICOR1, NEDGE) - CRDEDG (ICOR1, INOD2)) ELSE DETP2 = DETP2/(BLKCRD (ICOR1, NEDGE) - CRDEDG (ICOR1,& INOD2)) END IF END DO DPOL2 = DPOL2 + DETP2 END DO DPOL2 = DPOL2*POLI1 DERIV (ICOR1) = DPOL1 + DPOL2 ELSE ! ! DERIVATIVES FOR EDGE NODES ! DPOL1 = POLI3*PLAN2 DO ICOR2 = 1, N_SPACE IF (ICOR2 /= ICOR1) THEN DPOL1 = DPOL1*(RST (ICOR2) - BLKCRD (ICOR2, NOPV1)) & /(CDRFN (ICOR2) - BLKCRD (ICOR2, NOPV1)) ELSE DPOL1 = DPOL1/(CDRFN (ICOR2) - BLKCRD (ICOR2, NOPV1)) END IF END DO DPLA2 = 0.d0 DPOL3 = 0.d0 IF (ICOR1 == ISRFN) THEN DPLA2 = POLI1*POLI3/(CDRFN (ISRFN) - BLKCRD (ISRFN, & NOPV2)) DO INOD1 = 1, NODATC (1) IF (ABS (CRDEDG (1, INOD1) - CDRFN (ISRFN)) > 0.0001) & THEN DETP3 = 1. DO INOD2 = 1, NODATC (1) IF (ABS (CRDEDG (1, INOD2) - CDRFN (ISRFN)) & > 0.0001) THEN IF (INOD2 /= INOD1) THEN DETP3 = DETP3*(RST (ISRFN) - CRDEDG (1, INOD2))& /(CDRFN (ISRFN) - CRDEDG (1, INOD2)) ELSE DETP3 = DETP3/(CDRFN (ISRFN) - CRDEDG (1, INOD2)) END IF END IF END DO DPOL3 = DPOL3 + DETP3 END IF END DO DPOL3 = DPOL3*POLI1*PLAN2 END IF DERIV (ICOR1) = DPOL1 + DPLA2 + DPOL3 END IF END DO END SUBROUTINE DERIV_SHAPE_L_Q_H SUBROUTINE SHAPE_12_Q (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS OF A 12 NODE QUADRILATERAL ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (12) REAL(DP), PARAMETER :: F1 = 1.d0/32, F2 = 9.d0/32 REAL(DP) :: R_P, R_M, S_P, S_M, R2, S2 ! R,S = LOCAL COORDS OF PT 4--11---7---3 ! H = ELEM SHAPE FUNCTIONS I I ! 8 S 10 ! ELEMENT SKETCH TO RIGHT I .R I ! LOCAL COORD OF NODES: 12 6 ! 1-(-1,-1) 3-(+1,+1) I I ! SIDES OF ORDER 1, 2, OR 3 1---5---9---2 R_P = 1 + R ; R_M = 1 - R ; S_P = 1 + S ; S_M = 1 - S R2 = R*R ; S2 = S*S H ( 1) = F1*R_M*S_M*(-10 + 9*(R2 + S2)) H ( 2) = F1*R_P*S_M*(-10 + 9*(R2 + S2)) H ( 3) = F1*R_P*S_P*(-10 + 9*(R2 + S2)) H ( 4) = F1*R_M*S_P*(-10 + 9*(R2 + S2)) H ( 5) = F2*(1 - R2)*S_M*(1 - 3*R) H ( 6) = F2*(1 - S2)*R_P*(1 - 3*S) H ( 7) = F2*(1 - R2)*S_P*(1 + 3*R) H ( 8) = F2*(1 - S2)*R_M*(1 + 3*S) H ( 9) = F2*(1 - R2)*S_M*(1 + 3*R) H (10) = F2*(1 - S2)*R_P*(1 + 3*S) H (11) = F2*(1 - R2)*S_P*(1 - 3*R) H (12) = F2*(1 - S2)*R_M*(1 - 3*S) END SUBROUTINE SHAPE_12_Q SUBROUTINE SHAPE_15_T (R, S, H) ! ****************************************************** ! SHAPE FUNCTIONS FOR 15-NODED TRIANGLE ! ****************************************************** Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: PT1 = 42.666666666666667d0 REAL(DP), PARAMETER :: PT2 = 10.666666666666667d0 REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (15) REAL(DP) :: CC, T1, T2, T3, T4, T5, T6, T7, T8, T9 ! S ! | ! 3 ! R,S = LOCAL COORDS OF PT | \ ! H = ELEM SHAPE FUNCTIONS | \ ! 6 11 ! ELEM_NODES = ELEM INCIDENCES LIST | \ ! ELEMENT SKETCH TO RIGHT 9 15 8 ! | \ ! 1@(0,0) 2@(1,0) 2@(0,1) 12 13 14 5 ! | \ ! 1---4--7---10---2 --R CC = 1.d0 - R - S T1 = R - 0.25d0 ; T2 = R - 0.5d0 ; T3 = R - 0.75d0 T4 = S - 0.25d0 ; T5 = S - 0.5d0 ; T6 = S - 0.75d0 T7 = CC - 0.25d0 ; T8 = CC - 0.5d0 ; T9 = CC - 0.75d0 H (1) = PT2*CC*T7*T8*T9 H (2) = PT2*R*T1*T2*T3 H (3) = PT2*S*T4*T5*T6 H (4) = PT1*CC*R*T7*T8 H (5) = PT1*R*S*T1*T2 H (6) = PT1*S*CC*T4*T5 H (7) = 64.d0*CC*R*T1*T7 H (8) = 64.d0*R*S*T1*T4 H (9) = 64.d0*S*CC*T4*T7 H (10) = PT1*CC*R*T1*T2 H (11) = PT1*R*S*T4*T5 H (12) = PT1*S*CC*T7*T8 H (13) = 128.d0*R*S*CC*T7 H (14) = 128.d0*R*S*T1*CC H (15) = 128.d0*R*S*CC*T4 END SUBROUTINE SHAPE_15_T SUBROUTINE SHAPE_16_Q (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! LAGRANGIAN 16 NODE QUAD IN NATURAL COORDINATES ! USING TENSOR PRODUCTS OF 1D BASIS ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (16) REAL(DP) :: HR (4), HS (4) ! 4 4--11---7---3 The active Lagrangian Q16 elements ! . I S I have nodes input in the order shown ! 3 8 16 . 15 10 to the left. ! . I ..R I ! 2 12 13 14 6 ! . I I ! 1 1---5---9---2 ! ! 1---2---3---4 -->R ! Evaluate the 1D interpolations CALL SHAPE_4_L (R, HR) CALL SHAPE_4_L (S, HS) ! Form tensor products H (1) = HR (1)*HS (1) H (2) = HR (4)*HS (1) H (3) = HR (4)*HS (4) H (4) = HR (1)*HS (4) H (5) = HR (2)*HS (1) H (6) = HR (4)*HS (2) H (7) = HR (3)*HS (4) H (8) = HR (1)*HS (3) H (9) = HR (3)*HS (1) H (10) = HR (4)*HS (3) H (11) = HR (2)*HS (4) H (12) = HR (1)*HS (2) H (13) = HR (2)*HS (2) H (14) = HR (3)*HS (2) H (15) = HR (3)*HS (3) H (16) = HR (2)*HS (3) END SUBROUTINE SHAPE_16_Q SUBROUTINE SHAPE_16_QS (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** ! SHAPE FUNCTIONS FOR SERENDIPITY QUAD WITH 16 NODES ! A BI-4TH ORDER ELEMENT ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: PT667 = 0.66666666666667d0 REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (16) REAL(DP) :: RR, SS, RS, RP, RM, S_P, SM ! R,S = LOCAL COORDS OF PT 4--15--11---7---3 ! H = ELEM SHAPE FUNCTIONS I I ! 8 S 14 ! I . I ! ELEMENT SKETCH TO RIGHT 12 +..R 10 ! I I ! 1@(-1,-1) 3@(+1,+1) 16 6 ! I I ! 1---5---9--13---2 RR = R*R ; SS = S*S ; RS = R*S RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S H (1) = RM*SM*(-R*(4.d0*RR - 1.d0) - S*(4.d0*SS - 1.d0) - 3.d0)/12.d0 H (5) = - PT667*R*SM*RM*RP*(1.d0 - 2.d0*R) H (9) = 0.5d0*RM*RP*(1.d0 - 4.d0*RR)*SM H (13) = PT667*R*SM*RM*RP*(1.d0 + 2.d0*R) H (2) = RP*SM*(R*(4.d0*RR - 1.d0) - S*(4.d0*SS - 1.d0) - 3.d0)/12.d0 H (6) = - PT667*S*RP*SM*S_P*(1.d0 - 2.d0*S) H (10) = 0.5d0*SM*S_P*(1.d0 - 4.d0*SS)*RP H (14) = PT667*S*RP*SM*S_P*(1.d0 + 2.d0*S) H (3) = RP*S_P*(R*(4.d0*RR - 1.d0) + S*(4.d0*SS - 1.d0) - 3.d0)/12.d0 H (7) = PT667*R*S_P*RM*RP*(1.d0 + 2.d0*R) H (11) = 0.5d0*RM*RP*(1.d0 - 4.d0*RR)*S_P H (15) = - PT667*R*S_P*RM*RP*(1.d0 - 2.d0*R) H (4) = RM*S_P*(-R*(4.d0*RR - 1.d0) + S*(4.d0*SS - 1.d0) - 3.d0)/12.d0 H (8) = PT667*S*SM*S_P*(1.d0 + 2.d0*S)*RM H (12) = 0.5d0*SM*S_P*(1.d0 - 4.d0*SS)*RM H (16) = - PT667*S*RM*SM*S_P*(1.d0 - 2.d0*S) END SUBROUTINE SHAPE_16_QS SUBROUTINE SHAPE_16_R (R, S, A, B, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! WARNING, this code is not fully checked. No known bugs. ! C1 RECTANGULAR ELEMENT IN UNIT COORDINATES ! USING TENSOR PRODUCTS OF 1D BASIS ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, A, B REAL(DP), INTENT(OUT) :: H (16) REAL(DP) :: HR (4), HS (4) ! DOF ARE W W,X W,Y W,XY AT EACH NODE (N_G_DOF=4) ! X // R, Y // S. S ! A = PHYSICAL LENGTH IN X 4 -------- 3 ! B = PHYSICAL LENGTH IN Y I I ! R,S = LOCAL UNIT COORDS I I ! 1@(0,0), 3@(1,1) 1 -------- 2 ->R ! Evaluate the 1D interpolations CALL SHAPE_C1_L (R, A, HR) CALL SHAPE_C1_L (S, B, HS) ! Form tensor products H (1) = HR (1) * HS (1) H (2) = HR (2) * HS (1) H (3) = HR (1) * HS (2) H (4) = HR (2) * HS (2) H (5) = HR (3) * HS (1) H (6) = HR (4) * HS (1) H (7) = HR (3) * HS (2) H (8) = HR (4) * HS (2) H (9) = HR (3) * HS (3) H (10) = HR (4) * HS (3) H (11) = HR (3) * HS (4) H (12) = HR (4) * HS (4) H (13) = HR (1) * HS (3) H (14) = HR (2) * HS (3) H (15) = HR (1) * HS (4) H (16) = HR (2) * HS (4) END SUBROUTINE SHAPE_16_R SUBROUTINE SHAPE_17_Q (R, S, H) ! ****************************************************************** ! SHAPE FUNCTIONS FOR A SERENDIPITY QUAD WITH 17 NODES ! ****************************************************************** Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: PT667 = 0.6666666666666667d0 REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (17) REAL(DP) :: RR, SS, RS, RP, RM, S_P, SM ! R,S = LOCAL COORDS OF PT 4--15--11---7---3 ! H = ELEM SHAPE FUNCTIONS I I ! 8 S 14 ! I I I ! ELEMENT SKETCH TO RIGHT 12 17-R 10 ! I I ! 1@(-1,-1) 3@(+1,+1) 16 6 ! 17@(0,0) I I ! 1---5---9--13---2 RR = R*R ; SS = S*S ; RS = R*S RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S H (1) = RM*SM*( - 4.d0*R*(RR - 1.d0) - 4.d0*S*(SS - 1.d0) + 3.d0*RS)/12.d0 H (5) = - PT667*R*SM*RM*RP*(1.d0 - 2.d0*R) H (9) = 0.5d0*RM*RP*( - S - 4.d0*RR)*SM H (13) = PT667*R*SM*RM*RP*(1.d0 + 2.d0*R) H (2) = RP*SM*(4.d0*R*(RR - 1.d0) - 4.d0*S*(SS - 1.d0) - 3.d0*RS)/12.d0 H (6) = - PT667*S*RP*SM*S_P*(1.d0 - 2.d0*S) H (10) = 0.5d0*SM*S_P*(R - 4.d0*SS)*RP H (14) = PT667*S*RP*SM*S_P*(1.d0 + 2.d0*S) H (3) = RP*S_P*(4.d0*R*(RR - 1.d0) + 4.d0*S*(SS - 1.d0) + 3.d0*RS)/12.d0 H (7) = PT667*R*S_P*RM*RP*(1.d0 + 2.d0*R) H (11) = 0.5d0*RM*RP*(S - 4.d0*RR)*S_P H (15) = - PT667*R*S_P*RM*RP*(1.d0 - 2.d0*R) H (4) = RM*S_P*( - 4.d0*R*(RR - 1.d0) + 4.d0*S*(SS - 1.d0) - 3.d0*RS)/12.d0 H (8) = PT667*S*SM*S_P*(1.d0 + 2.d0*S)*RM H (12) = 0.5d0*SM*S_P*( - R - 4.d0*SS)*RM H (16) = - PT667*S*RM*SM*S_P*(1.d0 - 2.d0*S) H (17) = RM*RP*SM*S_P END SUBROUTINE SHAPE_17_Q SUBROUTINE SHAPE_8_20_H (R, S, T, H, ELEM_NODES) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! ELEMENT INTERPOLATION FUNCTIONS FOR AN 8 TO 20 NODE HEXAHEDRON ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: H (20) INTEGER, INTENT(IN) :: ELEM_NODES (20) INTEGER :: I1 (20), I2 (20), K, K1, K2 REAL(DP) :: RP, RM, S_P, SM, TP, TM, RZ, SZ, TZ, HK DATA I1/8*0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4/ DATA I2/8*0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8/ ! R,S,T = LOCAL COORDINATES OF THE POINT -1 LE (R,S,T) LE +1 ! H = ELEMENT INTERPOLATION FUNCTIONS ! H(I) = 0 IF ELEM_NODES(I) = 0 ! ELEM_NODES = ARRAY OF ELEMENT INCIDENCES, ! IF ELEM_NODES(I)=0 THEN LOCAL NODE I IS NOT CONSIDERED IN ANALYSIS ! I1, I2 = CORNER NODES OF TWELVE EDGES ! ! A SKETCH OF THE LOCAL NODES... T 3 *----*----* 2 ! : / . 10 /: ! FACES DEFINED BY R.H.R. ABOUT : / . / : ! THE POSITIVE LOCAL AXES *---S / *19 / *18 ! +R,(8,5,1,4,16,17,12,20) / 11* . 9* : ! -R,(7,6,2,3,14,18,10,19) / / . / : ! +S,(6,2,1,5,18,9,17,13) R / 7*.../*....* 6 ! -S,(7,3,4,8,19,11,20,15) / 12 . /14 / ! +T,(3,4,1,2,11,12,9,10) 4 *----*----*1 / ! -T,(7,8,5,6,15,16,13,14) : . : / ! AND R FACE : *15 : *13 ! 7 V 7 V 20 * . 17* / ! 7 V 7 V : . : / ! T <<< S COORD2 <<<< COORD1 :. 16 :/ ! ARE THE FACE COORD. PERMUTATIONS 8 *----*----* 5 RP = 0.5d0*(1.d0 + R) ; S_P = 0.5d0*(1.d0 + S) ; TP = 0.5d0*(1.d0 + T) RM = 0.5d0*(1.d0 - R) ; SM = 0.5d0*(1.d0 - S) ; TM = 0.5d0*(1.d0 - T) RZ = 1.d0 - R*R ; SZ = 1.d0 - S*S ; TZ = 1.d0 - T*T H (1) = TP*S_P*RP H (2) = TP*S_P*RM H (3) = TP*SM*RM H (4) = TP*SM*RP H (5) = TM*S_P*RP H (6) = TM*S_P*RM H (7) = TM*SM*RM H (8) = TM*SM*RP ! QUADRATIC EDGE BUBBLES H (9) = TP*S_P*RZ*0.5d0 H (10) = TP*SZ*RM*0.5d0 H (11) = TP*SM*RZ*0.5d0 H (12) = TP*SZ*RP*0.5d0 H (13) = TM*S_P*RZ*0.5d0 H (14) = TM*SZ*RM*0.5d0 H (15) = TM*SM*RZ*0.5d0 H (16) = TM*SZ*RP*0.5d0 H (17) = TZ*S_P*RP*0.5d0 H (18) = TZ*S_P*RM*0.5d0 H (19) = TZ*SM*RM*0.5d0 H (20) = TZ*SM*RP*0.5d0 ! LOOP OVER TWELVE ELEMENT EDGES DO K = 9, 20 IF ( ELEM_NODES (K) == 0 ) THEN ! SET UNUSED EDGE BUBBLE TO ZERO H (K) = 0.d0 ELSE ! ENRICH THE TWO CORNERS ON THE EDGE HK = H (K) K1 = I1 (K) K2 = I2 (K) H (K1) = H (K1) - HK H (K2) = H (K2) - HK H (K) = HK + HK END IF END DO END SUBROUTINE SHAPE_8_20_H SUBROUTINE SHAPE_20_H (R, S, T, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! ELEMENT INTERPOLATION FUNCTIONS FOR A 20 NODE HEXAHEDRON ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: H (20) INTEGER :: I1 (20), I2 (20), K, K1, K2 REAL(DP) :: RP, RM, S_P, SM, TP, TM, RZ, SZ, TZ, HK DATA I1/8*0, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4/ DATA I2/8*0, 2, 3, 4, 1, 6, 7, 8, 5, 5, 6, 7, 8/ ! R,S,T = LOCAL COORDINATES OF THE POINT -1 LE (R,S,T) LE +1 ! H = ELEMENT INTERPOLATION FUNCTIONS ! I1, I2 = CORNER NODES OF TWELVE EDGES ! ! A SKETCH OF THE LOCAL NODES... T 3 *----*----* 2 ! : / . 10 /: ! FACES DEFINED BY R.H.R. ABOUT : / . / : ! THE POSITIVE LOCAL AXES *---S / *19 / *18 ! +R,(8,5,1,4,16,17,12,20) / 11* . 9* : ! -R,(7,6,2,3,14,18,10,19) / / . / : ! +S,(6,2,1,5,18,9,17,13) R / 7*.../*....* 6 ! -S,(7,3,4,8,19,11,20,15) / 12 . /14 / ! +T,(3,4,1,2,11,12,9,10) 4 *----*----*1 / ! -T,(7,8,5,6,15,16,13,14) : . : / ! AND R FACE : *15 : *13 ! 7 V 7 V 20 * . 17* / ! 7 V 7 V : . : / ! T <<< S COORD2 <<<< COORD1 :. 16 :/ ! ARE THE FACE COORD. PERMUTATIONS 8 *----*----* 5 RP = 0.5d0*(1.d0 + R) ; S_P = 0.5d0*(1.d0 + S) ; TP = 0.5d0*(1.d0 + T) RM = 0.5d0*(1.d0 - R) ; SM = 0.5d0*(1.d0 - S) ; TM = 0.5d0*(1.d0 - T) RZ = 1.d0 - R*R ; SZ = 1.d0 - S*S ; TZ = 1.d0 - T*T H (1) = TP*S_P*RP H (2) = TP*S_P*RM H (3) = TP*SM*RM H (4) = TP*SM*RP H (5) = TM*S_P*RP H (6) = TM*S_P*RM H (7) = TM*SM*RM H (8) = TM*SM*RP ! QUADRATIC EDGE BUBBLES H (9) = TP*S_P*RZ*0.5d0 H (10) = TP*SZ*RM*0.5d0 H (11) = TP*SM*RZ*0.5d0 H (12) = TP*SZ*RP*0.5d0 H (13) = TM*S_P*RZ*0.5d0 H (14) = TM*SZ*RM*0.5d0 H (15) = TM*SM*RZ*0.5d0 H (16) = TM*SZ*RP*0.5d0 H (17) = TZ*S_P*RP*0.5d0 H (18) = TZ*S_P*RM*0.5d0 H (19) = TZ*SM*RM*0.5d0 H (20) = TZ*SM*RP*0.5d0 ! LOOP OVER TWELVE ELEMENT EDGES DO K = 9, 20 ! ENRICH THE TWO CORNERS ON THE EDGE HK = H (K) K1 = I1 (K) K2 = I2 (K) H (K1) = H (K1) - HK H (K2) = H (K2) - HK H (K) = HK + HK END DO END SUBROUTINE SHAPE_20_H SUBROUTINE SHAPE_2_L (R, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** ! SHAPE FUNCTIONS OF A 2 NODE LINE ELEMENT ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R REAL(DP), INTENT(OUT) :: H (2) ! R IS NATURAL COORD. R=-1 1------------2 R=1 H (1) = 0.5d0*(1.d0 - R) H (2) = 0.5d0*(1.d0 + R) END SUBROUTINE SHAPE_2_L SUBROUTINE SHAPE_3_T (S, T, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS FOR A THREE NODE UNIT TRIANGLE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: S, T REAL(DP), INTENT(OUT) :: H (3) ! S,T = LOCAL COORDINATES OF THE POINT 3 T ! H = SHAPE FUNCTIONS . . . ! NODAL COORDS 1-(0,0) 2-(1,0) 3-(0,1) 1..2 0..S H (1) = 1.d0 - S - T H (2) = S H (3) = T END SUBROUTINE SHAPE_3_T !SUBROUTINE SHAPE_4_12_Q (R, S, H, ELEM_NODES) !! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* !! SHAPE FUNCTIONS OF 4 TO 12 NODE QUADRILATERAL !! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* !Use Precision_Module ! IMPLICIT NONE ! REAL(DP), INTENT(IN) :: R, S ! REAL(DP), INTENT(OUT) :: H (12) ! INTEGER, INTENT(IN) :: ELEM_NODES (12) ! INTEGER :: IP (4), JP (4), NEXT (4), I, K ! REAL(DP) :: P, Q, TEMP ! DATA IP, JP, NEXT/1, 0, - 1, 0, 0, 1, 0, - 1, 2, 3, 4, 1/ ! !! R,S = LOCAL COORDS OF PT 4--11---7---3 !! H = ELEM SHAPE FUNCTIONS : S : !! ELEM_NODES = ELEM INCIDENCES LIST 8 : 10 !! ELEMENT SKETCH TO RIGHT : *..R : !! LOCAL COORD OF NODES: 12 6 !! 1-(-1,-1) 3-(+1,+1) : : !! SIDES OF ORDER 1, 2, OR 3 1---5---9---2 ! !!--> GENERATE FOUR NODE BILINEAR QUADRILATERAL ! CALL SHAPE_4_Q (R, S, H) ! !!--> LOOP OVER SIDES ! DO I = 1, 4 ! H (I + 4) = 0.d0 ! H (I + 8) = 0.d0 !! IS SIDE HIGHER THAN LINEAR ? ! IF ( ELEM_NODES (I + 4) > 0 ) THEN ! K = NEXT (I) !! FIND PT RELATIVE TO SIDE (CYCLIC COORDINATES) ! P = R*IP (I) + S*JP (I) ! Q = - R*JP (I) + S*IP (I) ! TEMP = (1.d0 - Q)*0.5d0 !!--> IS SIDE QUADRATIC OR CUBIC ? ! IF ( ELEM_NODES (I + 8) > 0 ) THEN !! CUBIC ! H (I + 8) = TEMP*(1.d0 - 3.d0*P - P*P + 3.d0*P*P*P)*9.d0/16.d0 ! H (I + 4) = TEMP*(1.d0 + 3.d0*P - P*P - 3.d0*P*P*P)*9.d0/16.d0 !! CORRECT CORNER POINTS FOR CUBIC ! H (I) = H (I) - H (I + 4)*2.d0/3.d0 - H (I + 8)/3.d0 ! H (K) = H (K) - H (I + 8)*2.d0/3.d0 - H (I + 4)/3.d0 ! ELSE !! QUADRATIC ! H (I + 4) = TEMP*(1.d0 - P*P) !! CORRECT CORNER POINTS FOR QUADRATIC ! H (I) = H (I) - H (I + 4)*0.5d0 ! H (K) = H (K) - H (I + 4)*0.5d0 ! END IF ! END IF ! END DO !END SUBROUTINE SHAPE_4_12_Q SUBROUTINE SHAPE_4_L (X, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! CALCULATE SHAPE FUNCTION OF A 4 NODE LINE ELEMENT ! IN NATURAL COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: X REAL(DP), INTENT(OUT) :: H (4) ! H = ELEMENT SHAPE FUNCITON ! X = LOCAL COORDIATE OF POINT, -1 <= X <= +1 ! LOCAL NODE COORD. ARE 1----2---3----4 H (1) = (1.d0 - X)*(3.d0*X + 1.d0)*(3.d0*X - 1.d0)/16.d0 H (2) = 9.d0*(1.d0 + X)*(X - 1.d0)*(3.d0*X - 1.d0)/16.d0 H (3) = 9.d0*(1.d0 + X)*(1.d0 - X)*(3.d0*X + 1.d0)/16.d0 H (4) = (1.d0 + X)*(3.d0*X + 1.d0)*(3.d0*X - 1.d0)/16.d0 END SUBROUTINE SHAPE_4_L SUBROUTINE SHAPE_4_P (R, S, T, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS OF A 4 NODE TETRAHEDRON ELEMENT ! (PYRAMID) IN UNIT COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module ! s IMPLICIT NONE ! | REAL(DP), INTENT(IN) :: R, S, T ! 3 REAL(DP), INTENT(OUT) :: H (4) ! | . ! . 1---2 --r H (1) = 1.d0 - R - S - T ! / . H (2) = R ! 4 H (3) = S ! / H (4) = T ! t END SUBROUTINE SHAPE_4_P SUBROUTINE SHAPE_4_Q (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS OF A 4 NODE PARAMETRIC QUAD ! IN NATURAL COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (4) REAL(DP) :: RP, RM, S_P, SM ! (R,S) = A POINT IN THE NATURAL COORDS 4--3 ! H = LOCAL INTERPOLATION FUNCTIONS I I ! H(I) = 0.25d0*(1+R*R(I))*(1+S*S(I)) I I ! R(I) = LOCAL R-COORDINATE OF NODE I 1--2 ! LOCAL COORDS, 1=(-1,-1) 3=(+1,+1) RP = 1.d0 + R ; RM = 1.d0 - R S_P = 1.d0 + S ; SM = 1.d0 - S H (1) = 0.25d0*RM*SM H (2) = 0.25d0*RP*SM H (3) = 0.25d0*RP*S_P H (4) = 0.25d0*RM*S_P END SUBROUTINE SHAPE_4_Q SUBROUTINE SHAPE_4_R (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! C0 RECTANGULAR ELEMENT IN UNIT COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (4) ! DOF ARE W AT EACH NODE (N_G_DOF=1) ! X // R, Y // S. S ! A = PHYSICAL LENGTH IN X 4 -------- 3 ! B = PHYSICAL LENGTH IN Y I I ! R,S = LOCAL UNIT COORDS I I ! 1@(0,0), 3@(1,1) 1 -------- 2 ->R H (1) = 1.d0 - R - S + R*S H (2) = R - R*S H (3) = R*S H (4) = S - R*S END SUBROUTINE SHAPE_4_R SUBROUTINE SHAPE_4_T (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** ! UNIT TRIANGLE WITH CENTROID (4TH) NODE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (4) ! R,S = UNIT COORDINATES, SEE SHAPE_4_T ! H = INTERPOLATION FUNCTIONS ! NODE 1@(0,0), 2@(1,0), 3@(0,1), 4@(1/3,1/3) H (4) = 27.d0*R*S*(1.d0 - R - S) H (1) = 1.d0 - R - S - H (4)/3.d0 H (2) = R - H (4)/3.d0 H (3) = S - H (4)/3.d0 END SUBROUTINE SHAPE_4_T SUBROUTINE SHAPE_7_T (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** ! QUADRATIC UNIT TRIANGLE PLUS CENTROID NODE ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: PT44 = 0.44444444444444d0 REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (7) ! R,S = UNIT COORDINATES, SEE SHAPE_6_T ! H = INTERPOLATION FUNCTIONS ! NODE 7@(1/3,1/3), OTHERS LIKE SHAPE_6_T H (7) = 27.d0*R*S*(1.d0 - R - S) H (1) = 1.d0 - 3.d0*R - 3.d0*S + 2.d0*R*R + 4.d0*R*S + 2.d0*S*S + H (7)/9.d0 H (2) = 2.d0*R*R - R + H (7)/9.d0 H (3) = 2.d0*S*S - S + H (7)/9.d0 H (4) = 4.d0*(R - R*R - R*S) - H (7)*PT44 H (5) = 4.d0*R*S - H (7)*PT44 H (6) = 4.d0*(S - R*S - S*S) - H (7)*PT44 END SUBROUTINE SHAPE_7_T SUBROUTINE SHAPE_8_H (R, S, T, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS OF 8 NODE PARAMETRIC HEXAHEDRON ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S, T REAL(DP), INTENT(OUT) :: H (8) REAL(DP) :: RP, RM, S_P, SM, TP, TM ! R,S,T = LOCAL COORDS OF PT ! R 3 *---------* 2 ! H = ELEM SHAPE FUNCTIONS ! | / . /| ! NODES ORDERED BY RHR ! | / . / | ! ABOUT THE R-AXIS ! *---T / . / | ! LOCAL COORD:1=(1,1,1) ! / / . / | ! 4=(1,1,-1) 7=(-1,-1,-1) ! / / . / | ! S / 7*.../.....* 6 RP = 1.d0 + R ! / . / / RM = 1.d0 - R ! 4 *---------*1 / S_P = 1.d0 + S ! | . | / SM = 1.d0 - S ! | . | / TP = 1.d0 + T ! | . | / TM = 1.d0 - T ! | . | / ! |. |/ ! 8 *---------* 5 H (1) = 0.125d0*RP*S_P*TP H (2) = 0.125d0*RP*SM*TP H (3) = 0.125d0*RP*SM*TM H (4) = 0.125d0*RP*S_P*TM H (5) = 0.125d0*RM*S_P*TP H (6) = 0.125d0*RM*SM*TP H (7) = 0.125d0*RM*SM*TM H (8) = 0.125d0*RM*S_P*TM END SUBROUTINE SHAPE_8_H SUBROUTINE SHAPE_9_Q (R, S, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS FOR 9-NODED QUAD ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, S REAL(DP), INTENT(OUT) :: H (9) REAL(DP) :: RP, RM, S_P, SM ! R,S = LOCAL COORDS OF PT 4-----7-----3 ! H = ELEM SHAPE FUNCTIONS I S I ! I . I ! ELEMENT SKETCH TO RIGHT 8 9..R 6 ! 1@(-1,-1) 3@(+1,+1) I I ! I I ! 1-----5-----2 RM = R - 1.D0 ; SM = S - 1.D0 RP = R + 1.D0 ; S_P = S + 1.D0 H (1) = 0.25D0*S*SM*R*RM H (2) = 0.25D0*S*SM*R*RP H (3) = 0.25D0*S*S_P*R*RP H (4) = 0.25D0*S*S_P*R*RM H (5) = - 0.5D0*S*SM*RP*RM H (6) = - 0.5D0*S_P*SM*R*RP H (7) = - 0.5D0*S*S_P*RP*RM H (8) = - 0.5D0*S_P*SM*R*RM H (9) = S_P*SM*RP*RM END SUBROUTINE SHAPE_9_Q SUBROUTINE SHAPE_C0_L (R, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** ! C0 LINE ELEM IN UNIT COORDINATES ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R REAL(DP), INTENT(OUT) :: H (2) ! H = SHAPE FUNCTIONS (R=0) 1----2 (R=1) ->R ! DOF ARE FUNCTION VALUES AT EACH NODE H (1) = 1.d0 - R H (2) = R END SUBROUTINE SHAPE_C0_L SUBROUTINE SHAPE_C1_L (R, A, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS FOR CUBIC HERMITE IN UNIT COORDINATES ! ( A C1 ELEMENT ) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: H (4) ! A = PHYSICAL LENGTH OF ELEMENT 1----------2 ---> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! H = SHAPE FUNCTIONS ARRAY ! DOF ARE FUNCTION AND SLOPE, WRT X, AT EACH NODE (N_G_DOF=2) ! D()/DX = D()/DR DR/DX = 1/A * D()/DR H(1) = 1.d0 - 3.0*R*R + 2.0*R*R*R H(2) = (R - 2.0*R*R + R*R*R)*A H(3) = 3.0*R*R - 2.0*R*R*R H(4) = (R*R*R - R*R)*A END SUBROUTINE SHAPE_C1_L SUBROUTINE SHAPE_C1_L2_NC (R, A, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS FOR CUBIC HERMITE IN NATURAL COORDINATES ! ( A C1 ELEMENT ) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: H (4) ! A = PHYSICAL LENGTH OF ELEMENT 1-----+-----2 ---> R ! R = LOCAL COORDINATE OF POINT R=-1 0 R=1 ! H = SHAPE FUNCTIONS ARRAY ! DOF ARE FUNCTION AND SLOPE, WRT X, AT EACH NODE (N_G_DOF=2) ! D()/DX = D()/DR DR/DX = 2/A * D()/DR H(1) = (2 - 3*R + R*R*R)/4 H(2) = (1 - R - R*R + R*R*R)*A/8 H(3) = (2 + 3*R - R*R*R)/4 H(4) = (-1 - R + R*R + R*R*R)*A/8 END SUBROUTINE SHAPE_C1_L2_NC SUBROUTINE SHAPE_C2_L (R, A, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** ! SHAPE FUNCTIONS FOR FIFTH ORDER HERMITE LINE ELEMENT ! ( A C2 ELEMENT IN UNIT COORDINATES) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-** Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: R, A REAL(DP), INTENT(OUT) :: H (6) REAL(DP) :: R3 ! A = PHYSICAL LENGTH OF ELEMENT 1----------2 -> R ! R = LOCAL COORDINATE OF POINT R=0 R=1 ! H = SHAPE FUNCTIONS ARRAY ! DOF ARE VALUE, SLOPE, CURVATURE AT EACH END (WRT X) ! D()/DX = D()/DR DR/DX = 1/A*D()/DR R3 = R*R*R H (1) = 1.d0 - 10.d0*R3 + 15.d0*R3*R - 6.d0*R3*R*R H (2) = (R - 6.d0*R3 + 8.d0*R3*R - 3.d0*R3*R*R)*A H (3) = (R*R - 3.d0*R3 + 3.d0*R3*R - R3*R*R)*A*A/2.d0 H (4) = 10.d0*R3 - 15.d0*R3*R + 6.d0*R3*R*R H (5) = (7.d0*R3*R - 3.d0*R3*R*R - 4.d0*R3)*A H (6) = (R3 - 2.d0*R3*R + R3*R*R)*A*A/2.d0 END SUBROUTINE SHAPE_C2_L SUBROUTINE SHAPE_C3_L (R, A, H) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* ! SHAPE FUNCTIONS FOR SEVENTH ORDER HERMITE LINE ELEMENT ! ( A C3 ELEMENT IN UNIT COORDINATES) ! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* Use Precision_Module IMPLICIT NONE RE