PROGRAM rboef
!-------------------------------------------------------------------------
! Beams on random elastic foundations
! Used in the paper submitted to Geotechnique in October 2010
! "Reliability analysis of beams on random elastic foundations"
! by D.V. Griffiths, J. Paiboon, J. Huang and G.A. Fenton
!--------------------------------------------------------------------------
 USE main
 USE geom
 USE gaf95
 IMPLICIT NONE
 INTEGER,PARAMETER::iwp=SELECTED_REAL_KIND(15)
 INTEGER::fixed_freedoms,i,iel,k,loaded_nodes,ndof=4,nels,neq,nod=2,      &
   nodof=2,nn,nprops=2,np_types,nr,nlen,j,maxfld,nrv,kseed,nsim,          &
   jr,nseed,iterm,nrfy
 CHARACTER(LEN=80)::trash,argv
 REAL(iwp)::penalty=1.0e20_iwp,zero=0.0_iwp,h1,total=0.0_iwp,             &
   average=0.0_iwp
 REAL thy,yl,onn,zro,for,hlf,threept5,eh
 LOGICAL debug,lxfld,lsquare
 CHARACTER varfnc*6
!-----------------------dynamic arrays-------------------------------------
 INTEGER,ALLOCATABLE::etype(:),g(:),g_g(:,:),kdiag(:),nf(:,:),no(:),      &
   node(:),num(:),sense(:)
 REAL(iwp),ALLOCATABLE::action(:),eld(:),ell(:),km(:,:),kv(:),loads(:),   &
   mm(:,:),prop(:,:),value(:),loadsi(:),trans(:),rot(:),dtots(:),         &
   stord(:,:),mn(:),sd(:)
 REAL, ALLOCATABLE::a(:,:),R(:,:),gfld(:,:),ac(:,:,:),ktots(:)
 data zro/0.0/, hlf/0.5/, onn/1.0/, threept5/3.5/
 iterm=33
 lsquare=.FALSE.
 debug=.FALSE.
 lxfld=.FALSE.
!-----------------------input and initialisation--------------------------
 CALL getname(argv,nlen)
 OPEN(10,FILE=argv(1:nlen)//'.dat')
 OPEN(11,FILE=argv(1:nlen)//'.res')
 READ(10,*)trash
 WRITE(11,'(A)')trash
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)h1
 WRITE(11,'(A)')"Physical length of soil to be modelled"
 WRITE(11,'(f8.2)')h1
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)nels
 WRITE(11,'(A)')"Number of finite elements in model"
 WRITE(11,'(I5)')nels
 WRITE(11,*)
 DO i=1,24
  nrfy=2**i
  IF(nrfy>=nels)EXIT
 ENDDO
 yl=h1/DBLE(nels)*DBLE(nrfy)
 WRITE(11,'(A)')"Length of random field"
 WRITE(11,'(f8.2)')yl
 WRITE(11,*)
 WRITE(11,'(A)')"Number of random cells in column"
 WRITE(11,'(I5)')nrfy
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)nrv
 WRITE(11,'(A)')"Number of random variables"
 WRITE(11,'(I5)')nrv
 WRITE(11,*)
 maxfld=3*nrfy/2
!-------------------------------------------------------------------------
 ALLOCATE(gfld(maxfld,nrv),a(7,nrv),r(nrv,nrv),ac(2,nrfy,nrv))
 READ(10,*)trash
 DO i=1,nrv
  READ(10,*)a(:,i)
 ENDDO
 WRITE(11,'(A)')"nrv statistical properties"
 WRITE(11,'(A)')"     mean          SD         dist         L           U           m           s"
 DO i=1,nrv
  WRITE(11,'(2F12.2,F10.1,4F12.2)')a(:,i)
 ENDDO
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)r
 IF(SUM(ABS(r))>DBLE(nrv))lxfld=.TRUE.
 WRITE(11,'(A)')"Correlation matrix"
 DO i=1,nrv
  WRITE(11,'(6f8.2)')r(i,:)
 END DO
 WRITE(11,*)
 DO i=1,nrv
  CALL propst(a(:,i))
 ENDDO
 READ(10,*)trash
 READ(10,*)thy
 WRITE(11,'(A)')"Spatial correlation length"
 WRITE(11,'(F20.8)')thy
 WRITE(11,*)
 READ(10,*)trash
 READ(10,*)kseed,nsim
 kseed = iseed(kseed)
 WRITE(11,'(A)')'Starting seed'
 WRITE(11,'(I5)')kseed
 WRITE(11,*)
 WRITE(11,'(A,I2,A)')'Number of simulations'
 WRITE(11,'(I5)')nsim
 WRITE(11,*)
 READ(10,*)trash
 READ(10,'(A)')varfnc
 WRITE(11,'(A)')"Covariance function name"
 WRITE(11,'(A)')varfnc
 WRITE(11,*)
 eh=yl/nrfy
!---------- adjust parameters for linear case ----------------------------
 DO i=1,nrv
   IF(a(3,i)>threept5)THEN
     DO iel=1,nrfy
       ac(1,iel,i)=a(1,i)*(eh*hlf+(iel-1)*eh)+a(2,i)
       ac(2,iel,i)=a(4,i)*ac(1,iel,i)
       ac(2,iel,i)=alog(onn+ac(2,iel,i)*ac(2,iel,i)/(ac(1,iel,i)*ac(1,iel,i)))
       ac(1,iel,i)=alog(ac(1,iel,i))-hlf*ac(2,iel,i)
       ac(2,iel,i)=sqrt(ac(2,iel,i))
     END DO
   END IF
 END DO
 ALLOCATE(trans(nsim),rot(nsim))
!---------- FINITE ELEMENT PART-------------------------------------------
 np_types=nels
 nn=nels+1
 ALLOCATE(g(ndof),num(nod),nf(nodof,nn),etype(nels),ell(nels),eld(ndof),  &
   km(ndof,ndof),mm(ndof,ndof),action(ndof),g_g(ndof,nels),               &
   prop(nprops,np_types),ktots(nels))
 ell=h1/(nels)
 nf=1
 nr=0
 CALL formnf(nf)
 neq=MAXVAL(nf)  
 ALLOCATE(kdiag(neq),loads(0:neq),loadsi(0:neq),dtots(neq))
!-----------------------loop the elements to find global array sizes------
 kdiag=0
 elements_1: DO iel=1,nels
   num=(/iel,iel+1/)
   CALL num_to_g(num,nf,g)
   g_g(:,iel)=g
   CALL fkdiag(kdiag,g)
 END DO elements_1
 DO i=2,neq
   kdiag(i)=kdiag(i)+kdiag(i-1)
 END DO
!-------------Read Loads--------------------------------------------------
 loadsi=zero
 READ(10,*)trash
 READ(10,*)loaded_nodes
 ALLOCATE(kv(kdiag(neq)),no(loaded_nodes),stord(nsim,loaded_nodes),       &
   mn(loaded_nodes),sd(loaded_nodes))
 READ(10,*)(no(i),loadsi(nf(:,no(i))),i=1,loaded_nodes)
 WRITE(11,'(A,I2,A)')'  Node    Force        Moment'
 DO i=1,loaded_nodes
   WRITE(11,'(I5,2F12.4)')no(i),loadsi(nf(:,no(i)))
 END DO
!---------- Start Monte-Carlo simulations---------------------------------
 ktots=zero
 dtots=zero
 DO  jr = 1, nsim
  IF(jr/=1)nseed = iseed(kseed + jr - 1)
  CALL  sim1dm1(iterm,nrv,yl,nrfy,thy,varfnc,kseed,lxfld,                 &
                a,ac,R,gfld,nrv,maxfld,lsquare,debug)
  ktots=ktots+gfld(1:nels,2)
  IF(jr==1)THEN
    WRITE(11,'(/A)')"****k values at first simulation****"
    DO i=1,nels
      WRITE(11,'(I5,2F20.4)')i,gfld(i,2)
    END DO
    WRITE(11,*)
    WRITE(11,'(/A)')"****tip displacements at each simulation****"
  END IF
!-----------------------global stiffness matrix assembly------------------
 loads=loadsi
 kv=zero   
 elements_2: DO iel=1, nels
   prop(1,iel)=gfld(iel,1)  
   CALL beam_km(km,prop(1,iel),ell(iel))
   mm=zero
   prop(2,iel)=gfld(iel,2)    
   IF(nprops>1)CALL beam_mm(mm,prop(2,iel),ell(iel))
   g=g_g(:,iel)
   CALL fsparv(kv,km+mm,g,kdiag)
 END DO elements_2
!-----------------------equation solution --------------------------------
   CALL sparin(kv,kdiag)
   CALL spabac(kv,loads,kdiag)
   dtots=dtots+loads(1:201:2)
   loads(0)=zero
   WRITE(11,'(I5,12f20.8)')jr,loads(nf(1,no))
   stord(jr,:)=loads(nf(1,no))
 END DO
!-------------------Stop Monte-Carlo simulations--------------------------
 ktots=ktots/nsim
 dtots=dtots/nsim
 WRITE(11,'(/A)')"****mean k values after nsim simulations****"
 DO i=1,nels
   WRITE(11,'(I5,2F20.4)')i,ktots(i)
 END DO
 WRITE(11,'(/A)')"****mean deflections after nsim simulations****"
 DO i=1,nn
   WRITE(11,'(I5,2F20.8)')i,dtots(i)
 END DO
 DO i=1,loaded_nodes
   mn(i)=SUM(stord(:,i))/nsim
   sd(i)=SQRT(DOT_PRODUCT(stord(:,i)-mn(i),stord(:,i)-mn(i))/DBLE((nsim-1)))
 END DO
!
 WRITE(11,*)
 WRITE(11,'(/A)')"  ****************Tip Deflection Stats****************"
 WRITE(11,'(/A)')"  Node   Mean        SD"
 DO i=1,loaded_nodes
   WRITE(11,'(I5,2F12.8)')no(i),mn(i),sd(i)
 END DO
STOP
END PROGRAM rboef