SUBROUTINE ABQMAIN C C PROGRAM FMAXIMP C==================================================================== C C INPUT FILE -- `FNAME'.fil, max_round_input.dat C C OUTPUT FILE -- `IMPFNAME` C C==================================================================== C C The use of ABA_PARAM.INC eliminates the need to have different C versions of the code for single and double precision. C ABA_PARAM.INC defines an appropriate IMPLICIT REAL statement and C sets the value of NPRECD to 1 or 2, depending on whether the C machine uses single or double precision. C C==================================================================== C ARRAY = Real array containing values read from results file C (.fil). Equivalenced to JRRAY. C JRRAY = Integer array containing values read from results file C (.fil). Equivalenced to ARRAY. C LRUNIT = Array containing unit number and format of results file C LRUNIT(1,*) --> Unit number of input file. C LRUNIT(2,*) --> Format of input file. C DISP = Contains the eigenvector for a particular eigenmode C COORD = Original coordinate data C INODE = Original node label C IDOF = DOF for the element C JEIGNO = Array of mode shapes used for calculating the perturbed C mesh C FNAME = Name of the results file C NODEMAX = Number of nodes in the model C IELMAX = Number of elements in the model C==================================================================== C INCLUDE 'ABA_PARAM.INC' DIMENSION ARRAY(513), JRRAY(NPRECD,513), LRUNIT(2,1) EQUIVALENCE (ARRAY(1), JRRAY(1,1)) C=================================================================== C ITOTAL must be greater than or equal to the number of nodes in the C model; MEIGEN must be greater than or equal to the number of eigenmodes. C=================================================================== PARAMETER (ITOTAL = 8000,MEIGEN=20) C C==================================================================== C DIMENSION DISP(3,ITOTAL,MEIGEN), COORD(3,ITOTAL),FCTR(MEIGEN) DIMENSION ALPHA(MEIGEN), NPAIR(MEIGEN) DIMENSION OLD_COORD(3,ITOTAL) DIMENSION INODE(ITOTAL), IDOF(30), JEIGNO(MEIGEN),EIGVAL(MEIGEN) CHARACTER FNAME*80,OUTFILE*80,INPFILE*80 PARAMETER (INPFILE = 'max_round_input.dat') C C==================================================================== C Define flags and counters. C C==================================================================== ICYCLE = 0 I1901 = 0 I101 = 0 I = 1 K = 1 C C==================================================================== C Define file access variables. C C==================================================================== NRU = 1 LRUNIT(1,NRU) = 8 LRUNIT(2,NRU) = 2 LOUTF = 0 C C==================================================================== C Open input and output files. C C==================================================================== OPEN(UNIT=16,FILE=INPFILE,STATUS='OLD') C READ(16,'(A)', IOSTAT = J ) OUTFILE IF (J .NE. 0) GOTO 950 C OPEN(UNIT=15,FILE=OUTFILE,STATUS='UNKNOWN',IOSTAT = J) IF (J .NE. 0) THEN WRITE(*,900) OUTFILE GOTO 950 ENDIF C C==================================================================== C Get the name of the results (.fil) file. C C==================================================================== WRITE(*,2000) WRITE(6,*) ' Enter the name of the results file (w/o .fil):' READ(16,'(A)', IOSTAT = J ) FNAME IF (J .NE. 0) GOTO 950 READ(16,*) N_CIRCUM READ(16,*) N_LENGTH C C==================================================================== C Access ABAQUS libraries to set up input file. C C==================================================================== CALL INITPF (FNAME, NRU, LRUNIT, LOUTF) C JUNIT = LRUNIT(1,NRU) C CALL DBRNU (JUNIT) C C==================================================================== C Read a record from the input file. C On the first pass through the file obtain the number of nodes for C a diagnostic check. C==================================================================== CALL DBFILE (0, ARRAY, JRCD) DO WHILE (JRCD .EQ. 0) IF (JRRAY(1,2) .EQ. 1980) IEIGNO = JRRAY(1,3) IF (JRRAY(1,2) .EQ. 1921 ) THEN NODEMAX = JRRAY(1,8) IELMAX = JRRAY(1,7) ICYCLE = ICYCLE +1 ENDIF CALL DBFILE (0, ARRAY, JRCD) ENDDO C CALL DBFILE (2, ARRAY, JRCD) C=================================================================== C User is given a choice of eigenmodes for calculating the perturbed C mesh. C C=================================================================== READ(16,*)ACT_IMP C WRITE(*,2010) NODEMAX, IELMAX C 5 READ(16,*,ERR = 950) JEIGNO(I) IF (JEIGNO(I) .EQ. 0) GOTO 10 I=I+1 GOTO 5 C 10 CONTINUE NMODES = I-1 IF (NMODES.GT.MEIGEN)THEN WRITE (6,*) 'INCREASE MEIGEN TO' ,NMODES STOP END IF C CALL DBFILE (0, ARRAY, JRCD) C DO WHILE (JRCD .EQ. 0) C C==================================================================== C If this is the first pass through the file and the current record C is the nodal coordinate record, then read the original nodal C coordinates and the node numbers. Make sure that the third C coordinate exists before saving it. C C==================================================================== IF (JRRAY(1,2) .EQ. 1901 .AND. ICYCLE .LE. 1) THEN I1901 = I1901 + 1 INODE(I1901) = JRRAY(1,3) COORD(1,I1901) = ARRAY(4) COORD(2,I1901) = ARRAY(5) COORD(3,I1901) = 0.0D0 IF (JRRAY(1,1) .GE. 6) COORD(3,I1901) = ARRAY(6) C DO II = 1, 3 OLD_COORD(II,I1901) = COORD(II,I1901) END DO C C==================================================================== C If this is the first pass through the file and the current record C is the active degree of freedom record, save the active d.o.f. C If the d.o.f. is active in the model, IDOF(XX) equals the C position of d.o.f. XX in the output arrays. If the d.o.f. is not C active, IDOF(XX) is zero for d.o.f. XX (i.e., for planar models C IDOF(1) = 1, IDOF(2) = 2, IDOF(3) = 0, IDOF(4) = 0, IDOF(5) = 0, C IDOF(6) = 3, etc.). ITRANS equals the number of translational C d.o.f.'s in the model. C C==================================================================== ELSE IF (JRRAY(1,2) .EQ. 1902) THEN DO 15 IXX = 1, JRRAY(1,1)-2 IDOF(IXX) = JRRAY(1,IXX+2) 15 CONTINUE ITRANS = 3 IF (IDOF(3) .EQ. 0) ITRANS = 2 C C==================================================================== C If the current record is the modal record, save the current C eigenvalue number and eigenvalue. C C==================================================================== C ELSE IF (JRRAY(1,2) .EQ. 1980) THEN IEIGNO = JRRAY(1,3) DO J = 1, NMODES IF (JEIGNO(J) .EQ. IEIGNO) THEN K = J EIGVAL(K) = ARRAY(4) END IF ENDDO C C==================================================================== C If the current record is the displacement record and the current C eigenvalue was requested, read the displacement data. The data C will be in the coordinate system specified in the C `*NODE FILE,GLOBAL=' option. If nodal transformations were C performed and GLOBAL=NO was used, the displacements will be in C the local system. If nodal transformations were used and C GLOBAL=YES, the results will be in the global system. In all C other cases the results will be in the global system. Also, C make sure that degrees of freedom are active in the model before C saving them in the appropriate array location. C C==================================================================== C ELSE IF (JRRAY(1,2) .EQ. 101 .AND. IEIGNO .EQ. JEIGNO(K)) THEN I101 = I101 + 1 DISP(1,I101,K) = ARRAY(4) DISP(2,I101,K) = ARRAY(5) IF (IDOF(3) .NE. 0) DISP(3,I101,K) = ARRAY(IDOF(3)+3) IF (INODE(I101) .EQ. 0) INODE(I101) = JRRAY(1,3) IF( I101 .EQ. NODEMAX ) THEN WRITE(6,16) JEIGNO(K) 16 FORMAT(/,2X,'Nodal coordinate data being computed for', 1 ' eigenvalue . . .',I5) C C==================================================================== C FACTOR should be entered as a perturbation factor in terms of a C percentage value multiplied by a geometric parameter C (e.g., shell thickness). C==================================================================== C WRITE(6,2020) READ(16,*,IOSTAT = J) FCTR(K) FACTOR = FCTR(K) IF (J .NE. 0) GOTO 950 ICYCLE = ICYCLE + 1 I101 = 0 C C==================================================================== ENDIF ENDIF C CALL DBFILE(0, ARRAY, JRCD) ENDDO C C==================================================================== C Next line is added for diagnostics. C C==================================================================== IF (NODEMAX .NE. I1901) THEN WRITE(*,2025) NODEMAX,I1901 GOTO 950 ENDIF C IF (ICYCLE .LE. 1) THEN WRITE(*,30) 30 FORMAT(2X,'. . . NO EIGENVECTORS WERE FOUND . . .',/ 1 2X,'The input file for the buckling analysis must contain:', 2 /,2X,'--> *NODE FILE <--', 3 /,2X,'--> U <--') GOTO 950 ENDIF C==================================================================== C Determine if there are repeated eigenvalues. C C==================================================================== TOL = 1D-6 NP = 0 DO I = 1, NMODES NPAIR(I) = 0 ENDDO C DO I=1, NMODES-1 IF (NPAIR(I) .EQ. 0) THEN EIG1 = EIGVAL(I) DO J = I+1, NMODES EIG2 = EIGVAL(J) DIFF = ABS ((EIG1 - EIG2)/EIG1) IF (DIFF.GT.TOL) GO TO 100 NP = NP+1 NPAIR(I) = NP NPAIR(J) = -NP GO TO 101 100 END DO ENDIF 101 END DO C==================================================================== C Determine the scale factors for the eigenvectors. C C Each of the eigenvectors corresponding to the repeated C eigenvalues will be scaled in such a way that the maximum radial C imperfection is introduced along the X-axis when both the eigenvectors C are used to introduce the imperfection. C C The maximum radial imperfection for the unique eigenvectors and the C linear combination of the repeated eigenvectors is scaled to one. C==================================================================== C DO I = 1, NMODES C IF (NPAIR(I).GE.0) THEN C RAD_MAX = 0.D0 ALPHA(I)= 1.D0 TAN_MAX = 0.D0 C C=================================================================== C Find the node with the maximum radial displacement. C=================================================================== C DO K = 1, NODEMAX C=================================================================== C Determine the local radial and tangential directions of node K. C=================================================================== C RADIUS = SQRT((OLD_COORD(1,K)**2)+(OLD_COORD(2,K)**2)) C UNIT_R1 = OLD_COORD(1,K)/RADIUS UNIT_R2 = OLD_COORD(2,K)/RADIUS C UNIT_T1 = -UNIT_R2 UNIT_T2 = UNIT_R1 C C================================================================== C Project the displacement vector onto the local directions to C obtain the radial and tangential displacements. C================================================================== C RAD_DISP=DISP(1,K,I)*(UNIT_R1)+ * DISP(2,K,I)*(UNIT_R2) C TAN_DISP=DISP(1,K,I)*(UNIT_T1)+ * DISP(2,K,I)*(UNIT_T2) C IF (ABS(RAD_DISP).GT.ABS(RAD_MAX))THEN RAD_MAX=RAD_DISP TAN_MAX=TAN_DISP N_RMAX= K END IF C END DO C C================================================================== C Determine the scale factors. C C For unique eigenvalues, scale so that max UR=1.0. C================================================================== C IF (NPAIR(I).EQ.0) THEN IF (ABS(RAD_MAX).GT.0) THEN ALPHA(I)=ALPHA(I)/RAD_MAX END IF C================================================================== C For repeated eigenvalues... C================================================================== ELSE C C================================================================== C Calculate the corresponding node on the x-axis (NodeX) and save U. C================================================================== C MAXROW = (N_RMAX/N_CIRCUM) C IF (MOD(N_RMAX,N_CIRCUM) .EQ. 0) THEN NODEX = (N_CIRCUM*(MAXROW-1))+1 ELSE NODEX = N_CIRCUM*MAXROW+1 ENDIF C X1 = DISP(1,NODEX,I) Y1 = DISP(2,NODEX,I) C C================================================================== C Save U for this node in the corresponding eigenvector C for the repeated eigenvalue. C================================================================== C DO J=I+1, NMODES IF (NPAIR(J) .EQ. -1*NPAIR(I))THEN X2 = DISP(1, NODEX, J) Y2 = DISP(2, NODEX, J) ALPHA(J)=1.D0 GO TO 150 ENDIF END DO 150 CONTINUE C================================================================== C Check the determinant (>0) and calculate the scale factors. C================================================================== DETERM = (X1*Y2-Y1*X2) IF (ABS(DETERM).GT.0.D0) THEN ALPHA(I) = (Y2*RAD_MAX-X2*TAN_MAX)/DETERM ALPHA(J) = (X1*TAN_MAX-Y1*RAD_MAX)/DETERM ELSE WRITE(6,*) 'DETERMINANT IS ZERO; WILL NOT SCALE' END IF IF (ABS(RAD_MAX).GT.0.D0) THEN ALPHA(I)=ALPHA(I)/RAD_MAX ALPHA(J)=ALPHA(J)/RAD_MAX END IF END IF ENDIF END DO C==================================================================== C Compute the perturbed mesh. This section assumes that nodal C displacements are in the GLOBAL coordinate system. If they are C not, the correct transformations should be applied prior to C perturbing the mesh. The user should supply this coding. Also, C only the translational degrees of freedom should be used for the C perturbation (X = Xo + u). C==================================================================== C DO I = 1, NMODES FACTOR = FCTR(I)*ALPHA(I) CALL NODEGEN(COORD,DISP(1,1,I),I1901,FACTOR) END DO C C==================================================================== C Compute the max out-of-round imperfection and close the file. C==================================================================== C DMAX_IMP=0.D0 C DO K = 1, NODEMAX C R_OLD = SQRT(OLD_COORD(1,K)**2 + OLD_COORD(2,K)**2) R_NEW = SQRT( COORD(1,K)**2 + COORD(2,K)**2) TEMP=ABS(R_NEW - R_OLD) IF (TEMP.GT.DMAX_IMP) THEN DMAX_IMP = TEMP NN=INODE(K) END IF ENDDO C C=================================================================== C Final result. C=================================================================== C WRITE(*,*) 'MAX IMPERFECTION = ',DMAX_IMP,' AT NODE ',NN DO K = 1, NMODES FACTOR = FCTR(K)*ALPHA(K) WRITE(15,120)JEIGNO(K), FACTOR*(ACT_IMP/DMAX_IMP) END DO 120 FORMAT(I6,',',E16.9) C CLOSE (15) CLOSE (16) STOP ' . . . PROGRAM FINISHED SUCCESSFULLY . . . ' C C==================================================================== C 900 FORMAT(//, 1 /,2X,'TROUBLE OPENING FILE',1X,A) 950 WRITE(*,1000) 1000 FORMAT(//, 1 /,2X,' . . . TROUBLE READING DATA . . . ', 2 /,2X,' . . . PROGRAM STOPPED . . . ',/) 2000 FORMAT(//, 1 /,' +----------------------------------------------+', 2 /,' | |', 3 /,' | P R O G R A M --- F M A X I M P |', 4 /,' | |', 5 /,' +----------------------------------------------+',//) 2010 FORMAT(//, 1 /,2X,'Number of nodes in model . . . . . . . ',I5, 2 /,2X,'Number of elements in model . . . . . . . ',I5) 2015 FORMAT(//, 1 /,2X,'Number of mode shapes available . . . . . . .',I5, 2 //,2X,'Enter the mode shape(s) to be used in calculating', 3 /,2x,'the perturbed mesh (zero when finished):') 2020 FORMAT(//, 1 /,2X,'Enter the imperfection factor to be introduced ', 2 /,2X,'into the geometry for this eigenmode: ') 2025 FORMAT(//, 1 /,2X,'. . . TROUBLE READING COORDINATE DATA . . . ', 2 /,2X,'Number of coordinates in model . . . . . .',I5, 3 /,2X,'Number of coordinates read . . . . . . . .',I5) STOP END C==================================================================== C SUBROUTINE NODEGEN(COORD,DISP,I1901,FACTOR) C C==================================================================== C PURPOSE: Defines new coordinate data based upon a fraction of the C eigenvector obtained in a buckling analysis. C C INPUT: C C COORD = Current coordinate data C DISP = Displacement data (eigenvector) C I1901 = Total number of nodes C FACTOR = Imperfection factor (e.g., percentage of shell C thickness) C==================================================================== C INCLUDE 'ABA_PARAM.INC' DIMENSION COORD(3,*),DISP(3,*) C DO I = 1, I1901 COORD(1,I) = COORD(1,I) + FACTOR*DISP(1,I) COORD(2,I) = COORD(2,I) + FACTOR*DISP(2,I) COORD(3,I) = COORD(3,I) + FACTOR*DISP(3,I) ENDDO RETURN END