Skip to content
Snippets Groups Projects
surface.f90 10.2 KiB
Newer Older
  • Learn to ignore specific revisions
  • 
    MODULE definitions
    
       type structure
          integer,dimension(:),pointer::octree
          integer noctree,nnode,nleaves
          integer nelemr,nlsf 
          double precision,dimension(:),pointer::x,y,z
          double precision,dimension(:),pointer::u,v,w,temp,pressure,strain
          double precision,dimension(:),pointer::nodal_pressure
          double precision,dimension(:),pointer::e2d
          double precision,dimension(:),pointer::crit
          integer,dimension(:,:),pointer::icon,iconr
          double precision,dimension(:,:),pointer:: lsf
          integer,dimension(:),pointer::on(:)
          logical, dimension(:), pointer :: whole_leaf_in_fluid
       end type structure
    
       type sheet
          integer nsurface,nt
          double precision,dimension(:),pointer::r,s,x,y,z,xn,yn,zn,u,v,w
          integer,dimension(:,:),pointer::icon
          integer nhull,levelt,itype,material
       end type sheet
    
    END MODULE definitions
    
    
    !------------------------------------------------------------------------------|
    !------------------------------------------------------------------------------|
    !                                                                              |
    !              ||===\\                                                         | 
    !              ||    \\                                                        |
    !              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
    !              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
    !              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
    !              ||===//     \\==//    \\==\\  \\==\\  ||                        |
    !                                                                              |
    !------------------------------------------------------------------------------|
    !------------------------------------------------------------------------------|
    !                                                                              |
    !              ROUTINE    Nov. 2006                                            |
    !                                                                              |
    !------------------------------------------------------------------------------|
    !------------------------------------------------------------------------------|
    
    program data_postprocessing
    
    !------------------------------------------------------------------------------|
    !(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
    !------------------------------------------------------------------------------|
    
    !------------------------------------------------------------------------------|
    !(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
    !------------------------------------------------------------------------------|
    use definitions
    
    implicit none
    
    type (sheet),dimension(:),allocatable::surface
    type(structure),SAVE :: ov
    integer :: ie,mer,nn,ner,i,i1,i2,i3,i4,i5,kfx,kfy,kfz,kft,iunit,ilsf,nne,is,k,npcl,j,np
    integer :: output_u_field
    integer :: output_v_field
    integer :: output_w_field
    integer :: output_disp_field
    integer :: output_press_field
    integer :: output_temp_field
    integer :: output_e2d_field
    integer :: output_crit_field
    integer :: output_strain_field
    integer :: output_lsf_field
    integer :: output_lsf_ptcles
    integer :: output_surface_icon
    integer :: output_strain_tensor
    integer :: output_ps
    integer :: output_cubes
    double precision :: eps,dil,current_time,activation_time,zmin,xx,yy,zz
    logical,dimension(:),allocatable::influid,do_it
    integer myicon(100),nstep,ndir,iter,ii,lsf,nstep0,nstep1
    character clsf*3,c4*4,cc4*4,dir*128
    character cs
    double precision veg,hhh,lll,ztop
    
    double precision, allocatable, dimension (:) :: ov_nodepressure, ov_nodee2d, ov_nodecrit
    double precision, allocatable, dimension (:) :: countnode 
    double precision,dimension(:,:),allocatable::s,n1,n2,n3
    double precision,dimension(:,:,:),allocatable::strain
    
    integer iproc,nproc,ierr,nstepi
    
    !==============================================================================
    !==============================================================================
    
    nne=8
    
    print*,'Enter directory (default is OUT)'
    read (*,'(a)') dir
    ndir=len_trim(dir)
    if (ndir.eq.0) then
      ndir=3
      dir(1:ndir)='OUT'
    endif
    
    print*,'Enter time steps >'
    read*,nstep0,nstep1,nstepi
    
    print*,'Enter top surface position'
    read*,ztop
    
    open (57,file='surface.dat',status='unknown')
    open (47,file='velocity-x.dat',status='unknown')
    open (48,file='velocity-z.dat',status='unknown')
    open (49,file='ball.dat',status='unknown')
    
    do nstep=nstep0,nstep1,nstepi
    
    write (c4,'(i4)') nstep
    if (nstep.lt.1000) c4(1:1)='0'
    if (nstep.lt.100) c4(1:2)='00'
    if (nstep.lt.10) c4(1:3)='000'
    
    !==============================================================================
    !==============================================================================
    
    
    write(*,*) '**************************************************************************'
    write(*,*) '--------------------------------------------------------------------------'
    write(*,*) '--------------------------------------------------------------------------'
    write(*,*) '------------------- VTK file production in progress ----------------------'
    write(*,*) '--------------------------------------------------------------------------'
    write(*,*) '--------------------------------------------------------------------------'
    if (cs=='y') then
       write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.txt'
       open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.txt',status='old')
    else
       write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'.txt'
       open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'.txt',status='old')
    end if
    write(*,*) '--------------------------------------------------------------------------'
    
    read (7,*) ov%noctree,ov%nnode,ov%nleaves,ner,ov%nlsf,npcl,current_time
    nn=ov%nnode
    allocate(ov%on(nn))
    allocate(ov%x(nn),ov%y(nn),ov%z(nn))
    allocate(ov%u(nn),ov%v(nn),ov%w(nn))
    allocate(ov%pressure(ov%nleaves),ov%strain(nn))
    allocate(ov%e2d(ov%nleaves))
    allocate(ov%whole_leaf_in_fluid(ov%nleaves))
    allocate(ov%crit(ov%nleaves))
    allocate(ov%temp(nn))
    allocate(ov%lsf(nn,ov%nlsf))
    allocate(influid(nn))
    allocate(surface(ov%nlsf))
    allocate(ov%nodal_pressure(nn))
    
    read(7,*)(ov%x(i),ov%y(i),ov%z(i),   &
              ov%u(i),ov%v(i),ov%w(i),   &
              ov%lsf(i,1:ov%nlsf),       &
              ov%temp(i),                &
              ov%nodal_pressure(i),       &
              ov%strain(i),              &
              kfx,kfy,kfz,kft,i=1,nn)
    
    
    ! icon array
    allocate(ov%icon(8,ov%nleaves))
    read(7,*) (ov%icon(1:8,ie),            &
    !           ov%pressure(ie),            &
               ov%crit(ie),                &
               ov%e2d(ie),                 &
               ov%whole_leaf_in_fluid(ie), &
               ie=1,ov%nleaves)
    
    ! octree information
    allocate(ov%octree(ov%noctree))
    read(7,*) (ov%octree(i),i=1,ov%noctree)
    
    ! bad face infromation
    mer=ner   ! mer=1
    if(ner.gt.mer)mer=ner
    allocate(ov%iconr(9,mer))
    read(7,*) (ov%iconr(1:9,ie),ie=1,ner)
    
    ov%nelemr=ner
    
    ! void information
    read (7,*) (i1,i2,i3,i4,influid(i),i=1,nn)
    ov%on=0
    do i=1,nn
       if (influid(i)) ov%on(i)=1
    enddo
    
    read(7,*) (i5,i=1,ner)
    
    !==============================================================================
    !=====[measure work, strain, compute principal directions]====================
    !==============================================================================
    
    allocate (s (3,ov%nleaves))
    allocate (n1(3,ov%nleaves))
    allocate (n2(3,ov%nleaves))
    allocate (n3(3,ov%nleaves))
    allocate (strain(3,3,ov%nleaves))
    
    ! rescaling of the values for ploiting purposes
    
    s=s/maxval(abs(s))*.01
    
    strain=strain/maxval(abs(strain))*.01
    
    !==============================================================================
    !=====[surface information]====================================================
    !==============================================================================
    
    do is=1,ov%nlsf
       read(7,*) surface(is)%nsurface,activation_time,surface(is)%nt
    print*,surface(is)%nsurface
       allocate (surface(is)%x(surface(is)%nsurface), &
                    surface(is)%y(surface(is)%nsurface), &
                    surface(is)%z(surface(is)%nsurface), &
                    surface(is)%xn(surface(is)%nsurface), &
                    surface(is)%yn(surface(is)%nsurface), &
                    surface(is)%zn(surface(is)%nsurface), &
                    surface(is)%r(surface(is)%nsurface), &
                    surface(is)%s(surface(is)%nsurface), &
                    surface(is)%u(surface(is)%nsurface), &
                    surface(is)%v(surface(is)%nsurface), &
                    surface(is)%w(surface(is)%nsurface), &
                    surface(is)%icon(3,surface(is)%nt))
      
       read (7,*) (surface(is)%r(i),surface(is)%s(i), &
                   surface(is)%x(i),surface(is)%y(i),surface(is)%z(i), &
                   surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i), &
                   surface(is)%u(i),surface(is)%v(i),surface(is)%w(i), &
                   i=1,surface(is)%nsurface)
    
       read (7,*) (surface(is)%icon(1:3,i),i=1,surface(is)%nt)
    
    if (is.eq.1) then
    write (57,'(a1,1x,2e15.8)') '&',current_time
    write (47,'(a)') '&'
    write (48,'(a)') '&'
    do i=1,surface(is)%nsurface
    if ((surface(is)%y(i)-.4999)*(surface(is)%y(i)-.5001).le.0.) then
    write (57,'(f14.8,2x,3e16.8)') surface(is)%x(i),surface(is)%z(i)-ztop,surface(is)%u(i),surface(is)%v(i)
    write (47,'(f14.8,2x,e15.8)') surface(is)%x(i),surface(is)%u(i)
    write (48,'(f14.8,2x,e15.8)') surface(is)%x(i),surface(is)%w(i)
    endif
    enddo
    endif
    
    
    if (is.eq.2) then
    write (49,'(3e15.8)') current_time,sum(surface(2)%z)/surface(2)%nsurface,sum(surface(2)%w)/surface(2)%nsurface
    endif
    
       deallocate (surface(is)%x, &
                    surface(is)%y, &
                    surface(is)%z, &
                    surface(is)%xn, &
                    surface(is)%yn, &
                    surface(is)%zn, &
                    surface(is)%r, &
                    surface(is)%s, &
                    surface(is)%icon)
      
    end do
    
    close(7)
    
    deallocate(ov%on)
    deallocate(ov%x,ov%y,ov%z)
    deallocate(ov%u,ov%v,ov%w)
    deallocate(ov%pressure,ov%strain)
    deallocate(ov%e2d)
    deallocate(ov%crit)
    deallocate(ov%temp)
    deallocate(ov%lsf)
    deallocate(influid)
    deallocate(surface)
    deallocate(ov%icon)
    deallocate(ov%octree)
    deallocate(ov%iconr)
    deallocate (s)
    deallocate (n1)
    deallocate (n2)
    deallocate (n3)
    deallocate (strain)
    
    enddo
    
    close (57)
    
    end