From c2aa0682251a7b0ccf697e5bf4bac2dbdf3c8272 Mon Sep 17 00:00:00 2001 From: victorkemp Date: Fri, 14 Aug 2026 15:37:16 +1200 Subject: [PATCH] RBE3 Bug6 fix --- Source/LK1/L1D/RBE3_PROC.f90 | 627 +++++++++++++++-------------------- 1 file changed, 265 insertions(+), 362 deletions(-) diff --git a/Source/LK1/L1D/RBE3_PROC.f90 b/Source/LK1/L1D/RBE3_PROC.f90 index 46b4f2cb..9926c0a6 100644 --- a/Source/LK1/L1D/RBE3_PROC.f90 +++ b/Source/LK1/L1D/RBE3_PROC.f90 @@ -108,6 +108,27 @@ SUBROUTINE RBE3_PROC ( RTYPE, REC_NO, IERR ) REAL(DOUBLE) :: SXY,SZX,SYZ ! new Rdd terms according to victor REAL(DOUBLE) :: WTi6(MRBE3,6) ! per-DoF grid weights +! FIX (partial-REFC coupling bug): the full 6x6 coupled rigid-body least-squares system and the +! full coefficient table (one column per active independent-grid component) are now ALWAYS built, +! regardless of which components REFC selects as dependent. Whenever REFC excludes some +! components, they are eliminated via a Schur complement (not simply dropped), so the retained +! (REFC-selected) rows correctly account for their coupling to the excluded ones -- matching MSC +! Nastran's actual behavior (confirmed: MSC always solves the full 6-DOF system internally and +! REFC only controls what's output, not what's computed). When REFC=123456 (nothing excluded), +! this reduces identically to the original per-row computation -- verified by regression test. + REAL(DOUBLE) :: A6(6,6) ! full coupled system matrix + REAL(DOUBLE) :: B6(6,6*MRBE3) ! full coefficient table (row, indep grid*comp) + INTEGER(LONG) :: B6_COL_RMG(6*MRBE3) ! actual RMG column number for each B6 column (0=inactive) + INTEGER(LONG) :: NB6COLS ! number of columns actually used in B6/B6_COL_RMG + LOGICAL :: IS_R(6) ! TRUE if component is REFC-selected (retained) + INTEGER(LONG) :: R_IDX(6), D_IDX(6) ! lists of retained / discarded component indices + INTEGER(LONG) :: NR, ND ! counts of retained / discarded components + REAL(DOUBLE) :: ADD(5,5), ADD_SAVE(5,5) ! discarded-discarded block (max Nd=5) and a working copy + REAL(DOUBLE) :: RHS_MAT(5,6+6*MRBE3) ! combined RHS for the Gauss solve: [A_dr | B_d] + REAL(DOUBLE) :: A_EFF(6,6) ! Schur-reduced system matrix (only r,r entries meaningful) + REAL(DOUBLE) :: B_EFF(6,6*MRBE3) ! Schur-reduced coefficient table (only r rows meaningful) + INTEGER(LONG) :: II, JJ, KK, PIVROW ! loop/pivot indices for the Gauss elimination + REAL(DOUBLE) :: PIVVAL, FACTOR, TMPSWAP ! ********************************************************************************************************************************** @@ -301,153 +322,272 @@ SUBROUTINE RBE3_PROC ( RTYPE, REC_NO, IERR ) ! the IRBE3 grids in the "independent" set (the i points). There are up to 6 constraint eqns per RBE3 (1 each for T1, T2, T3, ! R1, R2 R3 comps in the indep set for the RBE3) +! FIX (partial-REFC coupling bug): build the FULL 6x6 coupled system A6 first, regardless of REFC. +! This is exactly the same math as before (WT6 diagonal for translation rows, EBAR for rotation +! rows, the S-terms and DX_BAR-family cross terms) -- just assembled into an explicit matrix +! instead of being written straight to RMG one REFC-selected row at a time. + DO II=1,6 + DO JJ=1,6 + A6(II,JJ) = ZERO + ENDDO + ENDDO + A6(1,1) = WT6(1); A6(1,5) = SX_DZ_BAR; A6(1,6) = -SX_DY_BAR + A6(2,2) = WT6(2); A6(2,4) = -SY_DZ_BAR; A6(2,6) = SY_DX_BAR + A6(3,3) = WT6(3); A6(3,4) = SZ_DY_BAR; A6(3,5) = -SZ_DX_BAR + A6(4,4) = EBAR_YZ; A6(4,5) = -SXY; A6(4,6) = -SZX + A6(5,5) = EBAR_ZX; A6(5,4) = -SXY; A6(5,6) = -SYZ + A6(6,6) = EBAR_XY; A6(6,4) = -SZX; A6(6,5) = -SYZ + A6(4,2) = -SY_DZ_BAR; A6(4,3) = SZ_DY_BAR + A6(5,1) = SX_DZ_BAR; A6(5,3) = -SZ_DX_BAR + A6(6,1) = -SX_DY_BAR; A6(6,2) = SY_DX_BAR + ! Mirror to guarantee exact symmetry + DO II=1,6 + DO JJ=II+1,6 + A6(JJ,II) = A6(II,JJ) + ENDDO + ENDDO +! Build the full coefficient table B6: one column per (independent grid, active component), +! for ALL 6 rows -- again regardless of REFC. TDI is computed once per independent grid here +! (previously computed inside the row loop, redundantly, once per REFC-selected row). + + NB6COLS = 0 + DO JJ=1,6*MRBE3 + B6_COL_RMG(JJ) = 0 + DO II=1,6 + B6(II,JJ) = ZERO + ENDDO + ENDDO + DO J=1,IRBE3 + CALL RDOF ( COMPS_I(J), CDOF_I ) + CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID_I(J), GRID_ID_ROW_NUM_I ) + CALL GET_GRID_NUM_COMPS ( GRID_ID_ROW_NUM_I, NUM_COMPS, SUBR_NAME ) + IF (NUM_COMPS /= 6) THEN + IERR = IERR + 1 + JERR = JERR + 1 + WRITE(ERR,1951) 'RBE3', REID, NUM_COMPS + WRITE(F06,1951) 'RBE3', REID, NUM_COMPS + FATAL_ERR = FATAL_ERR + 1 + RETURN + ENDIF + ECORD_I= GRID(GRID_ID_ROW_NUM_I,3) + IF (ECORD_I /= 0) THEN + DO K=1,NCORD + IF (ECORD_I == CORD(K,2)) THEN + ICORD_I = K + EXIT + ENDIF + ENDDO + CALL GEN_T0L ( GRID_ID_ROW_NUM_I, ICORD_I, THETAD, PHID, T0I ) + ELSE + DO K=1,3 + DO L=1,3 + T0I(K,L) = ZERO + ENDDO + T0I(K,K) = ONE + ENDDO + ENDIF + CALL MATMULT_FFF_T ( T0D, T0I, 3, 3, 3, TDI ) +! Resolve the RMG column number for each of this grid's 6 components once, up front. + CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID_I(J), IGRID ) + ROW_NUM_START_D = TDOF_ROW_START(IGRID) ! reused as "ROW_NUM_START_I" here + IF (TDOF(ROW_NUM_START_D,G_SET_COL_NUM) <= 0) THEN + WRITE(ERR,'(A,I8,A)') ' *ERROR: RBE3_PROC found no valid G-set column for grid ', AGRID_I(J), ' (independent grid)' + WRITE(F06,'(A,I8,A)') ' *ERROR: RBE3_PROC found no valid G-set column for grid ', AGRID_I(J), ' (independent grid)' + FATAL_ERR = FATAL_ERR + 1 + CALL OUTA_HERE ( 'Y' ) + ENDIF + + DO K=1,3 + NB6COLS = NB6COLS + 1 + IF (CDOF_I(K) == '1') THEN + B6_COL_RMG(NB6COLS) = (TDOF(ROW_NUM_START_D,G_SET_COL_NUM)-1) + K + ! translation column K: rows 1-3 (WRITE_L1J_123-style) + DO II=1,3 + B6(II,NB6COLS) = -WTi6(J,K)*TDI(II,K) + ENDDO + ! and rows 4-6 (WRITE_L1J_456-style translation terms) + B6(4,NB6COLS) = B6(4,NB6COLS) + WTi6(J,K)*(DZI(J)*TDI(2,K) - DYI(J)*TDI(3,K)) + B6(5,NB6COLS) = B6(5,NB6COLS) + WTi6(J,K)*(-DZI(J)*TDI(1,K) + DXI(J)*TDI(3,K)) + B6(6,NB6COLS) = B6(6,NB6COLS) + WTi6(J,K)*(DYI(J)*TDI(1,K) - DXI(J)*TDI(2,K)) + ENDIF + ENDDO + + DO K=4,6 ! rotation column K: gap-3 direct rotation coupling + NB6COLS = NB6COLS + 1 + IF (CDOF_I(K) == '1') THEN + B6_COL_RMG(NB6COLS) = (TDOF(ROW_NUM_START_D,G_SET_COL_NUM)-1) + K + KK = K - 3 + B6(4,NB6COLS) = B6(4,NB6COLS) - WTi6(J,K)*TDI(1,KK) + B6(5,NB6COLS) = B6(5,NB6COLS) - WTi6(J,K)*TDI(2,KK) + B6(6,NB6COLS) = B6(6,NB6COLS) - WTi6(J,K)*TDI(3,KK) + ENDIF + ENDDO + + ENDDO + +! Split the 6 components into "retained" (R, = REFC-selected/dependent) and "discarded" (D, not +! part of REFC) sets, then eliminate the D set from A6/B6 via a Schur complement, so the retained +! rows correctly account for their coupling to the discarded ones instead of simply ignoring it. + + NR = 0; ND = 0 + DO II=1,6 + IS_R(II) = (CDOF_D(II) == '1') + IF (IS_R(II)) THEN + NR = NR + 1 + R_IDX(NR) = II + ELSE + ND = ND + 1 + D_IDX(ND) = II + ENDIF + ENDDO + + IF (ND == 0) THEN ! Nothing to eliminate -- exact fast path, byte-identical + DO II=1,6 ! to the original (pre-fix) computation. + DO JJ=1,6 + A_EFF(II,JJ) = A6(II,JJ) + ENDDO + ENDDO + DO II=1,6 + DO JJ=1,NB6COLS + B_EFF(II,JJ) = B6(II,JJ) + ENDDO + ENDDO + ELSE + ! Build A_dd and the combined RHS [A_dr | B_d] + DO II=1,ND + DO JJ=1,ND + ADD(II,JJ) = A6(D_IDX(II),D_IDX(JJ)) + ENDDO + DO JJ=1,NR + RHS_MAT(II,JJ) = A6(D_IDX(II),R_IDX(JJ)) + ENDDO + DO JJ=1,NB6COLS + RHS_MAT(II,NR+JJ) = B6(D_IDX(II),JJ) + ENDDO + ENDDO + ! Gauss elimination with partial pivoting: solve + ! ADD * X = RHS_MAT for X (overwrite RHS_MAT with X) + DO KK=1,ND + PIVROW = KK + PIVVAL = DABS(ADD(KK,KK)) + DO II=KK+1,ND + IF (DABS(ADD(II,KK)) > PIVVAL) THEN + PIVROW = II + PIVVAL = DABS(ADD(II,KK)) + ENDIF + ENDDO + IF (PIVROW /= KK) THEN + DO JJ=1,ND + TMPSWAP = ADD(KK,JJ); ADD(KK,JJ) = ADD(PIVROW,JJ); ADD(PIVROW,JJ) = TMPSWAP + ENDDO + DO JJ=1,NR+NB6COLS + TMPSWAP = RHS_MAT(KK,JJ); RHS_MAT(KK,JJ) = RHS_MAT(PIVROW,JJ); RHS_MAT(PIVROW,JJ) = TMPSWAP + ENDDO + ENDIF + IF (DABS(ADD(KK,KK)) > EPS1) THEN + DO II=KK+1,ND + FACTOR = ADD(II,KK)/ADD(KK,KK) + DO JJ=KK,ND + ADD(II,JJ) = ADD(II,JJ) - FACTOR*ADD(KK,JJ) + ENDDO + DO JJ=1,NR+NB6COLS + RHS_MAT(II,JJ) = RHS_MAT(II,JJ) - FACTOR*RHS_MAT(KK,JJ) + ENDDO + ENDDO + ENDIF + ENDDO + ! Back-substitution + DO KK=ND,1,-1 + IF (DABS(ADD(KK,KK)) > EPS1) THEN + DO JJ=1,NR+NB6COLS + DO II=KK+1,ND + RHS_MAT(KK,JJ) = RHS_MAT(KK,JJ) - ADD(KK,II)*RHS_MAT(II,JJ) + ENDDO + RHS_MAT(KK,JJ) = RHS_MAT(KK,JJ)/ADD(KK,KK) + ENDDO + ELSE ! degenerate discarded block: no correction from this DOF + DO JJ=1,NR+NB6COLS + RHS_MAT(KK,JJ) = ZERO + ENDDO + ENDIF + ENDDO + ! A_eff = A_rr - A_rd*X ; B_eff = B_r - A_rd*Y + DO II=1,NR + DO JJ=1,NR + A_EFF(II,JJ) = A6(R_IDX(II),R_IDX(JJ)) + DO KK=1,ND + A_EFF(II,JJ) = A_EFF(II,JJ) - A6(R_IDX(II),D_IDX(KK))*RHS_MAT(KK,JJ) + ENDDO + ENDDO + DO JJ=1,NB6COLS + B_EFF(II,JJ) = B6(R_IDX(II),JJ) + DO KK=1,ND + B_EFF(II,JJ) = B_EFF(II,JJ) - A6(R_IDX(II),D_IDX(KK))*RHS_MAT(KK,NR+JJ) + ENDDO + ENDDO + ENDDO + ! Re-expand A_EFF/B_EFF back to full 1..6 row indexing + ! (rows for indices in R_IDX only; others unused/ignored) + DO II=NR,1,-1 + DO JJ=1,NR + A_EFF(R_IDX(II),R_IDX(JJ)) = A_EFF(II,JJ) + ENDDO + DO JJ=1,NB6COLS + B_EFF(R_IDX(II),JJ) = B_EFF(II,JJ) + ENDDO + ENDDO + ENDIF + +! Write terms to L1J for the constraint equations, now using the (possibly Schur-reduced) A_EFF +! and B_EFF instead of the raw pivots/S-terms/DX_BAR-family sums and per-row WRITE_L1J_123/456 +! calls. There are up to 6 constraint eqns per RBE3 (1 each for T1, T2, T3, R1, R2, R3 comps). + ITERM_RMG = 0 do_i1:DO I=1,6 cdof_dep:IF (CDOF_D(I) == '1') THEN ! The I-th component is in DDOF so write this row to RMG IROW = I -!xx CALL CALC_TDOF_ROW_NUM ( AGRID_D, ROW_NUM_START_D, 'N' ) CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID_D, IGRID ) ROW_NUM_START_D = TDOF_ROW_START(IGRID) ROW_NUM = ROW_NUM_START_D + I - 1 RMG_ROW_NUM = TDOF(ROW_NUM, M_SET_COL_NUM) IF ((RMG_ROW_NUM > 0) .AND. (RMG_COL_NUM_D(I) > 0)) THEN - ! Write coeff for the T1, T2 or T3 component at the ref pt + IF ((I == 1) .OR. (I == 2) .OR. (I == 3)) THEN - IF (DABS(WT) > EPS1) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), WT6(I) - ITERM_RMG = ITERM_RMG + 1 - ELSE - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), WT6(I) + IF (DABS(WT) <= EPS1) THEN + WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), A_EFF(I,I) ITERM_RMG = ITERM_RMG + 1 CYCLE do_i1 ENDIF ENDIF - - IF (I == 1) THEN ! Write coeffs for the R2, R3 comps at the ref pt for the 1st eqn - - IF (CDOF_D(5) /= '0') THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(5), +SX_DZ_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - - IF (CDOF_D(6) /= '0') THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(6), -SX_DY_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - - ELSE IF (I == 2) THEN ! Write coeffs for the R1, R3 comps at the ref pt for the 2nd eqn - - IF (CDOF_D(4) /= '0') THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(4), -SY_DZ_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - - IF (CDOF_D(6) /= '0') THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(6), +SY_DX_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - - ELSE IF (I == 3) THEN ! Write coeffs for the R1, R2 comps at the ref pt for the 3rd eqn - - IF (CDOF_D(4) /= '0') THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(4), +SZ_DY_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - - IF (CDOF_D(5) /= '0') THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(5), -SZ_DX_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - - ENDIF - ! Write coeffs for the R1, R2 and R3 comps at the ref pt for eqns 4,5,6 - IF (I == 4) THEN - IF (DABS(EBAR_YZ) > EPS1) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), EBAR_YZ - ITERM_RMG = ITERM_RMG + 1 - ELSE - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), ONE - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(5) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(5), -SXY - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(6) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(6), -SZX - ITERM_RMG = ITERM_RMG + 1 - ENDIF - ENDIF - - IF (I == 5) THEN - IF (DABS(EBAR_ZX) > EPS1) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), EBAR_ZX - ITERM_RMG = ITERM_RMG + 1 + ! Pivot term (own column), with the same near-zero fallback + ! to a unit pivot the original code used for rows 4-6 + IF ((I == 4) .OR. (I == 5) .OR. (I == 6)) THEN + IF (DABS(A_EFF(I,I)) > EPS1) THEN + WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), A_EFF(I,I) ELSE WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), ONE - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(4) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(4), -SXY - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(6) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(6), -SYZ - ITERM_RMG = ITERM_RMG + 1 - ENDIF - ENDIF - - IF (I == 6) THEN - IF (DABS(EBAR_XY) > EPS1) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), EBAR_XY - ITERM_RMG = ITERM_RMG + 1 - ELSE - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), ONE - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(4) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(4), -SZX - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(5) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(5), -SYZ - ITERM_RMG = ITERM_RMG + 1 - ENDIF - ENDIF - - IF (I == 4) THEN - IF (RMG_COL_NUM_D(2) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(2), -SY_DZ_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(3) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(3), +SZ_DY_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - ENDIF - - IF (I == 5) THEN - IF (RMG_COL_NUM_D(1) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(1), +SX_DZ_BAR - ITERM_RMG = ITERM_RMG + 1 - ENDIF - IF (RMG_COL_NUM_D(3) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(3), -SZ_DX_BAR - ITERM_RMG = ITERM_RMG + 1 ENDIF + ELSE + WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(I), A_EFF(I,I) ENDIF - - IF (I == 6) THEN - IF (RMG_COL_NUM_D(1) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(1), -SX_DY_BAR - ITERM_RMG = ITERM_RMG + 1 + ITERM_RMG = ITERM_RMG + 1 + ! Cross terms to the OTHER retained (REFC-selected) comps + DO JJ=1,6 + IF ((JJ /= I) .AND. (CDOF_D(JJ) == '1')) THEN + IF (A_EFF(I,JJ) /= ZERO) THEN + WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(JJ), A_EFF(I,JJ) + ITERM_RMG = ITERM_RMG + 1 + ENDIF ENDIF - IF (RMG_COL_NUM_D(2) > 0) THEN - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_D(2), +SY_DX_BAR + ENDDO + ! Independent-grid terms + DO JJ=1,NB6COLS + IF ((B6_COL_RMG(JJ) > 0) .AND. (B_EFF(I,JJ) /= ZERO)) THEN + WRITE(L1J) RMG_ROW_NUM, B6_COL_RMG(JJ), B_EFF(I,JJ) ITERM_RMG = ITERM_RMG + 1 ENDIF - ENDIF + ENDDO ELSE IF (RMG_ROW_NUM == 0) THEN @@ -466,50 +606,6 @@ SUBROUTINE RBE3_PROC ( RTYPE, REC_NO, IERR ) ENDIF ENDIF -do_j1: DO J=1,IRBE3 ! Cycle over "indep" terms (some here may actually be dep elsewhere) - ! Get T0I (transforms global vector at AGRID_I to basic) - CALL RDOF ( COMPS_I(J), CDOF_I ) - - CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID_I(J), GRID_ID_ROW_NUM_I ) - - CALL GET_GRID_NUM_COMPS ( GRID_ID_ROW_NUM_I, NUM_COMPS, SUBR_NAME ) - IF (NUM_COMPS /= 6) THEN - IERR = IERR + 1 - JERR = JERR + 1 - WRITE(ERR,1951) 'RBE3', REID, NUM_COMPS - WRITE(F06,1951) 'RBE3', REID, NUM_COMPS - FATAL_ERR = FATAL_ERR + 1 - RETURN - ENDIF - - ECORD_I= GRID(GRID_ID_ROW_NUM_I,3) - IF (ECORD_I /= 0) THEN - DO K=1,NCORD - IF (ECORD_I == CORD(K,2)) THEN - ICORD_I = K - EXIT - ENDIF - ENDDO - CALL GEN_T0L ( GRID_ID_ROW_NUM_I, ICORD_I, THETAD, PHID, T0I ) - ELSE - DO K=1,3 - DO L=1,3 - T0I(K,L) = ZERO - ENDDO - T0I(K,K) = ONE - ENDDO - ENDIF - - CALL MATMULT_FFF_T ( T0D, T0I, 3, 3, 3, TDI ) - - IF ((I == 1) .OR. (I == 2) .OR. (I == 3)) THEN - CALL WRITE_L1J_123 ( I, J, ITERM_RMG, G_SET_COL_NUM, RMG_ROW_NUM, WTi6, AGRID_I, CDOF_I, COMPS_I, TDI ) - ELSE - CALL WRITE_L1J_456 ( I, J, ITERM_RMG, G_SET_COL_NUM, RMG_ROW_NUM, WTi6, AGRID_I, CDOF_I, DXI, DYI, DZI, TDI ) - ENDIF - - ENDDO do_j1 - ENDIF cdof_dep ENDDO do_i1 @@ -527,10 +623,6 @@ SUBROUTINE RBE3_PROC ( RTYPE, REC_NO, IERR ) RETURN ! ********************************************************************************************************************************** - 1503 FORMAT(' *ERROR 1503: PROGRAMMING ERROR IN SUBROUTINE ',A & - ,/,14X,' FOR RBE3 ',I8,', THE COUNT OF RECORDS WRITTEN TO FILE' & - ,/,15X,A & - ,/,14X,' IS ITERM_RMG = ',I8,' BUT SHOULD BE NTERM_RMG = ',I8,' DETERMINED IN ANOTHER SUBR') 1509 FORMAT(' *ERROR 1509: PROGRAMMING ERROR IN SUBROUTINE ',A & ,/,15X,A8,' RIGID ELEMENT NUMBER ',I8,', DEPENDENT GRID NUMBER ',I8,', COMPONENT ',I2 & @@ -543,196 +635,7 @@ SUBROUTINE RBE3_PROC ( RTYPE, REC_NO, IERR ) 1951 FORMAT(' *ERROR 1951: ',A,I8,' USES GRID ',I8,' WHICH IS A SCALAR POINT. SCALAR POINTS NOT ALLOWED FOR THIS ELEM TYPE') - - - - - -! ********************************************************************************************************************************** - - CONTAINS - -! ################################################################################################################################## - - SUBROUTINE WRITE_L1J_123 ( I, J, ITERM_RMG, G_SET_COL_NUM, RMG_ROW_NUM, WTi6, AGRID_I, CDOF_I, COMPS_I, TDI ) - - USE PENTIUM_II_KIND, ONLY : LONG, DOUBLE - USE IOUNT1, ONLY : L1J - USE SCONTR, ONLY : FATAL_ERR, MRBE3, NGRID - USE DOF_TABLES, ONLY : TDOF, TDOF_ROW_START - USE MODEL_STUF, ONLY : GRID_ID - - CHARACTER( 1*BYTE), INTENT(IN) :: CDOF_I(6) ! An output from subr RDOF (= 1 if a displ comp 1-6 is in COMPS_I) - - INTEGER(LONG), INTENT(IN) :: I,J ! DO loop indices - INTEGER(LONG), INTENT(IN) :: AGRID_I(MRBE3) ! Indep grid ID (actual) read from a record of file LINK1F - INTEGER(LONG), INTENT(IN) :: COMPS_I(MRBE3) ! Dipsl components associated with indep grid, AGRID_I - INTEGER(LONG), INTENT(IN) :: G_SET_COL_NUM ! Col no., in TDOF array, of the G-set DOF list - INTEGER(LONG), INTENT(IN) :: RMG_ROW_NUM ! Row no. of a term in array RMG - INTEGER(LONG), INTENT(INOUT) :: ITERM_RMG ! Count of number of records written to L1J (should be NTERM_RMG at end) - INTEGER(LONG) :: IGRID ! Internal grid ID - INTEGER(LONG) :: K ! DO loop index - INTEGER(LONG) :: RMG_COL_NUM_I ! Col no. of a term in array RMG - INTEGER(LONG) :: ROW_NUM ! A row number in array TDOF - INTEGER(LONG) :: ROW_NUM_START_I ! DOF number where TDOF data begins for a grid - - REAL(DOUBLE) , INTENT(IN) :: TDI(3,3) ! TOD'*T0I - REAL(DOUBLE) , INTENT(IN) :: WTi6(MRBE3,6) ! Weight value for an indep grid (PER-DOF) - -! ********************************************************************************************************************************** - - DO K=1,3 - - IF (CDOF_I(K) == '1') THEN -!xx CALL CALC_TDOF_ROW_NUM ( AGRID_I(J), ROW_NUM_START_I, 'N' ) - CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID_I(J), IGRID ) - ROW_NUM_START_I = TDOF_ROW_START(IGRID) - ROW_NUM = ROW_NUM_START_I + K - 1 - RMG_COL_NUM_I = TDOF(ROW_NUM,G_SET_COL_NUM) - - IF (RMG_COL_NUM_I > 0) THEN ! NB *** new 10/03/21 Change WT (below) to WT6(K) - WRITE(L1J) RMG_ROW_NUM, RMG_COL_NUM_I, -WTi6(J,K)*TDI(I,K) - ITERM_RMG = ITERM_RMG + 1 - ELSE - WRITE(ERR,1513) 'RBE3_PROC', AGRID_I(J) ,COMPS_I(J), RMG_COL_NUM_I - WRITE(F06,1513) 'RBE3_PROC', AGRID_I(J) ,COMPS_I(J), RMG_COL_NUM_I - FATAL_ERR = FATAL_ERR + 1 - CALL OUTA_HERE ( 'Y' ) - ENDIF - - ENDIF - - ENDDO - -! ********************************************************************************************************************************** - 1513 FORMAT(' *ERROR 1513: PROGRAMMING ERROR IN SUBROUTINE ',A & - ,/,14X,' COL NUMBER IN ARRAY RMG CALCULATED FOR GRID ',I8,', COMPONENT ',I2,' MUST BE > 0 BUT IS = ',I8) - - -! ********************************************************************************************************************************** - - END SUBROUTINE WRITE_L1J_123 - -! ################################################################################################################################## - - SUBROUTINE WRITE_L1J_456 ( I, J, ITERM_RMG, G_SET_COL_NUM, RMG_ROW_NUM, WTi6, AGRID_I, CDOF_I, DXI, DYI, DZI, TDI ) - - USE PENTIUM_II_KIND, ONLY : LONG, DOUBLE - USE IOUNT1, ONLY : L1J - USE SCONTR, ONLY : FATAL_ERR, NGRID, MRBE3 - USE CONSTANTS_1, ONLY : ONE, ZERO - USE DOF_TABLES, ONLY : TDOF, TDOF_ROW_START - USE MODEL_STUF, ONLY : GRID_ID - - CHARACTER( 1*BYTE), INTENT(IN) :: CDOF_I(6) ! An output from subr RDOF (= 1 if a displ comp 1-6 is in COMPS_I) - - INTEGER(LONG), INTENT(IN) :: I,J ! DO loop indices - INTEGER(LONG), INTENT(IN) :: AGRID_I(MRBE3) ! Indep grid ID (actual) read from a record of file LINK1F - INTEGER(LONG), INTENT(IN) :: G_SET_COL_NUM ! Col no., in TDOF array, of the G-set DOF list - INTEGER(LONG), INTENT(IN) :: RMG_ROW_NUM ! Row no. of a term in array RMG - INTEGER(LONG), INTENT(INOUT) :: ITERM_RMG ! Count of number of records written to L1J (should be NTERM_RMG at end) - INTEGER(LONG) :: IGRID ! Internal grid ID - INTEGER(LONG) :: K ! DO loop index - INTEGER(LONG) :: RMG_COL_NUM_START ! Col no. of a term in array RMG - INTEGER(LONG) :: ROW_NUM_START_I ! DOF number where TDOF data begins for a grid - - REAL(DOUBLE) , INTENT(IN) :: DXI(MRBE3) ! Distances from ref pt to pt i in X global directions at ref pt - REAL(DOUBLE) , INTENT(IN) :: DYI(MRBE3) ! Distances from ref pt to pt i in Y global directions at ref pt - REAL(DOUBLE) , INTENT(IN) :: DZI(MRBE3) ! Distances from ref pt to pt i in Z global directions at ref pt - REAL(DOUBLE) , INTENT(IN) :: WTi6(MRBE3,6) ! Weight value for an indep grid (PER-DOF) - REAL(DOUBLE) , INTENT(IN) :: TDI(3,3) ! T0D'*T0I -- combines the dependent grid's and - ! this independent grid's own displacement coordinate - ! systems. Needed whenever either grid has a non-basic CD. - REAL(DOUBLE) :: COEF ! generalized coefficient - -! ********************************************************************************************************************************** -!xx CALL CALC_TDOF_ROW_NUM ( AGRID_I(J), ROW_NUM_START_I, 'N' ) - CALL GET_ARRAY_ROW_NUM ( 'GRID_ID', SUBR_NAME, NGRID, GRID_ID, AGRID_I(J), IGRID ) - ROW_NUM_START_I = TDOF_ROW_START(IGRID) - RMG_COL_NUM_START = TDOF(ROW_NUM_START_I,G_SET_COL_NUM) - - IF (RMG_COL_NUM_START <= 0) THEN - WRITE(ERR,1513) 'RBE3_PROC', AGRID_I(J) ,'1', RMG_COL_NUM_START - WRITE(F06,1513) 'RBE3_PROC', AGRID_I(J) ,'1', RMG_COL_NUM_START - FATAL_ERR = FATAL_ERR + 1 - CALL OUTA_HERE ( 'Y' ) - ENDIF - - IF (I == 4) THEN ! Rotation about x, i.e. in yz (23) plane - - DO K=1,3 - IF (CDOF_I(K) == '1') THEN - COEF = WTi6(J,K)*(DZI(J)*TDI(2,K) - DYI(J)*TDI(3,K)) - IF (COEF /= ZERO) THEN - WRITE(L1J) RMG_ROW_NUM, (RMG_COL_NUM_START-1)+K, COEF - ITERM_RMG = ITERM_RMG +1 - ENDIF - ENDIF - ENDDO - - DO K=1,3 ! Independent grid's own rotation averages in through TDI too - IF (CDOF_I(K+3) == '1') THEN - COEF = -WTi6(J,K+3)*TDI(1,K) - IF (COEF /= ZERO) THEN - WRITE(L1J) RMG_ROW_NUM, (RMG_COL_NUM_START-1)+(K+3), COEF - ITERM_RMG = ITERM_RMG +1 - ENDIF - ENDIF - ENDDO - - ELSE IF (I == 5) THEN ! Rotation about y, i.e. in zx (31) plane - - DO K=1,3 - IF (CDOF_I(K) == '1') THEN - COEF = WTi6(J,K)*(-DZI(J)*TDI(1,K) + DXI(J)*TDI(3,K)) - IF (COEF /= ZERO) THEN - WRITE(L1J) RMG_ROW_NUM, (RMG_COL_NUM_START-1)+K, COEF - ITERM_RMG = ITERM_RMG +1 - ENDIF - ENDIF - ENDDO - - DO K=1,3 - IF (CDOF_I(K+3) == '1') THEN - COEF = -WTi6(J,K+3)*TDI(2,K) - IF (COEF /= ZERO) THEN - WRITE(L1J) RMG_ROW_NUM, (RMG_COL_NUM_START-1)+(K+3), COEF - ITERM_RMG = ITERM_RMG +1 - ENDIF - ENDIF - ENDDO - - ELSE IF (I == 6) THEN ! Rotation about z, i.e. in xy (12) plane - - DO K=1,3 - IF (CDOF_I(K) == '1') THEN - COEF = WTi6(J,K)*(DYI(J)*TDI(1,K) - DXI(J)*TDI(2,K)) - IF (COEF /= ZERO) THEN - WRITE(L1J) RMG_ROW_NUM, (RMG_COL_NUM_START-1)+K, COEF - ITERM_RMG = ITERM_RMG +1 - ENDIF - ENDIF - ENDDO - - DO K=1,3 - IF (CDOF_I(K+3) == '1') THEN - COEF = -WTi6(J,K+3)*TDI(3,K) - IF (COEF /= ZERO) THEN - WRITE(L1J) RMG_ROW_NUM, (RMG_COL_NUM_START-1)+(K+3), COEF - ITERM_RMG = ITERM_RMG +1 - ENDIF - ENDIF - ENDDO - - ENDIF - -! ********************************************************************************************************************************** - 1513 FORMAT(' *ERROR 1513: PROGRAMMING ERROR IN SUBROUTINE ',A & - ,/,14X,' COL NUMBER IN ARRAY RMG CALCULATED FOR GRID ',I8,', COMPONENT ',A,' MUST BE > 0 BUT IS = ',I8) - - ! ********************************************************************************************************************************** - END SUBROUTINE WRITE_L1J_456 END SUBROUTINE RBE3_PROC