C smallCycleFilter.f vers. 1.62 FAC Sep.28 2014 SAVE C Compile: gfortran -g -w -fbounds-check smallCycleFilter.f -o smallCycleFilter C Usage: smallCycleFilter scaleFactor peakError outputFile C scaleFactor is the multiplier for the history C peakError is the #PKERROR= window from flex10.env C Input data is one value per line. C Subroutine getPeaks() operates on list in COMMON memory. C Program eliminates small cycles that are less than a certain voltage in size(range) C Main function getPeakLoads() adapted from plateWeldflaw.f C Copyright (C) 2013 Al Conle C This file is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the license, or (at C your option) any later version. C This file is distributed in the hope that it will be useful, but C WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTA- C BILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public C License for more details. C You should have received a copy of the GNU General PUblic License along C with this program; if not, write to the Free Software Foundation, Inc., C 59 Temple Place -Suite 330, Boston, MA 02111-1307, USA. Try also their C web site: http://www.gnu.org/copyleft/gpl.html C Note that some subroutines included in this work are from GNU GPL licenced C program: http://fde.uwaterloo.ca/Fde/Calcs/saefcalc1.html C vers. 1.62 Correct iloadflag having no dimensions during read. C vers. 1.61 When no point eliminated only print out history. Apr01_2014 C vers. 1.6 Add code to detect the alignment between end of a block and C begin of next block C vers. 1.5 Abandon Delete pt. list for simple flag 1 or 0 Oct 24 2013 C Gets rid of "sort" of bad point numbers. Uses a flag system. C vers. 1.0 Abandon getPeakLoads() from plateWeldflaw.f. Try new method. C vers. 0.9 Initial adaption from plateWeldflaw.f subroutine getPeakLoads() C character*300 inp300,jnp300,Cwebpage ! used to read in lines as chars. character*1 inpone(300),jnpone(300) character*10 inpten(30) equivalence (inpone(1), inpten(1), inp300) equivalence (jnpone(1),jnp300) character*80 argv integer*4 iargc,nargc c Load history storage. If changing dimensions, also change same C in S/R getPeaks() (1 column of data only) C integer*4 sae(10000) C Stress membrane, bending, and total storage, and iloadflag() C iloadflag() = 1 if load ok, 0 not used, or intValue if repeat loads: real*4 ststot(1000000) !see also s/r getPeaks() storages integer*4 iloadflag(1000000) C Note! if you change above dimensions, then also change in s/r real*4 ststotmax,ststotmin,ststotwindow integer*4 nloads logical debugLoads common /LOADS/ ststot, iloadflag,ststotmax,ststotmin,ststotwindow, & nloads,debugLoads debugLoads=.false. ! use to debug load read and peakpick section. C Storage area limits. Change this if you change above real* and int* maxloads=1000000 ! max load storage C Test if -Wuninitilized warning works in gfortran: C i=ifix(x) C Works! C Initilize some variables. C--------------------------- Run time input data------------------ 184 continue write(6,185) write(0,185) 185 format("#smallCycleFilter.f vers. 1.62"/ & "#Usage: smallCycleFilter scale filter outfile"/) nargc = iargc() C Note that in HP fortran the arg numbers must be changed by -1 C and that iargc probably includes the "plateWeldflaw" as an arg. if( nargc .ne. 2)then write(0,*)" smallCycleFilter: usage ERROR" write(6,*)" smallCycleFilter: usage ERROR" write(0,186) write(6,186) 186 format( & "#Usage: smallCycleFilter Scale Filter outputFile"// & "# Where Scale is the multiplier applied to all history pts,"// & "# and Filter is the voltage range of small cycles to be", & " removed"// & "# Stopping now." & ) stop endif C The first arg is the history multiplication factor jvect=1 call getarg(jvect,argv) read(argv,*,err= 178)scaleValue write(6,*)"#scaleValue= ",scaleValue if(scaleValue .eq. 0.0)go to 178 jvect=jvect+1 call getarg(jvect,argv) read(argv,*,err=179)voltfilter write(0,*)"#Filter= ",voltfilter," volts" write(6,*)"#Filter= ",voltfilter," volts" C Both scaleValue and voltfilter have been read in successfully go to 190 C Bad arguments in command line 178 write(0,*)"# ERROR: bad scaleValue argument=",argv stop 179 write(0,*)"# ERROR: bad filter argument= ",argv stop 190 continue ststotwindow=voltfilter/scaleValue write(6,191)ststotwindow write(0,191)ststotwindow 191 format("#Convert filt to input units: ", & "Filter/ScaleFactor = ",e14.7) C------------------Get Load history file------------------------------ C Expected format: C #Comment 1 C #Comment line 2 C C # C 0.00 C 100.00 C -100.00 C ....etc C #END= 0 write(0,901) write(6,901) 901 format("# Getting Load history file from std.input ") C open(unit=10,file=histfile) ninput=0 nloads=0 900 continue !Loop back to here for next input line. read(5,"(a300)",end=950)inp300 ninput=ninput+1 C write(0,*)ninput if(inp300.eq." ")then ! Check for blank line C write(6,"(a1)")" " go to 900 endif if(inpone(1).ne."#")then C We may have a data line, or someone screwed up and put the # later in line. C See if 1st char is a # later in line do 905 i=1,300 if(inpone(i).eq." ") go to 905 C No? 1st non-blank found, if # its not a data line loc=i if(inpone(i).eq."#") then C Shift the whole mess over to begining of field inew=0 do 903 j=loc,300 inew=inew+1 inpone(inew)=inpone(j) 903 continue Cdebug write(0,*)"Shifted a comment: ",(inpone(j),j=1,inew) C Ok, now its a nice comment. Go play with it. go to 930 else C The first non blank is not a # go to 910 endif 905 continue go to 900 !assume garbage in line. 910 continue C 1st non blank is not a # Check if its a number. if(inpone(loc).ne."+" .and. & inpone(loc).ne."-" .and. & inpone(loc).ne."." .and. & inpone(loc).ne."0" .and. & inpone(loc).ne."1" .and. & inpone(loc).ne."2" .and. & inpone(loc).ne."3" .and. & inpone(loc).ne."4" .and. & inpone(loc).ne."5" .and. & inpone(loc).ne."6" .and. & inpone(loc).ne."7" .and. & inpone(loc).ne."8" .and. & inpone(loc).ne."9" )then C It must be a letter, not a number. Time to bomb out. write(0,915)ninput write(6,915)ninput 915 format(" ERROR stdin: input line no. ",I5, & " not a # and not a number. Fix data in load history file") stop endif C Ok, its a number nloads=nloads+1 if(nloads .gt. maxloads)then write(0,916)nloads write(6,916)nloads 916 format(" Error too many data points:",i5, & " recompile plateWeldflaw.f or reduce load history data") stop endif read(inp300,*,err=925)ststot1 if(nloads .eq. 1)then ststotmax=ststot1 ststotmin=ststot1 endif if(ststot1 .gt. ststotmax) ststotmax=ststot1 if(ststot1 .lt. ststotmin) ststotmin=ststot1 ststot(nloads)=ststot1 iloadflag(nloads) =1 ! in other programs this is read in. C write(6,920)ststot(nloads),nloads C 920 format("#history.org ",f7.1,1x,i6) go to 900 925 write(0,926)nloads,inp300 write(6,926)nloads,inp300 926 format("#READ ERROR: history file: data line no.= ",i9/ & "#LINE= ",a300/ "# Stopping now.") stop endif 930 continue C Yes, first char was a # C All comment lines are just written out, or deleted go to 900 950 continue ! end of file comes here -------------------------- if(nloads.eq.0)then write(6,951) write(0,951) 951 format("#Error: histfile: contains no data !") stop endif C All data is in. Close file C close(unit=10) write(6,952)nloads,ststotmax,ststotmin write(0,952)nloads,ststotmax,ststotmin 952 format("#history #Data input completed. nloads= ",i5/ & "#history #Max= ",f7.1/"#history #Min= ",f7.1/ & //) write(0,954)ststotwindow write(6,954)ststotwindow 954 format("# Eliminating non-reversal points from history..."/ &"#history #Elimination window = ",f7.2/ &"#history #Any 1/2 cylce smaller than this will be eliminated...") ivec=1 iret=0 call getPeaks(ivec,iret) if(iret.ne.0)then write(0,*)"#ERROR: getPeaks() iret=",iret," Stopping." stop endif 9000 continue 9005 continue stop end C====================================================================== SUBROUTINE getPeaks(ivec,iret) SAVE C ivec= function instruction, usually=1 C iret= return value C s/r to look at the total stress and remove any non-reversal points C from the load history. C Stress membrane, bending, and total storage: real*4 ststot(1000000) integer*4 iloadflag(1000000) ! list of points to be removed in 2nd pass of data C Use as logic flag. 1=keep, 0=delete real*4 ststotmax,ststotmin,ststotwindow real*4 stsold ! the last reversal point real*4 stshere ! the present position real*4 stsnext ! the next point being examined integer*4 nloads logical debugLoads common /LOADS/ ststot,iloadflag, ststotmax,ststotmin,ststotwindow, & nloads,debugLoads ngood=0 ! records the no. of good points n=0 ! points to the data being examined nbad=0 ! counts the no. of points being eliminated C Get the first point. n=n+1 stsold=ststot(1) nstsold=1 ngood=ngood+1 C check if the first and last point in the history are same if(ststot(1).eq.ststot(nloads) )then write(0,90) write(6,90) 90 format("#history #First and last pts in history are same."/ & "#history #Eliminating the 1st history point now") nbad=nbad+1 iloadflag(1)=0 write(6,*)"#Delete firstPoint: ",n,ststot(1) C Save a copy for rainflow file created near end of program ststot1=ststot(1) ! Not used. Adds an extra ramp per block! endif C Get the second point stspres to create the initial line. 100 continue n=n+1 stspres=ststot(n) !get the 2nd point nstspres=n if(stspres .eq. stsold)then ! remove the new point nbad=nbad+1 iloadflag(nstspres)=0 write(6,*)"#DeleteLd samePoint: ",nstspres,stspres go to 100 endif C 1st (stsold) and 2nd pts (stspres) have been read in---------------- C 2nd pt is not equal to 1st. Establish direction iupdown=+1 if(stspres .lt. stsold)iupdown=-1 C stspres is acceptable if(debugLoads)write(6,*)"#stsold,nstsold,stspres,nstspres", & stsold,nstsold,stspres,nstspres,iupdown C ok, first point and direction is set. Keep going 2000 continue n=n+1 if(n .gt. nloads)go to 9000 stsnext=ststot(n) nstsnext=n if(debugLoads)write(6,*)"#stsnext,nstsnext: ",stsnext,nstsnext, & iupdown C If its equal to the present one we can eliminate the next pt if(stsnext .eq. stspres)then nbad=nbad+1 iloadflag(nstsnext)=0 write(6,*)"#DeleteLd samePoint: ",nstsnext stsnext go to 2000 endif C No, its not equal. Could be a reversal or a continuation of ramp. if(stsnext .lt. stspres)go to 3000 ! Potentially it is a reversal point, downwards C If not, then things are going in the same direction if iupdown =1 C New is going upwards-------------------------------------------------- C If the same direction we can discard stspres if(iupdown .eq. +1)then C Same direction but further, discard pres pt if(debugLoads)write(6,*)"#discarding pres pt no.", & nstspres," iud= ",iupdown nbad=nbad+1 iloadflag(nstspres)=0 write(6,*)"#Delete non-revPoint: ",nstspres,stspres stspres=stsnext nstspres=nstsnext go to 2000 ! get next point endif C iupdown must have been = -1, thus previous direction was in compression. C If we get to this point then it is a potential reversal. C See if its size is bigger than the allowable window if((stsnext-stspres) .lt. ststotwindow)then ! too small, skip new pt nbad=nbad+1 iloadflag(nstsnext)=0 write(6,*)"#Delete tooSmallPoint: ",nstsnext,stsnext go to 2000 !go get a new one endif C Ok, it is bigger than the window. Must be a reversal if(debugLoads)write(6,*)"#Rev. C_T: old, pres, next ", & stsold,nstsold,stspres,nstspres,stsnext,nstsnext ngood=ngood+1 iupdown=+1 !change direction C Set up the reversal point as old stsold=stspres nstsold=nstspres stspres=stsnext nstspres=nstsnext go to 2000 C New is going downwards ------------------------------------------------------- 3000 continue C the new point is going downwards. Check the old direction if(iupdown .eq. -1)then ! same direction C Discard the old point, and replace it with the new pt. nbad=nbad+1 iloadflag(nstspres)=0 write(6,*)"#Delete non-revPoint: ",nstspres,stspres stspres=stsnext nstspres=nstsnext go to 2000 endif C iupdown was +1, previous direction was into tension C Potential change in direction. Is it bigger than the window if((stspres-stsnext) .lt. ststotwindow)then C too small, eliminate the new point nbad=nbad+1 iloadflag(nstsnext)=0 write(6,*)"#Delete tooSmallPoint: ",nstsnext,stsnext go to 2000 endif C Ok, reversal is bigger than the window. Its a reversal ngood=ngood+1 if(debugLoads)write(6,*)"#Rev.T_C: old, pres, next ", & stsold,nstsold,stspres,nstspres,stsnext,nstsnext iupdown=-1 stsold=stspres nstsold=nstspres stspres=stsnext nstspres=nstsnext go to 2000 9000 continue if(nbad.eq.0)then ! no points eliminated write(0,9002) write(6,9002) 9002 format("#history #No points needed to be eliminated."/) ngood=nloads go to 9820 ! print out the history endif write(0,9005)ngood,nbad write(6,9005)ngood,nbad 9005 format("#history #PPICK found ",i7,"good pts, and ",i7, & " to be eliminated. See list above.") C ------------------selction is done. Must now eliminate and update list 9400 continue C Now go throught the stress storage and take out those to be eliminated C Start at the first point to be eliminated iput=0 ! index points to next available (empty) storage point 9500 continue iput=iput+1 if(iloadflag(iput) .ne. 0)goto 9500 C write(6,*)"#empty at= ",iput !debug C yes flag is 0, this spot is open for a good pt. C Now find the next pt that needs to be saved jnextsavept=iput 9520 continue jnextsavept=jnextsavept+1 if(jnextsavept .gt. nloads)goto 9600 ! end of data found if(iloadflag(jnextsavept) .eq. 0)goto 9520 C write(6,*)"#nextFull at= ",jnextsavept !debug C Yes, we have the next point with iloadflag=1 C Move this next point to the available postition at iput ststot(iput)=ststot(jnextsavept) iloadflag(iput)=iloadflag(jnextsavept) iloadflag(jnextsavept)=0 !set the moved pt flag to 0 C For Pdprop crack propagation programs: C write(6,9590)XXX,YYY,ststot(i),iloadflag(iput),iput C 9590 format("#history #Filtered out",3(f7.1,1x),i6,1x,i6) goto 9500 9600 continue C End of data encountered. C Right now iput is still sitting on a flag=0 pt. iput=iput-1 ! point back to last non-zero flag point. C End-Begin Alignment Check----------------------------------- C Check if the first and last point in the history are same C Do this after filtering nstart=1 ! This is starting point of history, Change? if(ststot(1).eq.ststot(nloads) )then nstart=2 write(0,9604)ststot(1),ststot(nstart) write(6,9604)ststot(1),ststot(nstart) 9604 format("#history #1st and Last pts in Filtered hist are same."/ & "#history # ststot(1)= ",f8.2/ & "#history #Eliminating the 1st history point now"/ & "#history #New start pt: ststot()= ",f8.2) C Save a copy for rainflow file created near end of program ststot1=ststot(2) ! Not used. Adds an extra ramp per block! endif C There is still the problem of aligning the begin/end of the C filtered history. The PDlisting process does not like two C points in the same direction, one after the other. C The method is to look at the end slope, the joint slope, and C the begin slope. If any of these three have "two in a row" C with the same sign, we have an alignment problem. iwrap=1 ! =1 means Ok, =0 means problem iendDir= +1 if( (ststot(iput) - ststot(iput-1)) .lt. 0.)iendDir=-1 ijoinDir= +1 if( (ststot(nstart) - ststot(iput)) .lt. 0.)ijoinDir=-1 ibeginDir= +1 if( (ststot(nstart+1) - ststot(nstart)) .lt. 0.)ibeginDir=-1 write(0,*)"#ststot(nstart): ",ststot(nstart), nstart write(0,*)"#ststot(nstart+1): ",ststot(nstart+1), nstart+1 write(0,*)"#ststot(iput-1): ",ststot(iput-1), iput-1 write(0,*)"#ststot(iput): ",ststot(iput), iput write(0,*)"#iendDir: ",iendDir write(0,*)"#ijoinDir: ",ijoinDir write(0,*)"#ibeginDir: ",ibeginDir if(iendDir .eq. ijoinDir)iwrap=0 if(ijoinDir .eq. ibeginDir)iwrap=0 if(iwrap .eq. 0)then write(0,9608)ststot(iput-1),ststot(iput),ststot(nstart), & ststot(nstart+1) write(6,9608)ststot(iput-1),ststot(iput),ststot(nstart), & ststot(nstart+1) 9608 format("#history #Warning: After filtering: End of history and", & "Begining have an alignment problem:"/ & "#history # Last-1 Pt: ",f8.2/ & "#history # Last Pt: ",f8.2/ & "#history # 1st Pt: ",f8.2/ & "#history # 2nd Pt: ",f8.2) C Check if all three same direction if((iendDir.eq.ijoinDir).and.(ijointDir.eq.ibeginDir) )then write(0,9612) write(6,9612) 9612 format("#history # All in same direction. Eliminate 1st and", & " last points.") nstart=nstart +1 ! get rid of 1st pt. iput=iput-1 ! get rid of last pt. goto 9620 endif if(iendDir .eq. ijoinDir)then write(0,9616) write(6,9616) 9616 format("#history # Last and Join in same direction. ", & "Eliminating the Last point.") iput=iput-1 goto 9620 endif if(ijoinDir .eq. ibeginDir)then write(0,9618) write(6,9618) 9618 format("#history # Join and Begin are same direction.", & "Eliminating the First point.") nstart=nstart+1 goto 9620 endif endif C Thinking about this some more: It would probably achieve the C same effect by adding the 1st point onto the end of the pre-filter C history and then filtering. (too late now :) C Perhaps it is better to check "properly" as above anyway. 9620 continue write(0,9630)ngood,iput write(6,9630)ngood,iput 9630 format("#history #BugCheck: GoodPtsFound=",i8, & " GoodPtsSaved=",i8," should be nearly equal?.") C When nstart is not equal to 1 anymore we need to shift the C history such that the start point, presently = nstart, is pt 1 C We could probably skip this in the stand alone version. C Check results with a full print icount=0 do 9800 i=nstart,iput icount=icount+1 C write(6,9790)stsm(i),stsb(i),ststot(i),iloadflag(i),icount C 9790 format("#history #Filteredck ",3(f7.1,1x),i8,1x,i6) if(nstart .ne. 1)then C Shift the history to compensate for nstart > 1 ststot(icount) =ststot(i) iloadflag(icount)=iloadflag(i) endif 9800 continue C If history has been shifted for non-zero nstart we need nstart=1 ! if it was already =1 then no change. write(0,9610)ngood,icount write(6,9610)ngood,icount 9610 format("#history #BugCheck: GoodPtsFound=",i8, & " GoodPtsSaved=",i8," should be nearly equal.") 9820 continue C Stand alone version prints here. nstart now = 1. do 9850 i=1,icount write(6,9840)ststot(i) 9840 format(f12.5) 9850 continue write(6,9860) 9860 format("#END= 0") write(0,9861) 9861 format("#Warning: #END= 0 added to end of file") nloads=icount return end C================================================================================== SUBROUTINE subCutBlanks(idevice,inp1,idimen) SAVE C subCutBlanks.f s/r to remove the leading and trailing blanks for output C Used mainly to print long comment lines C Usage: call subCutBlanks(idevice,inp400,idimen) C idevice= output device no. in fortran C inp400= char storage vector C idimen= length of incoming char string C (or how many to look at ?) C Copyright (C) 2013 Al Conle C This file is free software; you can redistribute it and/or modify C it under the terms of the GNU General Public License as published by C the Free Software Foundation; either version 2 of the license, or (at C your option) any later version. C This file is distributed in the hope that it will be useful, but C WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTA- C BILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public C License for more details. C You should have received a copy of the GNU General PUblic License along C with this program; if not, write to the Free Software Foundation, Inc., C 59 Temple Place -Suite 330, Boston, MA 02111-1307, USA. Try also their C web site: http://www.gnu.org/copyleft/gpl.html integer*4 idevice,idim, ncharmax/400/ integer*4 ifirst,ilast character*1 inp1(400) ! if you change this, change also ncharmax above C Check for programming error: if(idimen .gt. ncharmax)then write(0,100)idimen,ncharmax 100 format("#Error: s/r subCutBlanks(): dimension of char varible= ", & i10," bigger than max allowed =",i5," stopping...") stop endif if(idimen .lt. 1 )then write(0,100)idimen 101 format("#Error: s/r subCutBlanks(): dimension of char varible= ", & " is less than 1, stopping...") stop endif C Find the first non-blank character do 200 i=1,ncharmax if(inp1(i) .eq. " ")goto 200 ifirst=1 goto 210 200 continue C if we end up here there are no chars within inp(1) and inp(idimen) C Just write out one blank for this line write(idevice,205) 205 format(" ") return 210 continue ! now hunt for the last char by starting at idimen and C working backwards. do 300 i=idimen,1,-1 if(inp1(i) .eq. " ")goto 300 ilast=i ! No?, then we have found an non blank go to 310 300 continue C If we end up here, the whole line was blanks, write out one blank write(idevice,205) return 310 continue ! ifirst and ilast are set, write out the stuff write(idevice,315)(inp1(i),i=ifirst,ilast) 315 format(400a1) return end