post.f90

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


!==============================================================================!
!==============================================================================!
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!==============================================================================!
!==============================================================================!
!                                                                              |
!              POST PROCESSING    Apr. 2008                                    |
!                                                                              |
!==============================================================================!
!==============================================================================!

program data_postprocessing

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
integer :: iunit,ilsf,nne,is,k,npcl,j,np,nl
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_surface_icon
integer :: output_strain_tensor
integer :: output_ps
integer :: output_cubes
integer :: output_rivers
integer :: output_regular
integer iproc,nproc,ierr,nnn,nnne,nnnep,levelmax,ix,iy,iz,nlinks,iordermin
integer iconmin,iconmax
integer,dimension(:),allocatable::nstrain,donor,order
integer,dimension(:),allocatable::levs,li
integer nz,ni,ijk
integer,dimension(:),allocatable::invoid,elvoid,rtf,ftr
integer myicon(100),nstep,ndir,iter,ii,lsf,levmax


logical,dimension(:),allocatable::influid,do_it,subset_leaves,instrain

character clsf*3,c4*4,cc4*4,dir*128
character cs

double precision :: eps,dil,current_time,activation_time,zmin,xx,yy,zz,maxe2d,dist
double precision, allocatable, dimension (:) :: ov_nodee2d, ov_nodecrit
double precision, allocatable, dimension (:) :: countnode 
double precision,dimension(:,:),allocatable::s,n1,n2,n3
double precision,dimension(:,:,:),allocatable::strain,strainn
double precision,dimension(:),allocatable::xi,yi,zi,ui,vi,wi,si,ei
double precision,dimension(:),allocatable::s11,s12,s13,s22,s23,s33,str11,str12,str13,str22,str23,str33
double precision,dimension(:),allocatable::azimuth1,azimuth3,dip1,dip3
double precision zmax,dz,l1,l2,l3,n11,n12,n13,n21,n22,n23,n31,n32,n33,con

double precision dxyz,x,y,z,xcut
integer icut 


!==============================================================================!
!==============================================================================!

call system('clear')

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 step >'
read*,nstep
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'

print*,'Is it a debug file ? (y or n)'
read (*,'(a)') cs
select case (cs)
case ('y')
   print*,'Which debug file: please enter grid iteration number'
   read*,iter
   write (cc4,'(i4)') iter
   if (iter.lt.1000) cc4(1:1)='0'
   if (iter.lt.100) cc4(1:2)='00'
   if (iter.lt.10) cc4(1:3)='000'
case ('n')
   print*,'it is not a debug file'
case default 
   print *,'error: please input y or n'
   stop
end select

!==============================================================================
!==============================================================================


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//'.bin'
   open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'_'//cc4//'.bin',status='old',form='unformatted')
else
   write(*,*) 'output file to be processed: ', '../'//dir(1:ndir)//'/time_'//c4//'.bin'
   open(unit=7,file='../'//dir(1:ndir)//'/time_'//c4//'.bin',status='old',form='unformatted')
end if
write(*,*) '--------------------------------------------------------------------------'
open(unit=77,file='input_of_outputs.txt',status='old')
read(77,*) output_u_field 
read(77,*) output_v_field 
read(77,*) output_w_field 
read(77,*) output_disp_field 
read(77,*) output_press_field
read(77,*) output_temp_field 
read(77,*) output_e2d_field 
read(77,*) output_crit_field 
read(77,*) output_strain_field 
read(77,*) output_lsf_field 
read(77,*) output_surface_icon
read(77,*) output_strain_tensor
read(77,*) output_ps
read(77,*) output_cubes
read(77,*) output_rivers
read(77,*) output_regular

write(*,*) 'output u      field       ->',(output_u_field==1)
write(*,*) 'output v      field       ->',(output_v_field==1)
write(*,*) 'output w      field       ->',(output_w_field==1)
write(*,*) 'output disp   field       ->',(output_disp_field==1)
write(*,*) 'output press  field       ->',(output_press_field==1)
write(*,*) 'output temp   field       ->',(output_temp_field==1)
write(*,*) 'output e2d    field       ->',(output_e2d_field==1)
write(*,*) 'output crit   field       ->',(output_crit_field==1)
write(*,*) 'output strain field       ->',(output_strain_field==1)
write(*,*) 'output lsf    field       ->',(output_lsf_field==1)
write(*,*) 'output surfaces icon      ->',(output_surface_icon==1)
write(*,*) 'output strain tensor      ->',(output_strain_tensor==1)
write(*,*) 'output principal stresses ->',(output_ps==1)
write(*,*) 'output elemental fields   ->',(output_cubes==1)
write(*,*) 'output rivers             ->',(output_rivers==1)
write(*,*) 'output regular            ->',(output_regular==1)
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(invoid(nn))
allocate(rtf(nn),ftr(nn),elvoid(nn))
allocate(surface(ov%nlsf))
allocate(ov%nodal_pressure(nn))
allocate(ov%icon(8,ov%nleaves))
allocate(ov%octree(ov%noctree))

!========================
!=====[nodal values]=====
!========================

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]=====
!======================

read(7) (ov%icon(1:8,ie),            &
         ov%crit(ie),                &
         ov%e2d(ie),                 &
         ov%whole_leaf_in_fluid(ie), & 
         ie=1,ov%nleaves)

!==============================
!=====[octree information]=====
!==============================

read(7) (ov%octree(i),i=1,ov%noctree)

!================================
!=====[bad face information]=====
!================================
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) (invoid(i),elvoid(i),ftr(i),rtf(i),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))

call find_volumes (ov%icon,ov%octree,ov%noctree,ov%nleaves,ov%lsf,ov%nlsf,ov%nnode,influid, &
                   ov%x,ov%y,ov%z,ov%u,ov%v,ov%w,s,n1,n2,n3,strain)

! rescaling of the values for plotting 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
   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))
  
if (cs=='y') then
   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), &
               i=1,surface(is)%nsurface)
else
   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)
end if

   read (7) (surface(is)%icon(1:3,i),i=1,surface(is)%nt)

end do

close(7)


!==============================================================================
!==============================================================================
!=====[C -> N]=================================================================
!==============================================================================
!==============================================================================

allocate(countnode(ov%nnode))
allocate(ov_nodee2d(ov%nnode))
allocate(ov_nodecrit(ov%nnode))

ov_nodee2d       = 0.d0
ov_nodecrit      = 0.d0
countnode        = 0.d0
 
do i=1,ov%nleaves
   do j=1,8
      k=ov%icon(j,i)
      ov_nodee2d(k)       = ov_nodee2d(k)      + ov%e2d(i)    
      ov_nodecrit(k)      = ov_nodecrit(k)     + ov%crit(i)    
      countnode(k)        = countnode(k)       + 1                     
   end do
end do

do i=1,ov%nnode
   if (countnode(i).gt.0.d0) then
      ov_nodee2d(i)      = ov_nodee2d(i)      / countnode(i)
      ov_nodecrit(i)     = ov_nodecrit(i)     / countnode(i)
   end if
end do

deallocate(countnode)

!==============================================================================
!==============================================================================

write(*,'(a,i10)') 'octree length              ',ov%noctree
write(*,'(a,i10)') 'nb of nodes                ',ov%nnode
write(*,'(a,i10)') 'nb of leaves               ',ov%nleaves
write(*,'(a,i10)') 'nb of bad faces            ',ner
write(*,'(a,i10)') 'nb of lsf                  ',ov%nlsf
write(*,*) '--------------------------------------------------------------------------'
write(*,'(a,2f30.20)') 'x                 range : ',minval(ov%x(1:nn)),      maxval(ov%x(1:nn))
write(*,'(a,2f30.20)') 'y                 range : ',minval(ov%y(1:nn)),      maxval(ov%y(1:nn))
write(*,'(a,2f30.20)') 'z                 range : ',minval(ov%z(1:nn)),      maxval(ov%z(1:nn))
write(*,'(a,2f30.20)') 'u                 range : ',minval(ov%u(1:nn)),      maxval(ov%u(1:nn))
write(*,'(a,2f30.20)') 'v                 range : ',minval(ov%v(1:nn)),      maxval(ov%v(1:nn))
write(*,'(a,2f30.20)') 'w                 range : ',minval(ov%w(1:nn)),      maxval(ov%w(1:nn))
write(*,'(a,2f30.20)') 'pressure (nodes)  range : ',minval(ov%nodal_pressure), maxval(ov%nodal_pressure)
write(*,'(a,2f30.20)') 'e2d               range : ',minval(ov%e2d),          maxval(ov%e2d)
write(*,'(a,2f30.20)') 'e2d (nodes)       range : ',minval(ov_nodee2d),      maxval(ov_nodee2d)
write(*,'(a,2f30.20)') 'crit              range : ',minval(ov%crit),         maxval(ov%crit)
write(*,'(a,2f30.20)') 'crit (nodes)      range : ',minval(ov_nodecrit),     maxval(ov_nodecrit)
write(*,'(a,2f30.20)') 'strain            range : ',minval(ov%strain(1:nn)), maxval(ov%strain(1:nn))
do is=1,ov%nlsf
write(*,*) '--------------------------------------------------------------------------'
write(*,'(a,i3)') 'surface',is
write(*,'(a,i7,a)') 'surface counts ',surface(is)%nsurface, ' points'
write(*,'(a,i7,a)') 'surface counts ',surface(is)%nt, ' triangles'
write(*,'(a,2f30.20)') 'surf%x range : ',minval(surface(is)%x),maxval(surface(is)%x)
write(*,'(a,2f30.20)') 'surf%y range : ',minval(surface(is)%y),maxval(surface(is)%y)
write(*,'(a,2f30.20)') 'surf%z range : ',minval(surface(is)%z),maxval(surface(is)%z)
end do

nnn=0
nnne=0
do i=1,nn
if (invoid(i).eq.0) nnn=nnn+1
enddo
do i=1,ov%nleaves
if (elvoid(i).eq.0) nnne=nnne+1
enddo

!==============================================================================
!=====[producing cubes.vtk]====================================================
!==============================================================================
if (output_cubes==1) then

   !=====[selection of leaves we are interested in]====
   allocate(subset_leaves(ov%nleaves))
   maxe2d=maxval(ov%e2d)
   subset_leaves=.false.
   do i=1,ov%nleaves
      xx=ov%x(ov%icon(1,i))
      yy=ov%y(ov%icon(1,i))
      zz=ov%z(ov%icon(1,i))
      if (.not.(xx.gt..5d0 .and. yy.gt. .5d0) .and.  ov%e2d(i)>0.06 *maxe2d) subset_leaves(i)=.true.
   end do
   !subset_leaves=ov%whole_leaf_in_fluid

   nl=count(subset_leaves)
   !np=8*nl
   np=8*ov%nleaves

   iunit=30
   open(unit=iunit,file='cubes.vtk')
   write(iunit,'(a)')'# vtk DataFile Version 3.0'
   write(iunit,'(a)')'velocities'
   write(iunit,'(a)')'ASCII'
   write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
   write(iunit,'(a7,i10,a6)')'POINTS ',np,' float'

   do i=1,ov%nleaves
      !if (subset_leaves(i)) then
         do k=1,8 
           xx=ov%x(ov%icon(k,i))
           yy=ov%y(ov%icon(k,i))
           zz=ov%z(ov%icon(k,i))
           write(iunit,'(3e11.4)') xx,yy,zz
         end do
      !end if
   enddo

   write(iunit,'(A6, 2I10)') 'CELLS ',nl,(1+8)*nl
   do ie=1,ov%nleaves
      if (subset_leaves(ie)) then
         write(iunit,'(9I10)') 8 , ((ie-1)*8+i-1,i=1,8)
      end if
   end do

   write(iunit,'(A11, I10)') 'CELL_TYPES ',nl
   do ie=1,ov%nleaves
      if (subset_leaves(ie)) then
         write(iunit,'(I2)') 11
      end if 
   end do

   write(iunit,'(a11,i10)')'POINT_DATA ',np

   write(iunit,'(a)')'SCALARS e2d float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
         do k=1,8 
            write(iunit,'(e11.4)') ov%e2d(i)
         end do
   end do
   
   write(iunit,'(a)')'SCALARS w float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)')sum(ov%w(ov%icon(:,i)))/8.d0
      end do
   end do
   write(iunit,'(a)')'SCALARS v float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)')sum(ov%v(ov%icon(:,i)))/8.d0
      end do
   end do
   write(iunit,'(a)')'SCALARS u float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)')sum(ov%u(ov%icon(:,i)))/8.d0
      end do
   end do

   write(iunit,'(a)')'SCALARS uvw float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,ov%nleaves
      do k=1,8 
      write(iunit,'(e11.4)') sqrt((sum(ov%u(ov%icon(:,i)))/8.d0)**2+(sum(ov%v(ov%icon(:,i)))/8.d0)**2+(sum(ov%w(ov%icon(:,i)))/8.d0)**2)
      end do
   end do

   write(*,*) '--------------------------------------------------------------------------'
   write(*,*) 'opla I am done producing cubes.vtk'
end if




!==============================================================================
!======[producing strain.vtk]==================================================
!==============================================================================

if (output_strain_tensor/=0) then

nnnep=0
levelmax=2
allocate (instrain(ov%nleaves))
instrain=.FALSE.
do i=1,ov%nleaves
   if (elvoid(i).eq.0) then
   xx=sum(ov%x(ov%icon(:,i)))/8.d0
   yy=sum(ov%y(ov%icon(:,i)))/8.d0
   zz=sum(ov%z(ov%icon(:,i)))/8.d0
   xx=xx*(2.d0**levelmax)
   yy=yy*(2.d0**levelmax)
   zz=zz*(2.d0**levelmax)
   ix=int(xx)
   iy=int(yy)
   iz=int(zz)
   if (abs(xx-float(ix)-0.5d0).lt.1.d-6 .and. &
       abs(yy-float(iy)-0.5d0).lt.1.d-6 .and. &
       abs(zz-float(iz)-0.5d0).lt.1.d-6) instrain(i)=.TRUE.
   if (instrain(i)) nnnep=nnnep+1
   endif
enddo 
!print*,'here',nnne,nnnep

iunit=30
open(unit=iunit,file='strain.vtk')
write(iunit,'(a)')'# vtk DataFile Version 3.0'
write(iunit,'(a)')'sigma1'
write(iunit,'(a)')'ASCII'
write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
write(iunit,'(a7,i10,a6)')'POINTS ',nnnep,' float'

do i=1,ov%nleaves
   if (elvoid(i).eq.0) then
   xx=sum(ov%x(ov%icon(:,i)))/8.d0
   yy=sum(ov%y(ov%icon(:,i)))/8.d0
   zz=sum(ov%z(ov%icon(:,i)))/8.d0
   if (instrain(i)) write(iunit,'(3f16.11)') xx,yy,zz
   endif
enddo 

write(iunit,'(a11,i10)')'POINT_DATA ',nnnep

write(iunit,'(a)')'TENSORS strain float'
do i=1,ov%nleaves
   if (elvoid(i).eq.0) then
   xx=sum(ov%x(ov%icon(:,i)))/8.d0
   yy=sum(ov%y(ov%icon(:,i)))/8.d0
   zz=sum(ov%z(ov%icon(:,i)))/8.d0
   if (instrain(i)) then
   write(iunit,'(3e11.4)') strain(1,:,i)
   write(iunit,'(3e11.4)') strain(2,:,i)
   write(iunit,'(3e11.4)') strain(3,:,i)
   endif
   endif
enddo
close (iunit)

deallocate (instrain)

write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing strain.vtk '

end if

!==============================================================================
!======[produce compressive.vtk file]==========================================
!==============================================================================

if (output_ps/=0) then

allocate(do_it(ov%nleaves))

np=0
do i=1,ov%nleaves
   zz=sum(ov%z(ov%icon(:,i)))/8.d0  !criterion
   if (zz<.3d0) then
      np=np+1
      do_it(i)=.true.
   end if
end do

write(*,*) '--------------------------------------------------------------------------'
write(*,*) 'subset of points in compressive.vtk:',np

iunit=30
open(unit=iunit,file='compressive.vtk')
write(iunit,'(a)')'# vtk DataFile Version 3.0'
write(iunit,'(a)')'sigma1'
write(iunit,'(a)')'ASCII'
write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
write(iunit,'(a7,i10,a6)')'POINTS ',np,' float'

do i=1,ov%nleaves
   if (do_it(i)) then
      xx=sum(ov%x(ov%icon(:,i)))/8.d0
      yy=sum(ov%y(ov%icon(:,i)))/8.d0
      zz=sum(ov%z(ov%icon(:,i)))/8.d0
      write(iunit,'(3f16.11)') xx,yy,zz
   end if
enddo

write(iunit,'(a11,i10)')'POINT_DATA ', np

write(iunit,'(a)')'VECTORS sigma1head+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') n1(1:3,i)*s(1,i)
enddo
write(iunit,'(a)')'VECTORS sigma1head- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') -n1(1:3,i)*s(1,i)
enddo

write(iunit,'(a)')'VECTORS sigma3tail+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') n2(1:3,i)*s(2,i)
enddo

write(iunit,'(a)')'VECTORS sigma3tail- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') -n2(1:3,i)*s(2,i)
enddo

write(iunit,'(a)')'VECTORS sigma2head+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') n3(1:3,i)*min(0.d0,s(3,i))
enddo

write(iunit,'(a)')'VECTORS sigma2head- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') -n3(1:3,i)*min(0.d0,s(3,i))
enddo

write(iunit,'(a)')'VECTORS sigma2tail+ float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') n3(1:3,i)*max(0.d0,s(3,i))
enddo

write(iunit,'(a)')'VECTORS sigma2tail- float'
do i=1,ov%nleaves
   if (do_it(i)) & 
   write(iunit,'(3e11.4)') -n3(1:3,i)*max(0.d0,s(3,i))
enddo

close (iunit)

deallocate(do_it)

write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing compressive.vtk '

end if

!==============================================================================
!======[produce visu.vtk file]=================================================
!==============================================================================


!iunit=30
!open(unit=iunit,file='cut.dat')
!do i=1,nnn              ! loop over nodes
!   x=ov%x(rtf(i))
!   y=ov%y(rtf(i))
!   z=ov%z(rtf(i))
!   icut=3   ! icut between 0 and 64 for uniform level 6
!   ! example , cross section in plane given by x=value
!   xcut=icut*2.d0**(-6)    ! background uniform grid level is 6
!   dxyz=2.d0**(-7)      ! higher level grid in boxes is 7
!   if (abs(x-xcut)<1.d-6 ) then   ! select nodes that are only on the plane
!      if (mod(nint(y/dxyz),2)==0) then  ! les coordonnees des noeuds surla grille de niveau 7 divisees par dxy vont prendre les valeurs 0 a 128
!      if (mod(nint(z/dxyz),2)==0) then
!         write(iunit,*) x,y,z
!      end if
!      end if
!   end if
!enddo
!close(iunit)


iunit=30
open(unit=iunit,file='visu.vtk')
write(iunit,'(a)')'# vtk DataFile Version 3.0'
write(iunit,'(a)')'velocities'
write(iunit,'(a)')'ASCII'
write(iunit,'(a)')'DATASET UNSTRUCTURED_GRID'
write(iunit,'(a7,i10,a6)')'POINTS ',nnn,' float'

do i=1,nnn
write(iunit,'(3f16.11)')ov%x(rtf(i)),ov%y(rtf(i)),ov%z(rtf(i))
enddo

write(iunit,'(A6, 2I10)') 'CELLS ',nnne,(1+nne)*nnne
iconmin=nn
iconmax=0
do ie=1,ov%nleaves
   if (elvoid(ie).eq.0) then
   myicon(1:nne)=ftr(ov%icon(1:nne,ie))-1
     do k=1,nne
     iconmin=min(iconmin,myicon(k))
     iconmax=max(iconmax,myicon(k))
     enddo
   write(iunit,'(9I10)')nne,myicon(1:nne)
   endif
enddo
!print*,iconmin,iconmax

write(iunit,'(A11, I10)') 'CELL_TYPES ',nnne
do ie=1,nnne
   write(iunit,'(I2)')11 ! octree  (8 nodes)
end do

write(iunit,'(a11,i10)')'POINT_DATA ',nnn

if(output_u_field==1) then 
   write(iunit,'(a)')'SCALARS u float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%u(rtf(i))
   enddo
end if

if(output_v_field==1) then 
   write(iunit,'(a)')'SCALARS v float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%v(rtf(i))
   enddo
end if

if(output_w_field==1) then 
   write(iunit,'(a)')'SCALARS w float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%w(rtf(i))
   enddo
end if


if(output_disp_field==1) then 
   write(iunit,'(a)')'SCALARS disp float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') sqrt(ov%u(rtf(i))**2+ov%v(rtf(i))**2+ov%w(rtf(i))**2)
   enddo
end if

if(output_press_field==1) then 
   write(iunit,'(a)')'SCALARS p float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%nodal_pressure(rtf(i))
   enddo
end if


if(output_e2d_field==1) then 
   write(iunit,'(a)')'SCALARS e2d float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov_nodee2d(rtf(i))
   enddo
end if

if(output_crit_field==1) then 
   write(iunit,'(a)')'SCALARS crit float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov_nodecrit(rtf(i))
   enddo
end if


if(output_strain_field==1) then 
   write(iunit,'(a)')'SCALARS s float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%strain(rtf(i))
   enddo
end if

if(output_temp_field==1) then 
   write(iunit,'(a)')'SCALARS t float 1'
   write(iunit,'(a)')'LOOKUP_TABLE default'
   do i=1,nnn
      write(iunit,'(e11.4)') ov%temp(rtf(i))
   enddo
end if

if(output_lsf_field==1) then 
   do ilsf=1,ov%nlsf
      write (clsf,'(i2)') ilsf
      if (ilsf.lt.10) clsf(1:1)='0'
      write(iunit,'(a)')'SCALARS lsf'//clsf//' float 1'
      write(iunit,'(a)')'LOOKUP_TABLE default'
      do i=1,nnn
         write(iunit,'(e11.4)') ov%lsf(rtf(i),ilsf)
      enddo
   enddo
end if
write(iunit,'(a)')'VECTORS velo float'
do i=1,nnn
   write(iunit,'(3e11.4)')ov%u(rtf(i)),ov%v(rtf(i)),ov%w(rtf(i))
enddo

if (output_strain_tensor/=0) then
allocate (strainn(3,3,nnn),nstrain(nnn))
strainn=0.d0
nstrain=0
  do ie=1,nnne
    if (elvoid(ie).eq.0) then
    do k=1,nne
    strainn(:,:,ftr(ov%icon(k,ie)))=strainn(:,:,ftr(ov%icon(k,ie)))+strain(:,:,ie)
    nstrain(ftr(ov%icon(k,ie)))=nstrain(ftr(ov%icon(k,ie)))+1
    enddo
    endif
  enddo
do i=1,nnn
if (nstrain(i).ne.0) strainn(:,:,i)=strainn(:,:,i)/nstrain(i)
enddo
write(iunit,'(a)')'TENSORS strain float'
do i=1,nnn
   write(iunit,'(3e11.4)') strainn(1,:,i)
   write(iunit,'(3e11.4)') strainn(2,:,i)
   write(iunit,'(3e11.4)') strainn(3,:,i)
enddo
deallocate (strainn,nstrain)
endif

write(iunit,'(a)')'done visu.vtk'
close(iunit)
write(*,*) '--------------------------------------------------------------------------'
write(*,*)'opla I am done producing visu.vtk '



!==============================================================================
!======[produce norm<i>.vtk file]==============================================
!==============================================================================

if (output_surface_icon==1) then 
   do is=1,ov%nlsf
      write(cs,'(i1)') is
      open(unit=30,file='norm'//cs//'.vtk')
      write(30,'(a)')'# vtk DataFile Version 3.0'
      write(30,'(a)')'surface'
      write(30,'(a)')'ASCII'
      write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
      write(30,'(a7,i10,a6)')'POINTS ',surface(is)%nsurface,' float'
      do i=1,surface(is)%nsurface
         write(30,'(3f16.11)')surface(is)%x(i),surface(is)%y(i),surface(is)%z(i)
      enddo
      write(30,'(a6, 2I10)') 'CELLS ',surface(is)%nt,(1+3)*surface(is)%nt
      do i=1,surface(is)%nt
         write(30,'(9I10)')3,surface(is)%icon(1:3,i)-1
      enddo
      write(30,'(A11, I10)') 'CELL_TYPES ',surface(is)%nt
      do i=1,surface(is)%nt
         write(30,'(I2)')5 ! octree  (8 nodes)
      end do
      write(30,'(a11,i10)') 'POINT_DATA ',surface(is)%nsurface
      write(30,'(a)')'SCALARS z float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
      do i=1,surface(is)%nsurface
         write(30,'(e11.4)') surface(is)%z(i)
      enddo
if (is.eq.1) then
      write(30,'(a)')'SCALARS erosion_rate float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
      do i=1,surface(is)%nsurface
         write(30,'(e11.4)') -surface(is)%r(i)
      enddo
      write(30,'(a)')'SCALARS rock_uplift float 1'
      write(30,'(a)')'LOOKUP_TABLE default'
      do i=1,surface(is)%nsurface
         write(30,'(e11.4)') surface(is)%w(i)
      enddo
endif
      write(30,'(a)') 'VECTORS normal float'
      do i=1,surface(is)%nsurface
         write(30,'(3e15.4)') surface(is)%xn(i),surface(is)%yn(i),surface(is)%zn(i)
      end do
      write(30,'(a)') 'VECTORS vel float'
      do i=1,surface(is)%nsurface
         write(30,'(3e15.4)') surface(is)%u(i),surface(is)%v(i),surface(is)%w(i)
      end do
      close(30)
      write(*,*) '--------------------------------------------------------------------------'
      write(*,*)'opla I am done producing norm'//cs//'.vtk' 
   enddo
end if

!==============================================================================
!======[produce rivers.vtk file]==============================================
!==============================================================================

if (output_rivers==1) then
allocate (donor(surface(1)%nsurface),order(surface(1)%nsurface))
call rivers (surface(1)%icon,surface(1)%nt,surface(1)%x,surface(1)%y,surface(1)%z, &
             surface(1)%nsurface,donor,order)
iordermin=3
nlinks=0
do i=1,surface(1)%nsurface
   if (order(i).gt.iordermin.and.donor(i).gt.0) nlinks=nlinks+1
enddo
      open(unit=30,file='rivers.vtk')
      write(30,'(a)')'# vtk DataFile Version 3.0'
      write(30,'(a)')'rivers'
      write(30,'(a)')'ASCII'
      write(30,'(a)')'DATASET UNSTRUCTURED_GRID'
      write(30,'(a7,i10,a6)')'POINTS ',surface(1)%nsurface,' float'
      do i=1,surface(1)%nsurface
         write(30,'(3f16.11)')surface(1)%x(i),surface(1)%y(i),surface(1)%z(i)+.001
      enddo
      write(30,'(a6, 2I10)') 'CELLS ',nlinks,(1+2)*nlinks
      do i=1,surface(1)%nsurface
         if (order(i).gt.iordermin.and.donor(i).gt.0) write(30,'(9I10)')2,i-1,donor(i)-1
      enddo
      write(30,'(A11, I10)') 'CELL_TYPES ',nlinks
      do i=1,nlinks
         write(30,'(I2)')3 ! octree  (2 nodes)
      end do
      close(30)
      write(*,*) '--------------------------------------------------------------------------'
      write(*,*)'opla I am done producing rivers.vtk' 
deallocate (donor,order)
end if

!==============================================================================
!======[produce regular.vtk file]==============================================
!==============================================================================

if (output_regular==1) then

allocate (levs(ov%nleaves))
call octree_find_element_level (ov%octree,ov%noctree,levs,ov%nleaves)
levmax=maxval(levs)

zmax=maxval(surface(1)%z)
dz=1.d0/(2.d0**levmax)
nz=(zmax/dz)+1

allocate (str11(ov%nnode),str12(ov%nnode),str13(ov%nnode), &
          str22(ov%nnode),str23(ov%nnode),str33(ov%nnode),nstrain(ov%nnode))

nstrain=0
str11=0.;str12=0.;str13=0.
str22=0.;str23=0.;str33=0.
  do ie=1,ov%nleaves
    if (elvoid(ie).eq.0) then
    do k=1,nne
    i=ov%icon(k,ie)
    str11(i)=str11(i)+strain(1,1,ie)
    str12(i)=str12(i)+strain(1,2,ie)
    str13(i)=str13(i)+strain(1,3,ie)
    str22(i)=str22(i)+strain(2,2,ie)
    str23(i)=str23(i)+strain(2,3,ie)
    str33(i)=str33(i)+strain(3,3,ie)
    nstrain(i)=nstrain(i)+1
    enddo
    endif
  enddo
do i=1,ov%nnode
  if (nstrain(i).ne.0) then
  str11(i)=str11(i)/nstrain(i)
  str12(i)=str12(i)/nstrain(i)
  str13(i)=str13(i)/nstrain(i)
  str22(i)=str22(i)/nstrain(i)
  str23(i)=str23(i)/nstrain(i)
  str33(i)=str33(i)/nstrain(i)
  endif
enddo

!print*,minval(str11),maxval(str11)
!print*,minval(str12),maxval(str12)
!print*,minval(str13),maxval(str13)
!print*,minval(str22),maxval(str22)
!print*,minval(str23),maxval(str23)
!print*,minval(str33),maxval(str33)

ni=nz*(2**levmax+1)*(2**levmax+1)
allocate (xi(ni),yi(ni),zi(ni),ui(ni),vi(ni),wi(ni),si(ni),ei(ni),li(ni), &
          s11(ni),s12(ni),s13(ni),s22(ni),s23(ni),s33(ni))
do k=1,nz
do j=1,2**levmax+1
do i=1,2**levmax+1
ijk=ijk+1
xi(ijk)=(i-1)*dz
yi(ijk)=(j-1)*dz
zi(ijk)=(k-1)*dz
call octree_interpolate_many (12,ov%octree,ov%noctree,ov%icon,ov%nleaves,ov%nnode, &
                              xi(ijk),yi(ijk),zi(ijk), &
                              ov%u,ui(ijk),ov%v,vi(ijk),ov%w,wi(ijk), &