      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
