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 var(513,melem),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 C================================================================== C C GET THE NAME OF THE RESULTS FILE. C C================================================================== C isect1 = 0 intpt1 = 0 icomp = 1 C write(6,*) '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 return end