SUBROUTINE ABQMAIN C==================================================================== C This program must be compiled and linked with the command: C abaqus make job=fpert C Run the program using the command: C abaqus fpert C==================================================================== C C PURPOSE: C This program computes a perturbed mesh based on a user-specified C perturbation factor. The original coordinate data and C eigenvectors are read from an ABAQUS results (.fil) file. C PROMPTS: C 1. `Enter the name of the results file (w/o .fil):' C 2. `Enter the mode shape(s) to be used in calculating the C perturbed mesh (zero when finished):' C 3. `Enter the imperfection factor to be introduced into the C geometry for this eigenmode:' C C==================================================================== C C INPUT FILE -- `FNAME'.fil C OUTPUT FILE -- fpert002.015 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 = Described in Section 7.0.0 of the Verification manual C JRRAY = Described in Section 7.0.0 of the Verification manual C LRUNIT = Described in Section 7.0.0 of the Verification manual 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==================================================================== 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 C=================================================================== PARAMETER (ITOTAL = 8000) C C==================================================================== C DIMENSION DISP(6,ITOTAL), COORD(3,ITOTAL) DIMENSION INODE(ITOTAL), IDOF(30), JEIGNO(10) CHARACTER FNAME*80,OUTFILE*(*) PARAMETER (OUTFILE = 'fpert002.015') 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 output file. C 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(5,'(A)', IOSTAT = J ) FNAME IF (J .NE. 0) GOTO 950 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=================================================================== WRITE(*,2010) NODEMAX, IELMAX WRITE(*,2015) IEIGNO 5 READ(5,*,ERR = 950) JEIGNO(I) IF (JEIGNO(I) .EQ. 0) GOTO 10 I=I+1 GOTO 5 C 10 CONTINUE 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==================================================================== 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 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. C C==================================================================== C ELSE IF (JRRAY(1,2) .EQ. 1980) THEN IEIGNO = JRRAY(1,3) DO J = 1, I-1 IF (JEIGNO(J) .EQ. IEIGNO) K = J 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) = ARRAY(4) DISP(2,I101) = ARRAY(5) IF (IDOF(3) .NE. 0) DISP(3,I101) = ARRAY(IDOF(3)+3) IF (IDOF(4) .NE. 0) DISP(4,I101) = ARRAY(IDOF(4)+3) IF (IDOF(5) .NE. 0) DISP(5,I101) = ARRAY(IDOF(5)+3) IF (IDOF(6) .NE. 0) DISP(6,I101) = ARRAY(IDOF(6)+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(5,*,IOSTAT = J) FACTOR IF (J .NE. 0) GOTO 950 ICYCLE = ICYCLE + 1 I101 = 0 C 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 Subroutine NODEGEN is the actual mesh generator. C==================================================================== CALL NODEGEN(COORD,DISP,I1901,FACTOR) 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=================================================================== C Output the coordinates of the perturbed mesh and close the file. C=================================================================== WRITE(*,100) OUTFILE 100 FORMAT(//,2X,'The perturbed mesh data are being written to:', 1 1X,A,//) C DO K = 1, NODEMAX WRITE(15,110) INODE(K), (COORD(J,K),J = 1, ITRANS) 110 FORMAT(I6,3(',',1PE14.6)) ENDDO CLOSE (15) WRITE(*,120) 120 FORMAT(//,2X,' . . . PROGRAM FINISHED SUCCESSFULLY . . . ') RETURN 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 P E R T |', 4 /,' | P E R T U R B E D M E S H |', 5 /,' | G E N E R A T O R |', 6 /,' | |', 7 /,' +----------------------------------------------+',//) 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) RETURN 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 = Original 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(6,*) 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