plot_postscript.f90

Spectral Element Method for wave propagation and rupture dynamics.

        call plot_boundaries()
      endif

   !-- plot a nodal scalar field as colored GLL cells
    else
      call plot_scal(vfield(:,1))
      if (mesh) call plot_mesh(color,numbers, .false.)
      call plot_boundaries()
    endif

 !-- plot an element scalar field as colored elements
  else if (elem_field) then
    call plot_efield()
    if (mesh) call plot_mesh(color,numbers, .false. )
    call plot_boundaries()
  endif


  if (symbols) then 
    if(color) then ! Sources and receivers in color ?
      write(psunit,*) 'Colreceiv'
    else
      write(psunit,*) '0 setgray'
    endif
    if (present(src) .and. associated(src)) call plot_sources()
    if (present(rec) .and. associated(rec)) call plot_receivers()
  endif

  write(psunit,*) '%'
  write(psunit,*) 'grestore'
  write(psunit,*) 'showpage'

  close(psunit)
  if (echo_run) write(iout,fmtok)

  return
 
  contains 


!-----------------------------------------------------------------------
! Plot headers

  subroutine plot_header()

  write(psunit,10) stitle,version
  write(psunit,*) '/CM {28.5 mul} def'  ! convert from cm to points (1/72 inch)
  write(psunit,*) '/L {lineto} def'
  write(psunit,*) '/LR {rlineto} def'
  write(psunit,*) '/M {moveto} def'
  write(psunit,*) '/MR {rmoveto} def'
  write(psunit,*) '/MK {mark} def' ! start array construction
  write(psunit,*) '/ST {stroke} def'
  write(psunit,*) '/CP {closepath} def'
  write(psunit,*) '/RG {setrgbcolor} def'
  write(psunit,*) '/GF {gsave fill grestore} def'
  write(psunit,*) '/GG {0 setgray ST} def'
  write(psunit,*) '/GC {Colmesh ST} def'
  write(psunit,*) '/RF {setrgbcolor fill} def'
  write(psunit,*) '/SF {setgray fill} def'
  write(psunit,*) '/GS {gsave} def'
  write(psunit,*) '/GR {grestore} def'
  write(psunit,*) '/SLW {setlinewidth} def'
  write(psunit,*) '/SCSF {scalefont setfont} def'
  write(psunit,*) '%---- symbols '
  write(psunit,*) '/Point {2 0 360 arc CP 0 setgray fill} def'
  write(psunit,*) '/VDot {-0.75 -1.5 MR 1.5 0 LR 0 3. LR -1.5 0 LR CP fill} def'
  write(psunit,*) '/HDot {-1.5 -0.75 MR 3. 0 LR 0 1.5 LR -3. 0 LR CP fill} def'
  write(psunit,*) '/Cross {GS 0.05 CM SLW GS 3 3 MR -6. -6. LR ST GR'
  write(psunit,*) 'GS 3 -3 MR -6. 6. LR ST GR 0.01 CM SLW} def'
  write(psunit,*) '/SmallLine {M 0.07 CM 0 rlineto} def'
  write(psunit,*) '/Losange {GS 0.05 CM SLW 0 4.2 MR -3 -4.2 LR 3 -4.2 LR'
  write(psunit,*) '3 4.2 LR CP ST GR 0.01 CM SLW} def'
  write(psunit,*) '%---- color settings'
  write(psunit,*) '% vector fields in magenta'
  write(psunit,*) '/Colvects {0.01 CM SLW 1. 0. 1. RG} def'
  write(psunit,*) '% element mesh in chartre'
  write(psunit,*) '/Colmesh {0.02 CM SLW 0.5 1. 0. RG} def'
  write(psunit,*) '% source and receivers in cyan'
  write(psunit,*) '/Colreceiv {0. 1. 1. RG} def'
  write(psunit,*) '%---- macros'
  write(psunit,*) '% arrow'
  write(psunit,*) '/F {M LR gsave LR ST grestore LR ST} def'
  write(psunit,*) '% element contour'
  write(psunit,*) '/CO {M counttomark 2 idiv {L} repeat cleartomark CP} def'
  write(psunit,*) '% filled quad cell'
  write(psunit,*) '/FQ {M L L L CP SF} def'
  write(psunit,*) '/FQC {M L L L CP RF} def'
  write(psunit,*) '%'
  write(psunit,*) '.01 CM SLW'
  write(psunit,*) '/Times-Roman findfont'
  write(psunit,*) '.35 CM SCSF'
  write(psunit,*) '%'
  write(psunit,*) '/vshift ',-height/2,' CM def'
  write(psunit,*) '/Rshow { currentpoint stroke M'
  write(psunit,*) 'dup stringwidth pop neg vshift MR show } def'
  write(psunit,*) '/Cshow { currentpoint stroke M'
  write(psunit,*) 'dup stringwidth pop -2 div vshift MR show } def'
  write(psunit,*) '/fN {/Helvetica-Bold findfont ',height,' CM SCSF} def'
  write(psunit,*) '%'
  write(psunit,*) 'gsave newpath 90 rotate'
  write(psunit,*) '0 ',-sizez,' CM translate 1. 1. scale'
  write(psunit,*) '%'

  return

10   format('%!PS-Adobe-2.0',/,'%%',/,'%%Title: ',A,/, &
        '%%Creator: SEM2DPACK Version ',A,/, &
        '%%Author: Jean-Paul Ampuero',/, &
        '%%BoundingBox: 0 0 612 792',/,'%%')

  end subroutine plot_header


!-----------------------------------------------------------------------
! Plot legends

  subroutine plot_legend()

  write(psunit,*) '0 setgray'
  write(psunit,*) '/Times-Roman findfont'
  write(psunit,*) '.5 CM SCSF'
  write(psunit,*) '24. CM 1.2 CM M'
  write(psunit,610) usoffset,it
  write(psunit,*) '%'
  write(psunit,*) '24. CM 1.95 CM M'
  write(psunit,600) usoffset,time
  write(psunit,*) '%'
  write(psunit,*) '24. CM 2.7 CM M'
  write(psunit,640) usoffset,maxfield
  write(psunit,*) '%'
  write(psunit,*) '/Times-Roman findfont'
  write(psunit,*) '.6 CM SCSF'
  if (color) write(psunit,*) '.4 .9 .9 RG'
  write(psunit,*) '11 CM 1.1 CM M'
  write(psunit,*) '(X) show'
  write(psunit,*) '%'
  write(psunit,*) '1.4 CM 9.5 CM M'
  write(psunit,*) 'currentpoint gsave translate 90 rotate 0 0 moveto'
  write(psunit,*) '(Z) show'
  write(psunit,*) 'grestore'
  write(psunit,*) '%'
  write(psunit,*) '/Times-Roman findfont'
  write(psunit,*) '.7 CM SCSF'
  if (color) write(psunit,*) '.8 0 .8 RG'
  write(psunit,*) '25.35 CM 18.9 CM M'
  write(psunit,*) usoffset,' CM 2 div neg 0 MR'
  write(psunit,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
  write(psunit,*) '(',stitle,') show'
  write(psunit,*) 'grestore'
  write(psunit,*) '26.45 CM 18.9 CM M'
  write(psunit,*) usoffset,' CM 2 div neg 0 MR'
  write(psunit,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
  write(psunit,*) '(SEM2DPACK '//version//' - Spectral Element Method) show'
  write(psunit,*) 'grestore'

  return

 600  format(F0.3,' neg CM 0 MR (Time =',EN12.3,' s) show')
 610  format(F0.3,' neg CM 0 MR (Time step = ',I0,') show')
 620  format(F0.3,' neg CM 0 MR (Cut =',F0.2,' \%) show')
 640  format(F0.3,' neg CM 0 MR (Max =',EN12.3,') show')

  end subroutine plot_legend

!-----------------------------------------------------------------------
! Draw the velocity model in background

  subroutine plot_model(PS)

  use elastic, only : elast_type,ELAST_inquire,ELAST_cpminmax,ELAST_csminmax

  integer, intent(in) :: PS   ! 1=P, 2=S, default=none

  double precision :: cmax,cmin,celem(grid%ngll,grid%ngll),c
  integer :: i,ip,j,jp,e

  if (PS==1) then
    call ELAST_cpminmax(elast,cmin,cmax)
  else if (PS==2) then
    call ELAST_csminmax(elast,cmin,cmax)
  else
    return
  endif

 ! quit if quasi-homogeneous velocity model, variation < 1%
  if ((cmax-cmin)/(cmin+cmax) < 0.02d0) then
    write(psunit,*) '%'
    write(psunit,*) '% no background : delta_v/v = ', 2d0*(cmax-cmin)/(cmax+cmin)
    write(psunit,*) '%'
  endif

  write(psunit,*) '%'
  write(psunit,*) '% background element fill'
  write(psunit,*) '%'

  do e=1,grid%nelem

    if (PS==1) then
      call ELAST_inquire(elast,e,cp=celem)
    else
      call ELAST_inquire(elast,e,cs=celem)
    endif

    do i=1,grid%ngll-1,isubsamp
      ip = min(grid%ngll,i+isubsamp)
      
      do j=1,grid%ngll-1,isubsamp
        jp = min(grid%ngll,j+isubsamp)
 
        write(psunit,500) point_scaled( grid%coord(:,grid%ibool(i,j,e)) )
        write(psunit,499) point_scaled( grid%coord(:,grid%ibool(ip,j,e)) )
        write(psunit,499) point_scaled( grid%coord(:,grid%ibool(ip,jp,e)) )
        write(psunit,499) point_scaled( grid%coord(:,grid%ibool(i,jp,e)) )
    
        c = (celem(i,j)-cmin)/(cmax-cmin)
        c = min( c*0.7 + 0.2 , 1.d0 ) ! rescale to avoid dark gray
        c = 1.d0 - c ! inverse scale: white = cmin, gray = cmax
    
        write(psunit,604) c
 
      enddo
    enddo

  enddo

 499  format(F0.2,1x,F0.2,' L')
 500  format(F0.2,1x,F0.2,' M')
 604  format('CP ',F0.2,' SF')

  end subroutine plot_model

!-----------------------------------------------------------------------
!-- Draw spectral element mesh

  subroutine plot_mesh(colors,elem_numbers,fill_domains)

  logical, intent(in) :: colors,elem_numbers,fill_domains
  double precision :: coord(NDIME,grid%ngll,grid%ngll) &
                     ,point(NDIME)
  double precision, pointer :: coorg(:,:)
  integer :: e,i,j,is,ir,imat,icol,ngnod

  write(psunit,*) '%'
  write(psunit,*) '% spectral element mesh'
  write(psunit,*) '%'

  ngnod = FE_GetNodesPerElement(grid%fem)

  do e=1,grid%nelem

    write(psunit,*) '% elem ',e
    write(psunit,*) 'MK'

    if (ngnod == 4) then ! tracer des droites si elements Q4
   
     ! get the coordinates of the control nodes
      coorg => FE_GetElementCoord(grid%fem,e)
      do i=1,ngnod
        write(psunit,601) point_scaled( coorg(:,i) )
      enddo
      write(psunit,601) point_scaled( coorg(:,1) )
      deallocate(coorg)
    
    else ! tracer des courbes si elements Q9

     ! get the coordinates of the GLL points
      do j=1,grid%ngll
      do i=1,grid%ngll
        coord(:,i,j) = grid%coord(:,grid%ibool(i,j,e))
      enddo
      enddo

      is=1
      do ir=1,grid%ngll
        write(psunit,601) point_scaled( coord(:,ir,is) )
      enddo

      ir=grid%ngll
      do is=2,grid%ngll
        write(psunit,601) point_scaled( coord(:,ir,is) )
      enddo
        
      is=grid%ngll
      do ir=grid%ngll-1,1,-1
        write(psunit,601) point_scaled( coord(:,ir,is) )
      enddo
        
      ir=1
      do is=grid%ngll-1,2,-1
        write(psunit,601) point_scaled( coord(:,ir,is) )
      enddo
      
    endif

    write(psunit,*) 'CO'

    if (fill_domains) then
     ! Use a different color for each domain
      imat = grid%tag(e)
      icol = mod(imat - 1,maxcolors) + 1
      write(psunit,600) RGB(:,icol)
    endif

    if (colors) then
      write(psunit,*) 'GC'
    else
      write(psunit,*) 'GG'
    endif

    ! write the element number
    if (elem_numbers) then

      coorg => FE_GetElementCoord(grid%fem,e)
      point = 0.25d0 * SUM( coorg(:,1:4) , dim =2 )
      deallocate(coorg)
      point = point_scaled(point)
      if (colors) write(psunit,*) '1 setgray'
      write(psunit,500) point
      write(psunit,502) e !-- write element number

    endif

  enddo

 500  format(F0.2,1x,F0.2,' M')
 502  format('fN (',I0,') Cshow')
 600  format(F0.2,1x,F0.2,1x,F0.2,' RG GF')
 601  format(F0.2,1x,F0.2)

  end subroutine plot_mesh

!-----------------------------------------------------------------------
!--  Draw the boundaries
   
  subroutine plot_boundaries

  double precision :: point1(NDIME),point2(NDIME)
  integer :: i,bce,ideb,iend

  write(psunit,*) '%'
  write(psunit,*) '% mesh boundaries'
  write(psunit,*) '%'

  write(psunit,*) '0.05 CM SLW'

  do i = 1,size(grid%bounds(:))
    write(psunit,*) '% boundary tag ',grid%bounds(i)%tag
  ! set color cyclically from the palette
    write(psunit,'( 3(F0.2,1x),"RG" )') RGB(:, mod(grid%bounds(i)%tag-1,maxcolors)+1 )