SUBROUTINE ABQMAIN
C====================================================================
C This program must be compiled and linked with the command:
C abaqus make job=felbow
C Run the program using the command:
C abaqus felbow
C==================================================================
C
C PURPOSE:
C
C This program reads results stored in an ABAQUS results (.fil)
C file and creates an input file for use in the Visualization
C module of ABAQUS/CAE (ABAQUS/Viewer) that permits
C postprocessing and visualization of elbow element results.
C
C
C Input file names: `fname.fil', where fname is the root file name
C of the input file.
C
C Output file name: output.dat
C
C==================================================================
C
C VARIABLES USED BY THIS PROGRAM:
C
C ndi -- Number of direct components in stress/strain tensor.
C nshr -- Number of shear components in stress/strain tensor.
C ndip1 -- ndi + 1
C array -- Real array containing values read from results file.
C Equivalenced to jrray.
C jrray -- Integer array containing values read from results file.
C Equivalenced to array.
C fname -- Root file name of input file (w/o .fil extension).
C nru -- Number of results files (.fil) to be read.
C lrunit -- Array containing unit number and format of results files:
C lrunit(1,*) --> Unit number of input file.
C lrunit(2,*) --> Format of input file.
C loutf -- Format of output file:
C 0 --> Standard ASCII format.
C 1 --> ABAQUS results file ASCII format.
C 2 --> ABAQUS results file binary format.
C junit -- Unit number of file to be opened.
C jrcd -- Error check return code.
C .EQ. 0 --> No errors.
C .NE. 0 --> Errors detected.
C key -- Current record key identifier.
C jelnum -- Current element number.
C intpn -- Integration point number.
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
C and sets the value of NPRECD to 1 or 2, depending on whether
C the machine uses single or double precision. ABA_PARAM.INC is
C referenced from the \site (for NT) or /site (for Unix)
C ABAQUS release subdirectory when ABAQUS/Make is executed.
C
C==================================================================
C
include 'aba_param.inc'
parameter (melem=10000,mnode=10000,mkey=3000,mintpt=360)
dimension array(513),jrray(nprecd,513),lrunit(2,1)
equivalence (array(1),jrray(1,1))
C
C==================================================================
C
character fname*80
C
C SOME DIMENSION STATEMENTS MAY BE NEEDED IF YOU ARE DOING ADDITIONAL
C CALCULATIONS ON THE DATA
C
dimension coords(3,mintpt)
dimension iconn(3,melem),loc_el(melem)
dimension naxipt(melem),iglb_nod(mnode),loc_nod(mnode)
dimension islct(mkey),nintpt(melem),nodel(melem),iglb_el(melem)
dimension xyz(3,mnode),xpos(mnode),tang1(3),tang2(3)
dimension rot(3,3),xc(3),xpr(3,mintpt),scord(3,mintpt),xnorm(3)
dimension isctchk(melem),icrdchk(melem),ivarchk(melem)
C The following line places array 'var' on the stack. This array is very
C large (~40MB in size). This will cause problems on Windows, where the size
C of the stack is fixed and cannot be easily grown.
C
C dimension var(513,melem)
C
C It is advisable to place large data structures like this into
C dynamically-allocated memory, like so:
ALLOCATABLE ::var(:,:)
ALLOCATE(VAR(513,melem))
C
C==================================================================
C
C GET THE NAME OF THE RESULTS FILE.
C
C==================================================================
C
isect1 = 0
intpt1 = 0
icomp = 1
C
write(6,'(A$)') 'Enter the name of the input file (w/o .fil): '
read(5,'(A)') fname
C
100 write(6,*) 'Enter the postprocessing option:'
write(6,*) ' 1 - variation along the riser '
write(6,*) ' 2 - variation around the circumference'
write(6,*) ' of the elbow'
write(6,*) ' 3 - ovalization of elbow cross-section'
read(5,*) ipost
if (ipost .ne. 1 .and. ipost .ne. 2 .and.
& ipost .ne. 3) then
write(6,*)'Invalid entry - try again'
go to 100
endif
C
write(6,*)'Enter the first element number'
read(5,*) jel_1
C
if (ipost .ne. 1) then
jel_2 = jel_1
else
write(6,*)'Enter the last element number'
read(5,*) jel_2
endif
if (jel_2 .lt. jel_1) then
write(6,*)'last elem number must be greater than first'
return
endif
jdiff = jel_2 - jel_1 + 1
if (jdiff .gt. melem) then
write(6,*)'Too many elements - max is',melem
return
endif
C
if (ipost .lt. 3) then
write(6,*)'Enter the section point number'
read(5,*) isect1
endif
C
if (ipost .eq. 1) then
write(6,*)'Enter the integration point number'
read(5,*) intpt1
endif
C
if (ipost .lt. 3) then
998 write(6,999)
999 format(' Enter the key for the variable you wish to process:')
C
write(6,500)
500 format(2x,'TEMP ...... 2',/,
& 2x,'COORD...... 8',/,
& 2x,'S.......... 11',/,
& 2x,'SINV....... 12',/,
& 2x,'SF......... 13',/,
& 2x,'E.......... 21',/,
& 2x,'PE......... 22',/,
& 2x,'CE......... 23',/,
& 2x,'IE......... 24',/,
& 2x,'EE......... 25',/,
& 2x,'THE........ 88',/,
& 2x,'LE......... 89',/,
& 2x,'NE......... 90',/,
& 2x,'SP.........401',/,
& 2x,'EP.........403',/,
& 2x,'NEP........404',/,
& 2x,'LEP........405',/,
& 2x,'EEP........408',/,
& 2x,'IEP........409',/,
& 2x,'THEP.......410',/,
& 2x,'PEP........411',/,
& 2x,'CEP........412')
C
read(5,*) key
C
if (key .ne. 2 .and. key .ne. 8 .and.
& key .ne. 11 .and. key .ne. 12 .and.
& key .ne. 13 .and. key .ne. 21 .and.
& key .ne. 22 .and. key .ne. 23 .and.
& key .ne. 24 .and. key .ne. 25 .and.
& key .ne. 88 .and. key .ne. 89 .and.
& key .ne. 90 .and. key .ne. 401 .and.
& key .ne. 403 .and. key .ne. 404 .and.
& key .ne. 405 .and. key .ne. 408 .and.
& key .ne. 409 .and. key .ne. 410 .and.
& key .ne. 411 .and. key .ne. 412) then
C
write(6,*)'Invalid key - try again'
go to 998
C
endif
C
C Integration point coordinates available only at middle surface
C Overwrite the section point and set to the default
C
if (key .eq. 8) isect1 = 0
C
if (key .ge. 8) then
write(6,*)'Enter the attribute number'
read(5,*) icomp
endif
C
if (key .eq. 13) then
write(6,*)'Section force data written only once per element'
if (ipost .eq. 1) then
write(6,*)'Will use default integration & section point'
write(6,*)' '
intpt1 = 1
isect1 = 0
elseif (ipost .gt. 1) then
write(6,*)'Only use with option 1: Processing terminated'
return
endif
endif
C
endif
C
write(6,*)' Enter the step number'
read(5,*) istep1
write(6,*)' Enter the increment number'
read(5,*) iinc1
C
C SELECT THE RECORDS TO BE PROCESSED
C
C IF THE ISLCT ARRAY IS SET TO 1 THE DATA WILL BE PROCESSED
C IF THE ISLCT ARRAY IS SET TO 0 THE DATA WILL NOT BE PROCESSED
C
do ii = 1, mkey
islct(ii) = 0
end do
C
C ELEMENT DEFINITION
islct(1900)=0
C NODE DEFINITION
islct(1901)=0
C INCREMENT START
islct(2000)=0
C ELEMENT HEADER
islct(1)=0
C COORD
islct(8)=0
C
if (ipost .lt. 3) islct(key) = 1
keyprv=0
C
nels = 0
numnp = 0
C
itimchk = 0
ielchk = 0
do i = 1, melem
loc_el(i) = 0
isctchk(i) = 0
icrdchk(i) = 0
ivarchk(i) = 0
do j = 1, 513
var(j,i) = 0.0
end do
end do
C
C==================================================================
C
C OPEN THE OUTPUT FILE.
C
C==================================================================
open(unit=9,file='output.dat',status='unknown')
C
nru = 1
loutf = 0
lrunit(1,1) = 8
lrunit(2,1) = 2
C
call initpf(fname,nru,lrunit,loutf)
C
junit = 8
C
call dbrnu(junit)
C
C REWIND FILE
C
call dbfile(2,array,jrcd)
C
C
C==================================================================
C
C Read records from the results (.fil) file and process the data.
C Cover a maximum of 10 million records in the file.
C
C==================================================================
C
do 1000 k100 = 1, 100
do 1000 k1 = 1, 99999
C
C READ RECORD FROM .FIL FILE AND PROCESS DATA.
C
call dbfile(0, array, jrcd)
if (jrcd .ne. 0 .and. keyprv .eq. 2001) then
write(6,*)'End of file'
close(junit)
go to 1001
else if (jrcd .ne. 0) then
write(6,*)'Error reading input file'
return
endif
C
C LENGTH AND KEY RECORD NUMBER ARE ALWAYS NEEDED
C
lenf= jrray(1,1)
key = jrray(1,2)
C
C SPECIFIC RECORD NUMBERS FOLLOW...
C
C==================================================================
C
C ELEMENT DEFINITIONS
C
C==================================================================
if ( key .eq. 1900 ) then
jelnum=jrray(1,3)
if (jelnum .gt. melem) then
write(6,*)'Maximum global elem number too large'
write(6,*)'Increase melem'
return
endif
ctype= array(4)
if (jelnum .ge. jel_1 .and.
& jelnum .le. jel_2) then
ielchk = 1
nels = nels + 1
jel = nels
iglb_el(jel) = jelnum
loc_el(jelnum) = jel
nodel(jel) = lenf - 4
naxipt(jel) = 1
if (ctype .ne. 'ELBOW31' .and.
& ctype .ne. 'ELBOW31B' .and.
& ctype .ne. 'ELBOW31C' .and.
& ctype .ne. 'ELBOW32') then
write(6,*)'Invalid element type encountered'
return
endif
if (ctype .eq. 'ELBOW32')naxipt(jel) = 2
do kk=5,lenf
ii=kk-4
iconn(ii,jel)=jrray(1,kk)
end do
endif
keyprv = key
goto 9099
C
C==================================================================
C
C NODE DEFINITIONS
C
C==================================================================
else if ( key .eq. 1901 ) then
jnod = jrray(1,3)
if (jnod .gt. mnode) then
write(6,*)'Maximum global node number exceeded'
write(6,*)'Increase mnode'
return
endif
numnp = numnp + 1
iglb_nod(numnp) = jnod
loc_nod(jnod) = numnp
xyz(1,numnp) = array(4)
xyz(2,numnp) = array(5)
xyz(3,numnp) = array(6)
keyprv = key
goto 9099
C
C==================================================================
C
C CURRENT STEP AND INCREMENT NUMBER.
C
C==================================================================
else if ( key .eq. 2000 ) then
ttime=array(3)
stime=array(4)
sfreq=array(12)
istep=jrray(1,8)
iinc =jrray(1,9)
C
C SET THE FLAG ITIME IF THIS IS THE REQUESTED STEP/INC
C
if (istep .eq. istep1 .and.
& iinc .eq. iinc1) then
itime = 1
itimchk = 1
else
itime = 0
endif
c
keyprv = key
goto 9099
C
C==================================================================
C
C ELEMENT VARIABLES AT INTEGRATION POINTS WITHIN THE ELEMENT
C
C==================================================================
C
C ELEMENT DATA
C
else if ( key .eq. 1 ) then
C
C PROCESS THIS RECORD IF THIS IS THE CORRECT STEP/TIME
C
if (itime .eq. 1) then
jelnum = jrray(1,3)
jel = loc_el(jelnum)
ielem = 0
C
C PROCESS THIS ELEMENT IF IT IS IN THE LIST OF DESIRED ELEMS
C
if (jel .gt. 0) then
ielem = 1
intpn = jrray(1,4)
if (intpn .gt. mintpt) then
write(6,*)'Max number of int points exceeded'
write(6,*)'Increase mintpt'
return
endif
isect = jrray(1,5)
ilocn = jrray(1,6)
if (ilocn .ne. 0) then
write(6,*)'Element data must be at the int points'
return
endif
crbar = array(7)
ndi = jrray(1,8)
nshr = jrray(1,9)
ndir = jrray(1,10)
nsfc = jrray(1,11)
if (isect .eq. isect1)then
isctflg = 1
isctchk(jel) = 1
else
isctflg = 0
endif
if ( (ipost .eq.1 .and. intpn .eq. intpt1) .or.
& ipost .gt. 1) then
iptflg = 1
else
iptflg = 0
endif
endif
endif
keyprv = key
goto 9099
C
C==================================================================
C
C COORDINATES
C
C==================================================================
else if ( key .eq. 8 .and. ipost .gt. 1) then
C
C PROCESS IF THIS IS THE CORRECT STEP/INC, ELEMENT, SECTION
C POINT, AND INTEGRATION POINT
C
icheck = itime*ielem*iptflg
if (icheck .eq. 1 ) then
c if (ipost .eq. 3) nintpt(jel) = intpn
icrdchk(jel) = 1
if (ipost .gt. 1) nintpt(jel) = intpn
do kk=3,lenf
ii=kk-2
coords(ii,intpn)=array(kk)
end do
C
if (islct(key) .eq. 1 .and. ipost .lt. 3) then
ivarchk(jel) = 1
jvar = intpn
nintpt(jel) = intpn
do kk=3,lenf
ii=kk-2
var(ii,jvar)=array(kk)
end do
if (icomp .gt. ii) then
write(6,*)'Invalid component'
return
endif
endif
C
endif
keyprv = key
goto 9099
C
C==================================================================
C
C VARIABLE RECORD
C
C==================================================================
C PROCESS IF THIS IS THE DESIRED VARIABLE
else if (islct(key) .eq. 1 .and. ipost .lt. 3) then
C
C CONTINUE PROCESSING IF THIS IS THE CORRECT STEP/INC, ELEMENT,
C SECTION POINT, AND INTEGRATION POINT
C
icheck = itime*ielem*isctflg*iptflg
if (icheck .eq. 1 ) then
ivarchk(jel) = 1
if (ipost .eq. 1) then
jvar = jel
else
jvar = intpn
endif
nintpt(jel) = intpn
do kk=3,lenf
ii=kk-2
var(ii,jvar)=array(kk)
end do
if (icomp .gt. ii) then
write(6,*)'Invalid component'
return
endif
endif
keyprv = key
goto 9099
C
C==================================================================
C
C ANYTHING ELSE
C
C==================================================================
else
c write (6,9000) key
9000 format(I5,' *** NO CODING FOR THIS RECORD ***')
keyprv = key
9099 end if
C
C==================================================================
C
C END LOOPS
C
C==================================================================
1000 continue
1001 continue
C
C==================================================================
C POSTPROCESS DATA HERE
C==================================================================
do ii = 1, nels
nintpt(ii) = nintpt(ii)/naxipt(ii)
end do
C
C==================================================================
C VERIFY THAT APPROPRIATE DATA WERE WRITTEN TO RESULTS FILE
C==================================================================
c iquit = 0
if (ielchk .eq. 0) then
write(6,*)'Desired element not found'
return
endif
C
if (itimchk .eq. 0) then
write(6,*)'Desired step/increment not found'
return
endif
C
do iel = 1, nels
if (isctchk(iel) .eq. 0 .and. ipost .lt. 3) then
write(6,*)'Desired element and/or section point not found'
return
endif
if (icrdchk(iel) .eq. 0 .and. ipost .gt. 1) then
write(6,*)'Integration point coords not found'
return
endif
if (ivarchk(iel) .eq. 0 .and. ipost .lt. 3) then
write(6,*)'Desired element variable not found'
return
endif
end do
C
C****************************************
C
C PLOT VARIABLE ALONG ELBOW LENGTH
C
C****************************************
C
if (ipost .eq. 1) then
C
nn_old = iconn(nodel(1),1)
do iel = 1, nels
nn_new = iconn(1,iel)
if (iel .gt. 1 .and. nn_new .ne. nn_old) then
write(6,*)'Error in element connectivity'
return
endif
nnum_b = iconn(1,iel)
nnum_b = loc_nod(nnum_b)
nnum_e = iconn(nodel(iel),iel)
nnum_e = loc_nod(nnum_e)
if (iel .eq. 1) xpos(nnum_b) = 0.0
xlength = 0.d+0
do jj = 1, 3
delx = xyz(jj,nnum_e) - xyz(jj,nnum_b)
xlength = xlength + delx*delx
end do
xlength = sqrt(xlength)
xpos(nnum_e) = xpos(nnum_b) + xlength
x_mid = 0.5*(xpos(nnum_e) + xpos(nnum_b))
C
C EXTRACT VARIABLE INFORMATION AND WRITE TO FILE
C
data = var(icomp,iel)
write(9,2000)x_mid,data
2000 format(5x,e17.9,5x,e17.9)
C
nn_old=iconn(nodel(iel),iel)
C
end do
return
C
C
C****************************************
C
C PLOT VARIABLE AROUND CIRCUMFERENCE OF ELBOW
C
C****************************************
C
else if (ipost .eq. 2) then
C
xpos(1) = 0.0
C
locel = loc_el(jel_1)
iaxial = 1
C
if (naxipt(locel) .eq. 2) then
3000 write(6,*)'Enter the axial station'
read(5,*)iaxial
if (iaxial .ne. 1 .and. iaxial .ne. 2) then
write(6,*)'Try again - station must be 1 or 2'
go to 3000
endif
endif
C
do ipt = 1 , nintpt(locel)
jpt = ipt + (iaxial - 1)*nintpt(locel)
if (ipt .gt. 1) then
xlength = 0.0
do jj = 1, 3
delx = coords(jj,jpt) - coords(jj,jpt-1)
xlength = xlength + delx*delx
end do
xlength = sqrt(xlength)
xpos(ipt) = xpos(ipt-1) + xlength
endif
C
C EXTRACT VARIABLE INFORMATION AND WRITE TO FILE
C
data = var(icomp,jpt)
write(9,2000)xpos(ipt),data
C
end do
ipt = 1
jpt = ipt + (iaxial - 1)*nintpt(locel)
data = var(icomp,jpt)
xfinal = xpos(nintpt(locel)) + (xpos(2)- xpos(1))
write(9,2000)xfinal,data
return
C
C****************************************
C
C OVALIZATION PLOT
C
C****************************************
C
else if (ipost .eq. 3) then
C
locel = loc_el(jel_1)
iaxial = 1
C
if (naxipt(locel) .eq. 2) then
3500 write(6,*)'Enter the axial station'
read(5,*)iaxial
if (iaxial .ne. 1 .and. iaxial .ne. 2) then
write(6,*)'Try again - station must be 1 or 2'
go to 3500
endif
endif
C
C AVERAGE THE COORDS OF THE INTEGRATION POINTS TO GET THE CENTER
C
do i = 1, 3
xc(i) = 0.0
end do
C
xnint = real(nintpt(locel))
C
do ipt = 1 , nintpt(locel)
jpt = ipt + (iaxial - 1)*nintpt(locel)
C
do i = 1, 3
xc(i) = xc(i) + coords(i,jpt)/xnint
end do
C
end do
C
C SHIFT THE COORDINATES SO THAT CENTER OF SECTION IS AT ORIGIN
C
do ipt = 1 , nintpt(locel)
jpt = ipt + (iaxial - 1)*nintpt(locel)
do i = 1, 3
scord(i,jpt) = coords(i,jpt) - xc(i)
end do
C
end do
C
C DETERMINE TANGENT VECTOR 1 (FROM CENTER OF CROSS-SECTION
C TO INT PT 1)
C
jpt = 1 + (iaxial - 1)*nintpt(locel)
xmag = 0.0
do i = 1, 3
tang1(i) = coords(i,jpt) - xc(i)
xmag = xmag + tang1(i)*tang1(i)
end do
xmag = sqrt(xmag)
do i = 1, 3
tang1(i) = tang1(i)/xmag
end do
C
C GUESS TANGENT VECTOR 2 (FROM CENTER OF CROSS-SECTION
C TO INT PT 2)
C
jpt = 2 + (iaxial - 1)*nintpt(locel)
xmag = 0.0
do i = 1, 3
tang2(i) = coords(i,jpt) - xc(i)
xmag = xmag + tang2(i)*tang2(i)
end do
xmag = sqrt(xmag)
do i = 1, 3
tang2(i) = tang2(i)/xmag
end do
C
C DETERMINE THE NORMAL VECTOR, TANG1 X TANG2
C
xnorm(1) = tang1(2)*tang2(3) - tang1(3)*tang2(2)
xnorm(2) = - tang1(1)*tang2(3) + tang1(3)*tang2(1)
xnorm(3) = tang1(1)*tang2(2) - tang1(2)*tang2(1)
xmag = 0.0
do i = 1, 3
xmag = xmag + xnorm(i)*xnorm(i)
end do
xmag = sqrt(xmag)
do i = 1, 3
xnorm(i) = xnorm(i)/xmag
end do
C
C REDEFINE TANGENT VECTOR 2
C
tang2(1) = - tang1(2)*xnorm(3) + tang1(3)*xnorm(2)
tang2(2) = tang1(1)*xnorm(3) - tang1(3)*xnorm(1)
tang2(3) = - tang1(1)*xnorm(2) + tang1(2)*xnorm(1)
C
C DETERMINE THE ROTATION MATRIX
C
do i = 1, 3
rot(1,i) = tang1(i)
rot(2,i) = tang2(i)
rot(3,i) = xnorm(i)
end do
C
C TRANSFORM COORDINATES INTO LOCAL SYSTEM
C
do ipt = 1 , nintpt(locel)
jpt = ipt + (iaxial - 1)*nintpt(locel)
do i = 1, 3
xpr(i,ipt) = 0.0
do j = 1, 3
xpr(i,ipt) = xpr(i,ipt) + rot(i,j)*scord(j,jpt)
end do
end do
end do
C
C OUTPUT COORDS TO FILE
C
do ipt = 1, nintpt(locel)
write(9,4000)(xpr(i,ipt),i=1,2)
4000 format(5x,e17.9,5x,e17.9,5x,e17.9)
end do
C
ipt = 1
write(9,4000)(xpr(i,ipt),i=1,2)
endif
C
close (unit=9)
C
IF( ALLOCATED(var) ) DEALLOCATE(var)
return
end