plot_postscript.f90

Spectral Element Method for wave propagation and rupture dynamics.

  ! draw a straight segment for each boundary element
    do bce = 1,grid%bounds(i)%nelem
      ideb   = grid%bounds(i)%node( grid%bounds(i)%ibool(1,bce) )
      point1 = point_scaled( grid%coord(:,ideb) )
      iend   = grid%bounds(i)%node( grid%bounds(i)%ibool(grid%ngll,bce) )
      point2 = point_scaled( grid%coord(:,iend) )
      write(psunit,'( 2(2(F0.2,1x),A) )') point1,'M ',point2,'L ST'
    enddo
  enddo

  write(psunit,*) '0 setgray'
  write(psunit,*) '0.01 CM SLW'

  end subroutine plot_boundaries


!---------------------------------------------------------------------
! Plot a vector field defined at each node
!
  subroutine plot_vect()

  double precision :: point(NDIME),Uinterp(NDIME),factor
  double precision, pointer :: coorg(:,:)
  double precision, allocatable :: vloc(:,:)
  integer :: e,i,j,k,l,ipoin

  if (maxfield>0d0) then
    write(psunit,*) '%'
    write(psunit,*) '% vector field'
    write(psunit,*) '%'
  else
    write(psunit,*) '%'
    write(psunit,*) '% empty vector field'
    write(psunit,*) '%'
    return
  endif

 ! scale the vector field
  if (ScaleField > 0.d0) then
    factor = 1.d0/ScaleField
  else
    factor = 1.d0/maxfield
  endif

  if(color) then
    write(psunit,*) 'Colvects'
  else
    write(psunit,*) '0 setgray'
  endif

  if (interpol) then  !-- interpolate the vector field

    allocate(vloc(grid%ngll*grid%ngll,size(vfield,2)))
    do e=1,grid%nelem

      coorg => FE_GetElementCoord(grid%fem,e)

      k=1
      do i=1,grid%ngll
      do j=1,grid%ngll
        vloc(k,:) = vfield(grid%ibool(i,j,e),:)
        k=k+1
      enddo
      enddo

      do j=1,DisplayPts
      do i=1,DisplayPts
        point = matmul( coorg, fe_interp(:,i,j) )
        Uinterp = matmul(se_interp(:,i,j),vloc)
        call plot_single_vect(factor*Uinterp,point)
      enddo
      enddo

      deallocate(coorg)

    enddo
    deallocate(vloc)

  else !-- Plot vectors at mesh nodes

    do ipoin=1,grid%npoin
      call plot_single_vect(factor*vfield(ipoin,:), grid%coord(:,ipoin) )
    enddo

  endif

  write(psunit,*) '0 setgray'


  end subroutine plot_vect


!----------------------------------------------------------------------
! Plots scalar values defined element wise
  subroutine plot_efield

  double precision :: min_efield,max_efield,scale,x1
  double precision, pointer :: coorg(:,:)
  integer :: e

  min_efield = minval(efield)
  max_efield = maxval(efield)
  scale = max_efield - min_efield
  if (scale < abs(min_efield+max_efield)*1d-3) then
    write(psunit,*) '%'
    write(psunit,*) '% empty element-wise scalar field'
    write(psunit,*) '%'
    return    
  else
    write(psunit,*) '%'
    write(psunit,*) '% element-wise scalar field'
    write(psunit,*) '%'
    scale=1d0/scale
  endif

  do e=1,grid%nelem

    x1 = (efield(e)-min_efield)*scale
    !x1 = min( x1*0.7 + 0.2 , 1. ) ! rescale to avoid too dark a gray
    x1 = 1.d0 - x1 ! invert the scale: white = min, gray = max

    coorg => FE_GetElementCoord(grid%fem,e)
    write(psunit,'(9(F0.2,1x),"FQ")') x1 &
                     ,point_scaled( coorg(:,1) ) &
                     ,point_scaled( coorg(:,2) ) &
                     ,point_scaled( coorg(:,3) ) &
                     ,point_scaled( coorg(:,4) )
    deallocate(coorg)

  enddo
    
  end subroutine plot_efield

!----------------------------------------------------------------------
! Plots nodal scalar values 
! Actually, the average of the 4 neighbouring GLL nodes
  subroutine plot_scal(s)

  double precision, intent(in) :: s(:)

  double precision, allocatable :: spl(:,:), cpl(:,:,:)
  double precision :: stmp(DisplayPts*DisplayPts)
  double precision, pointer :: coorg(:,:)
  double precision :: val, scale
  integer :: e,i,j,k,iglob,npl

  if (maxfield>0d0) then
    write(psunit,*) '%'
    write(psunit,*) '% scalar field'
    write(psunit,*) '%'
  else
    write(psunit,*) '%'
    write(psunit,*) '% empty scalar field'
    write(psunit,*) '%'
    return
  endif

  if (ScaleField>0d0) then
    scale = 1d0/ScaleField
  else
    scale = 1d0/maxval(abs(s))
  endif

  if (interpol) then
    npl = DisplayPts
  else
    npl = grid%ngll
  endif
  allocate(spl(npl,npl), cpl(npl,npl,2) )

  do e=1,grid%nelem

    if (interpol) then
      coorg => FE_GetElementCoord(grid%fem,e)
      k=0
      do j=1,grid%ngll
      do i=1,grid%ngll
        k=k+1
        stmp(k) = s(grid%ibool(i,j,e))
      enddo
      enddo
      do j=1,npl
      do i=1,npl
        cpl(i,j,:) = matmul( coorg, fe_interp(:,i,j) )
        spl(i,j) = sum( se_interp(:,i,j)*stmp)
      enddo
      enddo

    else
      do j=1,grid%ngll
      do i=1,grid%ngll
        iglob = grid%ibool(i,j,e)
        cpl(i,j,:) = grid%coord(:,iglob)
        spl(i,j) = s(iglob)
      enddo
      enddo
    endif
    
    do j=1,npl-1
    do i=1,npl-1
      val = 0.25d0 *(spl(i,j)+spl(i+1,j)+spl(i+1,j+1)+spl(i,j+1)) *scale
      val = min(val,1d0)
      val = max(val,-1d0)
     ! at this point, val in [-1:1]

      if (color) then ! min = -1 = blue, max = 1 = red
        if (abs(val)<0.01d0) cycle   ! skip if white
        write(psunit,'(11(F0.2,1x),"FQC")') &
          1d0+min(val,0d0),1d0-abs(val),1d0-max(val,0d0) &
          ,point_scaled( cpl(i,j,:) ) &
          ,point_scaled( cpl(i+1,j,:) ) &
          ,point_scaled( cpl(i+1,j+1,:) ) &
          ,point_scaled( cpl(i,j+1,:) )

      else            ! min = 1 = white, max = 0 = gray
        if (val< -0.99d0) cycle   ! skip if white
        val = 0.5d0*(1d0-val) 
        write(psunit,'(9(F0.2,1x),"FQ")') val &
          ,point_scaled( cpl(i,j,:) ) &
          ,point_scaled( cpl(i+1,j,:) ) &
          ,point_scaled( cpl(i+1,j+1,:) ) &
          ,point_scaled( cpl(i,j+1,:) )

      endif

    enddo
    enddo
  enddo

  deallocate(spl, cpl )

  end subroutine plot_scal

!----------------------------------------------------------------------
! Places a vector VECT at location ORIGIN
  subroutine plot_single_vect(VECT,ORIGIN)

    double precision, dimension(NDIME), intent(in) :: VECT,ORIGIN

    double precision, parameter :: cutvect = 0.01d0 ! min arrow length / scale factor
    double precision, parameter :: sizemax = 1d0 ! Maximum arrow length (in cm)
    double precision, parameter :: ArrowHeadAngle = 20.d0 ! Angle between arrow body and head (degrees)
    double precision, parameter :: ArrowHBRatio = 0.4d0 ! Arrow head to body ratio

    double precision :: d,d1,d2,dummy,theta,thetaup,thetadown &
                       ,x2,z2,xa,za,xb,zb,convert
    
   ! pi conversion
    convert = 3.141592653589793d0/180.d0

    x2 = VECT(1)*sizemax
    z2 = VECT(2)*sizemax

    d = sqrt(x2*x2 + z2*z2)

! ignore small vectors
    if (d < cutvect*sizemax) return

    d1 = d * ArrowHBRatio
    d2 = d1 * cos(ArrowHeadAngle*convert)

    dummy = x2/d
    dummy = min(dummy,0.9999d0)
    dummy = max(dummy,-0.9999d0)
    theta = acos(dummy)

    if(z2 < 0.d0) theta = 360.d0*convert - theta
    thetaup = theta - ArrowHeadAngle*convert
    thetadown = theta + ArrowHeadAngle*convert

! plot vector
    x2 = x2 * centim
    z2 = z2 * centim
    xa = -d2*cos(thetaup)
    za = -d2*sin(thetaup)
    xa = xa * centim
    za = za * centim
    xb = -d2*cos(thetadown)
    zb = -d2*sin(thetadown)
    xb = xb * centim
    zb = zb * centim
    write(psunit,700) xb,zb,xa,za,x2,z2,point_scaled(ORIGIN)

  700  format(8(F0.1,1x),'F') 
    
  end subroutine plot_single_vect

!------------------------------------------------------------------
!
  function point_scaled(coord)

    double precision, intent(in) :: coord(NDIME)
    double precision :: point_scaled(NDIME)
   
   ! origin in the page (cm)
    double precision, parameter :: PageX0 = 2.2d0, PageZ0 = 2.2d0

    point_scaled(1) = (coord(1)-xmin)*rapp_page + PageX0
    point_scaled(2) = (coord(2)-zmin)*rapp_page + PageZ0
    point_scaled    = point_scaled * centim

  end function point_scaled

!------------------------------------------------------------------
!
  subroutine plot_sources

  double precision :: coord(NDIME)
  integer :: i

  write(psunit,*) '% sources'
  do i=1,size(src)
    call SO_inquire(src(i),coord=coord)
    write(psunit,510) point_scaled(coord)
    write(psunit,*) 'Cross'
  enddo

  510  format(F0.1,1x,F0.1,' M')
  end subroutine plot_sources

!------------------------------------------------------------------
!
  subroutine plot_receivers
  
  double precision, pointer :: posrec(:,:)
  integer :: nrec,i
  !-- limite pour afficher des points a la place des recepteurs
  integer, parameter :: ndots = 10

  call REC_inquire(rec,coord=posrec)
  nrec = size(posrec,2)
  write(psunit,*) '% start receiver line'
  do i=1,nrec
    write(psunit,510) point_scaled(posrec(:,i))
    if(nrec > ndots.and.i /= 1.and.i /= nrec) then
      write(psunit,*) 'VDot'
    else
      write(psunit,*) 'Losange'
    endif
  enddo
  write(psunit,*) '% end receiver line'

  510  format(F0.1,1x,F0.1,' M')
  end subroutine plot_receivers

end subroutine PLOT_PS

end module plot_postscript