! copyright 2005, J. E. Akin, all rights reserved. SUBROUTINE BARACENTRIC_TO_UNIT_DERIV (N_GEOM, N_PARM, DLB, DLU) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * ! CONVERT SIMPLEX BARACENTRIC LOCAL DERIVATIVES TO ! UNIT COORDINATE LOCAL DERIVATIVES ! * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module ! for DP IMPLICIT NONE INTEGER, INTENT(IN) :: N_GEOM, N_PARM REAL(DP), INTENT(IN) :: DLB (N_PARM+1, N_GEOM) REAL(DP), INTENT(OUT):: DLU (N_PARM, N_GEOM) INTEGER :: J, K REAL(DP) :: DLB_LAST ! DLB = LOCAL DERIVATIVES IN BARACENTRIC COORDINATES ! DLU = LOCAL DERIVATIVES IN UNIT COORDINATES ! N_GEOM = NUMBER OF NODES ! N_PARM = DIMENSION OF UNIT SIMPLEX SPACE ! LOOP OVER THE INTERPOLATION FUNCTIONS DO K = 1, N_GEOM DLB_LAST = DLB (N_PARM+1, K) ! SUBTRACT OFF THE LAST BARACENTRIC VALUE DO J = 1, N_PARM DLU (J, K) = DLB (J, K) - DLB_LAST END DO END DO END SUBROUTINE BARACENTRIC_TO_UNIT_DERIV SUBROUTINE FUND_MAG_ONE (N_GEOM, DELTA, COORD, RJ, A, FFM) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND 1ST FUNDAMENTAL MAGNITUDES ON A 3-D SURFACE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: N_GEOM REAL(DP), INTENT(IN) :: DELTA (2, N_GEOM) REAL(DP), INTENT(IN) :: COORD (N_GEOM, 3) REAL(DP), INTENT(OUT) :: A (2, 2), RJ (2, 3), FFM ! N_GEOM = NUMBER OF NODES DEFINING THE SURFACE ! DELTA = LOCAL DERIVATIVES OF THE SHAPE FUNCTIONS ! COORD = THREE SPATIAL COORDINATES OF THE NODES ! RJ = REDUCED JACOBIAN MATRIX ! A = 1ST FUND. MAG. (METRIC TENSOR) ! FFM = DET(A) ! INDEX NOTATION: I=1,2,3, R,S=1,2, _=SUBSCRIPTS ! A_RS = X_I,R * X_I,S RJ_RI = X_I,R ! FORM REDUCED JACOBIAN RJ = MATMUL (DELTA, COORD) ! FORM METRIC TENSOR A (1, 1) = RJ (1, 1) * RJ (1, 1) + RJ (1, 2) * RJ (1, 2) & + RJ (1, 3) * RJ (1, 3) A (1, 2) = RJ (1, 1) * RJ (2, 1) + RJ (1, 2) * RJ (2, 2) & + RJ (1, 3) * RJ (2, 3) A (2, 2) = RJ (2, 1) * RJ (2, 1) + RJ (2, 2) * RJ (2, 2) & + RJ (2, 3) * RJ (2, 3) FFM = A (1, 1) * A (2, 2) - A (1, 2) * A (1, 2) ! SURFACE AREA DA = SQRT(FFM)*DR*DS ! ANGLE R-S COS (ANG) = A12/SQRT( A11*A22 ) ! LENGTHS DL1 = SQRT(A11)*DR ; DL2 = SQRT(A22)*DS IF ( FFM < 0.0 ) PRINT *,'FATAL DATA ERROR IN FUND_MAG_ONE' END SUBROUTINE FUND_MAG_ONE FUNCTION INSIDE_2D_POLYGON (N_VERT, XY, COORD) RESULT (IN_CODE) !* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND IF PT (X,Y) IS INSIDE ANY 2-D POLYGON WITH N_VERT VERTICIES !* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! IEEE TRAN EL COMP, EC-16, NO-6, DEC, 67 Use Precision_Module IMPLICIT NONE INTEGER, INTENT(IN) :: N_VERT ! number of vertices REAL(DP), INTENT(IN) :: XY (2) ! coordinates of point REAL(DP), INTENT(IN) :: COORD (N_VERT, 2) ! coords of vertices INTEGER :: IN_CODE ! return status REAL(DP), PARAMETER :: TWO_PI = 6.28318530717958623_DP REAL(DP), PARAMETER :: ZERO = 2*TINY (1.0_DP) REAL(DP) :: ALFA (N_VERT + 1), BETA (N_VERT + 1), DIFF ! work INTEGER :: I ! FOR RIGHT HAND SYSTEM ORDER VERTICIES IN CCW MANNER ! IN_CODE = 0 IF OUTSIDE, IN_CODE = 1 IF INSIDE, = ?? if on edge IN_CODE = 0 DO I = 1, N_VERT ALFA (I) = ATAN2 ((COORD (I, 2) - XY (2)), & (COORD (I, 1) - XY (1)))/TWO_PI END DO ALFA (N_VERT + 1) = ALFA (1) ; BETA (1) = ALFA (1) DO I = 2, N_VERT + 1 DIFF = ALFA (I) - ALFA (I-1) IF ( ABS (DIFF - 0.5_DP) <= ZERO ) THEN IN_CODE = 1 ; RETURN ELSE BETA (I) = BETA (I-1) + DIFF IF ( DIFF > 0.5_DP ) BETA (I) = BETA (I) - 1.0_DP IF ( DIFF < (-0.5_DP) ) BETA (I) = BETA (I) + 1.0_DP END IF END DO IF ( ABS (BETA (N_VERT + 1) - (BETA (1) + 1.0_DP)) <= ZERO ) & IN_CODE = 1 END FUNCTION INSIDE_2D_POLYGON SUBROUTINE INSIDE_3T (XY, COORD, INSIDE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! IS POINT XY INSIDE A TRIANGLE WITH GIVEN COORDINATES ! (NODES ARE ASSUMED TO BE NUMBERED COUNTERCLOCKWISE) ! (SEE INSIDE_2D_POLYGON AS ALTERNATE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: XY (2), COORD (3, 2) INTEGER, INTENT(OUT) :: INSIDE REAL(DP) :: TWO_A ! COORD = PHYSICAL COORDINATES OF NODES ! INSIDE = 0 IF XY IS NOT IN ELEMENT, ELSE = 1 ! XY = PHSYICAL COORDINATES OF POINT INSIDE = 0 ! INITALIZE AS OUTSIDE (ON EDGE IS INSIDE) ! TESTING TRIANGULAR AREA ON 3 SIDES (+ IF IN,0 ON EDGE) TWO_A = COORD (1, 1) * (COORD (2, 2) - XY (2) ) & + COORD (2, 1) * (XY (2) - COORD (1, 2) ) & + XY (1) * (COORD (1, 2) - COORD (2, 2) ) IF ( TWO_A < 0.0 ) RETURN TWO_A = COORD (2, 1) * (COORD (3, 2) - XY (2) ) & + COORD (3, 1) * (XY (2) - COORD (2, 2) ) & + XY (1) * (COORD (2, 2) - COORD (3, 2) ) IF ( TWO_A < 0.0 ) RETURN TWO_A = COORD (3, 1) * (COORD (1, 2) - XY (2) ) & + COORD (1, 1) * (XY (2) - COORD (3, 2) ) & + XY (1) * (COORD (3, 2) - COORD (1, 2) ) IF ( TWO_A < 0.0 ) RETURN ! POINT IS INSIDE INSIDE = 1 END SUBROUTINE INSIDE_3T SUBROUTINE IS_INSIDE (XYZ, INSIDE) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! is xyz inside a line, or flat T3 or Q4 element ? ! * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants Use Elem_Type_Data ! for GEOMETRY (LT_GEOM, N_SPACE) REAL(DP), INTENT(IN) :: XYZ (N_SPACE) ! point COORD INTEGER, INTENT(OUT) :: INSIDE ! 0=outside, 1=inside or on edge REAL(DP) :: Q_134 (3, n_space) ! 2nd triangle of quadrilateral ! GEOMETRY = COORDINATES OF ELEMENT GEOMETRY NODES ! INSIDE = FLAG FOR POINT LOCATION ! LT_GEOM = NUMBER OF GEOMETRICAL NODES ! N_SPACE = DIMENSION OF PHYSICAL SPACE ! XYZ = COORDINATES OF TEST POINT inside = 0 ! outside element ? select case (N_SPACE) case (1) if ( xyz(1) > minval( GEOMETRY (:,1) ) .and. & xyz(1) < maxval( GEOMETRY (:,1) ) ) inside = 1 case (2) ! see INSIDE_2D_POLYGON as alternative call INSIDE_3T (xyz, GEOMETRY(1:3,:), inside) ! check triangle if ( inside == 1 ) return ! was there if ( LT_N == 4 .or. LT_N == 8 .or. & ! check quadrilateral LT_N == 9 ) then ! also see WHERE_IN_4Q for alternative q_134 (1,:) = GEOMETRY (1,:) ! build other triangle q_134 (2,:) = GEOMETRY (3,:) q_134 (3,:) = GEOMETRY (4,:) call INSIDE_3T (xyz, q_134, inside) ! check other triangle end if case default stop 'Extend is_inside to 3-D' end select END SUBROUTINE IS_INSIDE SUBROUTINE LCONTR_ (R_IN_OUT, S_IN_OUT, V, COORD, XYZ, N_GEOM, & N_SPACE, KOUNT, L_TYPE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE COORDINATES OF POINTS ON A CONTOUR ! ON AN ISOPARAMETRIC SURFACE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, PARAMETER :: MAX = 300 REAL(DP), PARAMETER :: TRY_DL = 0.005 INTEGER, INTENT(IN) :: N_GEOM, N_SPACE, L_TYPE REAL(DP), INTENT(IN) :: V (N_GEOM) REAL(DP), INTENT(IN) :: COORD (N_GEOM, N_SPACE) INTEGER, INTENT(OUT) :: KOUNT REAL(DP), INTENT(OUT) :: XYZ (MAX, N_SPACE) REAL(DP), INTENT(INOUT) :: R_IN_OUT, S_IN_OUT INTEGER :: I REAL(DP) :: R (MAX), S (MAX), DER (2), PT(N_SPACE), P_XYZ (3) REAL(DP) :: H (N_GEOM), DH (2, N_GEOM) REAL(DP) :: DL, DV_DR, DV_DS, GRAD, R_NEW, S_NEW EQUIVALENCE (DER (1), DV_DR), (DER (2), DV_DS) ! MAX = MAX. NO. CONTOUR POINTS ! DL = LOCAL COORD. LENGTH OF EACH SEGMENT ! KOUNT = NO. OF PTS. ON CONTOUR CURVE ! R,S = LOCAL COORDINATES ! XYZ = GLOBAL COORD. ARRAY FOR COUNTOUR PTS. ! PT = LOCAL COORDINATES OF A SINGLE POINT ! P_XYZ = GLOBAL COORDINATES OF A SINGLE POINT ! COORD = GLOBAL COORD. OF NODES OF ELEMENT ! N_GEOM = NUMBER OF NODES PER ELEMENT ! N_SPACE = DIMENSION OF GLOBAL SPACE, 2 OR 3 ! R_IN_OUT = R COORD OF 1ST (IN) OR LAST (OUT) PT ! S_IN_OUT = S COORD OF 1ST (IN) OR LAST (OUT) PT ! H = ELEMENT INTERPOLATION FUNCTIONS ! DH = LOCAL COORD. DERIVATIVES OF H ! V = NODAL VALUES OF QUANTITY TO BE CONTOURED ! L_TYPE = ELEM. TYPE, 0=QUADRILATERAL, 1=TRIANGLE KOUNT = 1 DL = TRY_DL IF ( L_TYPE == 0 ) DL = DL * 2.0 !--> LOCAL COORDINATE CALCULATIONS R (1) = R_IN_OUT ; S (1) = S_IN_OUT ! MARCH ALONG CONTOUR 5 CONTINUE !9 use do while loop ! FORM SHAPE FUNCTION LOCAL DERIVATIVES PT (1) = R (KOUNT) ; PT (2) = S (KOUNT) CALL SCALAR_DERIVS (PT, DH) ! FORM LOCAL DERIVATIVE OF VARIABLE, DER = DH*V DER = MATMUL ( DH, V ) GRAD = SQRT (DV_DR * DV_DR + DV_DS * DV_DS) ! FIND LOCAL COORD. OF SEGMENT END R_NEW = R (KOUNT) - DL * DV_DS / GRAD S_NEW = S (KOUNT) + DL * DV_DR / GRAD ! IS NEXT POINT IN THE ELEMENT IF ( L_TYPE == 0 ) THEN ! QUADRILATERAL (-1 TO +1): IS POINT OUTSIDE ? 10 IF ( ABS (R_NEW) > 1.0 .OR. ABS (S_NEW) > 1.0 ) GO TO 20 ELSE ! TRIANGLE (UNIT COORDINATES): OUTSIDE ? IF ( R_NEW < 0.0 .OR. S_NEW < 0.0 .OR. & (R_NEW + S_NEW) > 1.0) GO TO 20 END IF ! ADD POINT TO CONTOUR LIST KOUNT = KOUNT + 1 R (KOUNT) = R_NEW S (KOUNT) = S_NEW IF ( KOUNT < MAX ) GO TO 5 !--> LINE COMPLETED IN LOCAL COORDINATES 20 CONTINUE !9 use do while loop IF ( KOUNT == 1 ) THEN DL = - DL GO TO 5 END IF ! RETURN LOCAL COORD. OF LAST POINT R_IN_OUT = R (KOUNT) S_IN_OUT = S (KOUNT) !--> CONVERT TO GLOBAL COORDINATES DO I = 1, KOUNT ! EVALUATE SHAPE FUNCTIONS CALL SCALAR_SHAPES (PT, H) ! EVALUATE GLOBAL COORDINATES, XYZ = H*COORD P_XYZ = MATMUL ( H, COORD ) XYZ (I, 1) = P_XYZ (1) XYZ (I, 2) = P_XYZ (2) IF ( N_SPACE == 3 ) XYZ (I, 3) = P_XYZ (3) END DO !--> CONTOUR FINISHED, RETURN FOR PLOTTING END SUBROUTINE LCONTR_ SUBROUTINE SURFACE_LAME (GRAD, ROOT_E, F, ROOT_G, H) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! GET THE SURFACE LAME PARAMETERS, E F G AND H, ! AT A POINT ON A PARAMETRIC CURVE OR SURFACE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: GRAD (N_PARM, N_GEOM) REAL(DP), INTENT(OUT) :: ROOT_E, F, ROOT_G, H REAL(DP) :: E, G ! NOD_PER_EL = NUMBER OF NODES PER ELEMENT ! N_PARM = NUMBER OF PARAMETRIC DIMENSIONS <= N_GEOM ! N_GEOM = DIMENSION OF PHYSICAL SPACE ! DELTA = LOCAL DERIVATIVES OF NOD_PER_EL INTERPOLATION FUNCTIONS ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! GRAD = TANGENT MATRIX FROM TANGENT_VECTOR; ! = JACOBIAN IF N_PARM = N_GEOM ! ROOT_E = SQUARE ROOT OF E, E = X,R^2 + Y,R^2 + Z,R^2 ! ROOT_F = SQUARE ROOT OF F, F = X,S*X,R + Y,S*Y,R + Z,S*Z,R ! ROOT_G = SQUARE ROOT OF G, G = X,S^2 + Y,S^2 + Z,S^2 IF ( N_PARM > 2) STOP 'INVALID el_parm ARGUMENT, SURFACE_LAME' G = 0. ; F = 0. ; E = GRAD (1, 1)**2 IF ( N_GEOM > 1) E = E + GRAD (1, 2)**2 IF ( N_GEOM == 3) E = E + GRAD (1, 3)**2 IF ( N_PARM > 1) THEN G = GRAD (2, 1) **2 + GRAD (2, 2) **2 IF ( N_GEOM == 3) G = G + GRAD (2, 3) **2 F = GRAD (1, 1) * GRAD (2, 1) + GRAD (1, 2) * GRAD (2, 2) IF ( N_GEOM == 3) F = F + GRAD (1, 3) * GRAD (2, 3) END IF H = SQRT (E * G - F * F) ; ROOT_G = SQRT (G) ; ROOT_E = SQRT (E) END SUBROUTINE SURFACE_LAME SUBROUTINE SURFACE_METRICS (P_GRAD, METRIC, E, F, G, H) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! COMPUTE SURFACE METRIC MEASURES & LAME PARAMETERS ! METRIC = P_GRAD*P_GRAD_TRANPOSE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: P_GRAD (N_PARM, N_SPACE) REAL(DP), INTENT(OUT) :: METRIC (N_PARM, N_PARM), E, F, G, H ! N_SPACE = Dimension of physical space, >= N_PARM ! N_PARM = Dimension of parametric space ! P_GRAD = Parametric gradient of physical coordinates ! METRIC = Surface differential geometry metric measure ! METRIC = [ E F ], Lame Parameters E, F, G, H ! [ F G ] H**2 = E*G - F**2 = Det|Metric| > 0 METRIC = MATMUL (P_GRAD, TRANSPOSE(P_GRAD)) E = METRIC (1, 1) F = METRIC (1, N_PARM) G = METRIC (N_PARM, N_PARM) H = SQRT (E * G - F**2) END SUBROUTINE SURFACE_METRICS SUBROUTINE SURF_NORM_VEC (RJ, FFM, VECTOR) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND UNIT NORMAL VECTOR TO A 3-D SURFACE AT A POINT ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE REAL(DP), INTENT(IN) :: RJ (2, 3), FFM REAL(DP), INTENT(OUT) :: VECTOR (3) REAL(DP) :: ROOT ! VECTOR = GLOBAL COMPONENTS OF UNIT NORMAL VECTOR ! THE FOLLOWING QUANTITIES ARE OBTAINED IN SUBR FUND_MAG_ONE ! RJ = REDUCED JACOBIAN MATRIX ! A = 1ST FUND MAGNITUDES(METRIC TENSOR), FFM = DET(A) ROOT = SQRT (FFM) VECTOR (1) = (RJ (1, 2)*RJ (2, 3) - RJ (1, 3)*RJ (2, 2))/ ROOT VECTOR (2) = (RJ (1, 3)*RJ (2, 1) - RJ (1, 1)*RJ (2, 3))/ ROOT VECTOR (3) = (RJ (1, 1)*RJ (2, 2) - RJ (1, 2)*RJ (2, 1))/ ROOT END SUBROUTINE SURF_NORM_VEC SUBROUTINE TANGENT_VECTOR (DELTA, COORD, GRAD) ! * * * * * * * * * * * * * * * * * * * * * * * * * ! CALCULATE THE TANGENT VECTORS AT A LOCAL POINT ! ON A PARAMETRIC CURVE OR PARAMETRIC SURFACE ! * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants IMPLICIT NONE REAL(DP), INTENT(IN) :: DELTA (N_PARM, N_GEOM) REAL(DP), INTENT(IN) :: COORD (N_GEOM, N_SPACE) REAL(DP), INTENT(OUT) :: GRAD (N_PARM, N_SPACE) ! N_GEOM = NUMBER OF NODES PER ELEMENT ! N_PARM = NUMBER OF PARAMETRIC DIMENSIONS <= N_SPACE ! N_SPACE = DIMENSION OF PHYSICAL SPACE ! DELTA = LOCAL DERIVATIVES OF N_GEOM INTERPOLATION ! FUNCTIONS AT POINT OF INTEREST. ! COORD = SPATIAL COORDINATES OF ELEMENT'S NODES ! GRAD = TANGENT MATRIX = DELTA*COORD ! ROW 1 IS DR/DU, 2 DR/DV, 3 DR/DW ! FOR N_PARM = N_SPACE, GRAD = JACOBIAN GRAD = MATMUL ( DELTA, COORD ) END SUBROUTINE TANGENT_VECTOR SUBROUTINE WHERE_IN_4Q (XY, COORD, INSIDE, R, S) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! FIND LOCAL (R,S) COORDINATES OF POINT XY IN A Q4 ELEMENT ! REF: C. HUA, FE IN ANAL & DESIGN, 7, 159-166, 1990 ! (SEE INSIDE_2D_POLYGON AS ALTERNATE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use System_Constants !b !b Use Precision_Module IMPLICIT NONE REAL(DP), PARAMETER :: ZERO = 0.D0, HALF = 0.5D0, & ONE = 1.D0, FOUR = 4.D0, TOL = 1.D-7 REAL(DP), INTENT(IN) :: XY (2), COORD (4, 2) INTEGER, INTENT(OUT) :: INSIDE REAL(DP), INTENT(OUT) :: R, S REAL(DP) :: A, B, C, ROOT, R1, R2 REAL(DP) :: TOABC (3, 4), ABC (3, 2), A1, A2, B1, B2, C1, C2 REAL(DP) :: AB, AC, AD, BC, BD, DB, DC, D1, D2 ,TWO_A REAL(DP) :: ERROR_X, ERROR_Y ! NODES ARE ASSUMED TO BE NUMBERED COUNTERCLOCKWISE ! COORD = PHYSICAL COORDINATES OF NODES ! INSIDE = 0 IF XY IS NOT IN ELEMENT, ELSE = 1 ! R,S = NATURAL COORDINATES OF THE PT (-1 TO +1) ! XY = PHSYICAL COORDINATES OF POINT INSIDE ELEMENT ! DATA TOABC / 1.D0, -1.D0, -1.D0, -1.D0, 1.D0, -1.D0, & 1.D0, 1.D0, 1.D0, -1.D0, -1.D0, 1.D0 / EQUIVALENCE (ABC(1, 1), A1), (ABC(2, 1), B1), (ABC(3, 1), C1),& (ABC(1, 2), A2), (ABC(2, 2), B2), (ABC(3, 2), C2) ! INITALIZE AS OUTSIDE (ON EDGE IS INSIDE) INSIDE = 0 ; R = HUGE(R) ; S = HUGE(S) ! LOOP OVER 4 SIDES TESTING TRIANGULAR AREA (+ IF IN) TWO_A = COORD (1, 1) * (COORD (2, 2) - XY (2) ) & + COORD (2, 1) * (XY (2) - COORD (1, 2) ) & + XY (1) * (COORD (1, 2) - COORD (2, 2) ) IF ( TWO_A < -TOL) RETURN TWO_A = COORD (2, 1) * (COORD (3, 2) - XY (2) ) & + COORD (3, 1) * (XY (2) - COORD (2, 2) ) & + XY (1) * (COORD (2, 2) - COORD (3, 2) ) IF ( TWO_A < -TOL) RETURN TWO_A = COORD (3, 1) * (COORD (4, 2) - XY (2) ) & + COORD (4, 1) * (XY (2) - COORD (3, 2) ) & + XY (1) * (COORD (3, 2) - COORD (4, 2) ) IF ( TWO_A < -TOL) RETURN TWO_A = COORD (4, 1) * (COORD (1, 2) - XY (2) ) & + COORD (1, 1) * (XY (2) - COORD (4, 2) ) & + XY (1) * (COORD (4, 2) - COORD (1, 2) ) IF ( TWO_A < -TOL) RETURN ! POINT IS INSIDE, GET LOCAL COORDINATES INSIDE = 1 D1 = 4. * XY (1) - COORD (1, 1) - COORD (2, 1) & - COORD (3, 1) - COORD (4, 1) D2 = 4. * XY (2) - COORD (1, 2) - COORD (2, 2) & - COORD (3, 2) - COORD (4, 2) ABC = MATMUL ( TOABC, COORD ) AB = A1 * B2 - A2 * B1 ; AC = A1 * C2 - A2 * C1 AD = A1 * D2 - A2 * D1 ; DC = D1 * C2 - D2 * C1 ! CHECK CASES IF ( (A1 * A2 * AB * AC) /= ZERO .OR. & (A1 == ZERO .AND. (A2 * C1) /= ZERO) .OR. & (A2 == ZERO .AND. (A1 * B2) /= ZERO) ) THEN A = AB ; B = C1 * B2 - C2 * B1 - AD ; C = DC ROOT = SQRT (B * B - FOUR * A * C) R1 = HALF * (-B + ROOT) / A ; R2 = HALF * (-B - ROOT) / A R = R2 IF ( -ONE <= R1 .AND. ONE >= R1 ) R = R1 S = (AD-AB * R) / AC ELSEIF ( ( (A1 * A2) /= ZERO) .AND. (AB == ZERO) ) THEN R = A1 * DC / (B1 * AC + A1 * AD) ; S = AD / AC ELSEIF ( ( (A1 * A2) /= ZERO) .AND. (AC == ZERO) ) THEN DB = D1 * B2 - D2 * B1 R = AD / AB ; S = A1 * DB / (C1 * AB + A1 * AD) ELSE BC = B1 * C2 - B2 * C1 ; BD = B1 * D2 - B2 * D1 R = DC / (A1 * D2 + BC) ; S = BD / (A2 * D1 + BC) END IF ! CHECK FOR ERRORS ERROR_X = ABS (B1 * R + C1 * S - D1 + A1 * R * S) ERROR_Y = ABS (B2 * R + C2 * S - D2 + A2 * R * S) IF ( ERROR_X > TOL .OR. ERROR_Y > TOL ) WRITE (6, * ) & 'WARNING: WHERE_IN_4Q ERROR_S ARE ', ERROR_X, ERROR_Y N_WARN = N_WARN + 1 ! INCREMENT WARNING IF ( -ONE > R .OR. ONE < R ) STOP 'WHERE_IN_4Q ERROR IN R' IF ( -ONE > S .OR. ONE < S ) STOP 'WHERE_IN_4Q ERROR IN S' END SUBROUTINE WHERE_IN_4Q SUBROUTINE TRI_GEOM (CODE, COORD, VALUE) ! * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT GEOMETRIC PROPERTIES OF AN ARBITRARY TRIANGLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE INTEGER, INTENT (IN) :: CODE REAL(DP), INTENT (IN) :: COORD (3, 2) REAL(DP), INTENT (OUT) :: VALUE REAL(DP), PARAMETER :: PI = 3.1415926536d0 REAL(DP) :: XI, XJ, XK, XB, YI, YJ, YK, YB, AREA ! VALUE = INTEGRAL OF (X**M)*(Y**N)*DA (0<=M,N<=2) ! VALUE BY CODE: ! 1-AREA, 2-VOL REV ABOUT Y AXIS, ! 3-VOL REV ABOUT Y, 4-FIRST MOMENT ABOUT X, ! 5-FIRST MOMENT ABOUT Y, 6-SEC MOMENT ABOUT X, ! 7-SEC MOMENT WRT X-Y AXES, 8-SEC MOMENT ABOUT Y XI = COORD (1, 1) ; XJ = COORD (2, 1) ; XK = COORD (3, 1) YI = COORD (1, 2) ; YJ = COORD (2, 2) ; YK = COORD (3, 2) AREA = 0.5 * (XI * (YJ - YK) + XJ * (YK - YI) + XK * (YI - YJ) ) XB = (XI + XJ + XK) / 3.d0 ; YB = (YI + YJ + YK) / 3.d0 SELECT CASE (CODE) CASE (1) ! M=0, N=0, AREA VALUE = AREA CASE (2) ! M=0, N=0, VOLUME OF REVOLUTION ABOUT X-AXIS VALUE = 2.d0 * PI * AREA * YB IF ( (ABS (YI) + ABS (YJ) + ABS (YK) ) /= & ABS (YI + YJ + YK) ) WRITE (6, 5) ; 5 FORMAT & ('TRI_GEOM: INVALID DATA FOR VOLUME OF REVOLUTION') CASE (3) ! M=0, N=0, VOLUME OF REVOLUTION ABOUT Y-AXIS VALUE = 2.d0 * PI * AREA * XB IF ( (ABS (XI) + ABS (XJ) + ABS (XK) ) /= & ABS (XI + XJ + XK) ) WRITE (6, 5) CASE (4) ! M=0, N=1, FIRST MOMENT ABOUT X-AXIS VALUE = AREA * YB CASE (5) ! M=1, N=0, FIRST MOMENT ABOUT Y-AXIS VALUE = AREA * XB CASE (6) ! M=0, N=2, SECOND MOMENT ABOUT X-AXIS VALUE = AREA * (YI * YI + YJ * YJ + YK * YK & + 9.0 * YB * YB) / 12.d0 CASE (7) ! M=1, N=1, SECOND MOMENT W.R.T. X-Y AXES VALUE = AREA * (XI * YI + XJ * YJ + XK * YK & + 9.0 * XB * YB) / 12.d0 CASE (8) ! M=2, N=0, SECOND MOMENT ABOUT Y-AXIS VALUE = AREA * (XI * XI + XJ * XJ + XK * XK & + 9.0 * XB * XB) / 12.d0 CASE DEFAULT STOP 'TRI_GEOM: ARGUMENT OUT OF RANGE' END SELECT END SUBROUTINE TRI_GEOM SUBROUTINE TRI_CG_INERTIA (COORD, INERTIA) ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * ! EXACT CENTROIDAL INERTIA TENSOR OF AN ARBITRARY TRIANGLE ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * Use Precision_Module IMPLICIT NONE integer, save :: count = 0 REAL(DP), INTENT (IN) :: COORD (3, 2) REAL(DP), INTENT (OUT) :: INERTIA (2, 2) REAL(DP) :: XI, XJ, XK, XB, YI, YJ, YK, YB, AREA XI = COORD (1, 1) ; XJ = COORD (2, 1) ; XK = COORD (3, 1) YI = COORD (1, 2) ; YJ = COORD (2, 2) ; YK = COORD (3, 2) AREA = 0.5 * (XI * (YJ - YK) + XJ * (YK - YI) + XK * (YI - YJ) ) XB = (XI + XJ + XK) / 3.d0 ; YB = (YI + YJ + YK) / 3.d0 ! With respect to global origin INERTIA (1, 1) = AREA * (YI * YI + YJ * YJ + YK * YK & + 9.0 * YB * YB) / 12.d0 INERTIA (2, 1) = -AREA * (XI * YI + XJ * YJ + XK * YK & + 9.0 * XB * YB) / 12.d0 INERTIA (2, 2) = AREA * (XI * XI + XJ * XJ + XK * XK & + 9.0 * XB * XB) / 12.d0 if ( count == 0 ) then print *,'INERT ', AREA, XB, YB print *,'INERT ', INERTIA end if ! With respect to centroid INERTIA (1, 1) = INERTIA (1, 1) - AREA * YB * YB INERTIA (2, 1) = INERTIA (2, 1) + AREA * XB * YB INERTIA (2, 2) = INERTIA (2, 2) - AREA * XB * XB INERTIA (1, 2) = INERTIA (2, 1) if ( count == 0 ) then count = 1 print *,'CG ', INERTIA end if END SUBROUTINE TRI_CG_INERTIA