create_surfaces.f90

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              ||===\\                                                         | 
!              ||    \\                                                        |
!              ||     ||   //==\\   ||  ||   //==||  ||/==\\                   |
!              ||     ||  ||    ||  ||  ||  ||   ||  ||    ||                  |
!              ||    //   ||    ||  ||  ||  ||   ||  ||                        |
!              ||===//     \\==//    \\==\\  \\==\\  ||                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
!                                                                              |
!              CREATE_SURF    May. 2008                                        |
!                                                                              |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine create_surf(surface,is,debug)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! for a given surfaces passed as argument, this routine computes
! the 8 fields x,y,z,xn,yn,zn,r,s for each point of the surface
! according to its type and levelt, and the triangulation

!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use definitions
use randomodule
use constants

implicit none

type (sheet) :: surface
integer is
integer debug

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer :: levelt, surface_type, nproc,iproc,ierr,i,ij,j,ie,i1,i2,i3,k
integer :: seed, indix,err
integer :: ntmax,nhmax,npmax,nmax,nmode,nvmax,nnpnmax
integer :: nh,nohalt_hull,loc,nnn,nnlist,ntrilist
integer :: dmode,nt,nbmax,it,nhull,icircles,irect
integer,dimension(:),  allocatable :: hulltriangles,lt_work,ln_work
integer,dimension(:),  allocatable :: vis_tlist,vis_elist,add_tlist,nb
integer,dimension(:,:),allocatable :: vertices,neighbour
integer,dimension(:,:),allocatable :: nn
double precision :: sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10
double precision :: x1,x2,x3,y1,y2,y3,z1,z2,z3,xne,yne,zne,xyzn
double precision :: delta,eps,epsil,xtemp,ytemp
double precision :: delta_angle,delta_radius,angle,dist
double precision :: xx,yy,xini,xend,yini,yend,dx,dy
double precision,dimension(:),  allocatable :: field
double precision,dimension(:,:),allocatable :: points,centres
logical clockwise
character(len=4) :: cs
character ch
character*72  :: shift
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

INCLUDE 'mpif.h'

call mpi_comm_size (mpi_comm_world,nproc,ierr)
call mpi_comm_rank (mpi_comm_world,iproc,ierr)

shift=' '

if (surface%surface_type.lt.0) then
   !---------------------------------
   ! reading surface from ascii file
   !---------------------------------
   if (iproc.eq.0) then
      write (ch,'(i1)') -surface%surface_type
      write(*,*) shift//'surface-'//ch//'.txt'
      open(71,file='surface-'//ch//'.txt',status='old')
      read(71,*) surface%nsurface,surface%nt
      write(*,*) shift//'-> nb of points   :',surface%nsurface
      write(*,*) shift//'-> nb of triangles:',surface%nt
   end if
   call mpi_bcast (surface%nsurface,1,mpi_integer,0,mpi_comm_world,ierr)
   call mpi_bcast (surface%nt,1,mpi_integer,0,mpi_comm_world,ierr)
   deallocate (surface%x,surface%y,surface%z)
   deallocate (surface%xn,surface%yn,surface%zn)
   deallocate (surface%r,surface%s)
   deallocate (surface%icon)
   allocate (surface%x(surface%nsurface),surface%y(surface%nsurface),surface%z(surface%nsurface))
   allocate (surface%xn(surface%nsurface),surface%yn(surface%nsurface),surface%zn(surface%nsurface))
   allocate (surface%r(surface%nsurface),surface%s(surface%nsurface))
   allocate (surface%icon(3,surface%nt))
   if (iproc.eq.0) then
      write(*,*) surface%nsurface,surface%nt
      do i=1,surface%nsurface
         read (71,*) surface%x(i),surface%y(i),surface%z(i),surface%xn(i),surface%yn(i),surface%zn(i)
         surface%r(i)=surface%x(i)
         surface%s(i)=surface%y(i)
         if (surface%x(i)<0 .or. surface%x(i)>1 .or. surface%y(i)<0 .or. surface%y(i)>1 .or. surface%z(i)<0 .or. surface%z(i)>1 .and. iproc.eq.0)  &
            print *,surface%x(i),surface%y(i),surface%z(i) 
      enddo
      do i=1,surface%nt
         read (71,*) (surface%icon(j,i),j=1,3)
      enddo
      close (71)
   endif
   call mpi_bcast(surface%x,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%y,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%z,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%xn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%yn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%zn,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%r,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%s,surface%nsurface,mpi_double_precision,0,mpi_comm_world,ierr)
   call mpi_bcast(surface%icon,surface%nt*3,mpi_integer,0,mpi_comm_world,ierr)
else
   !----------------
   ! create surface
   !----------------

   seed=1234567

   sp01=surface%sp01
   sp02=surface%sp02
   sp03=surface%sp03
   sp04=surface%sp04
   sp05=surface%sp05
   sp06=surface%sp06
   sp07=surface%sp07
   sp08=surface%sp08
   sp09=surface%sp09
   sp10=surface%sp10
   surface_type=surface%surface_type
   levelt=surface%levelt

!   delta=2.d0**(-levelt) ! spacing between points on hull

   delta=1.d0 / (2.d0**levelt +1) ! opla

   epsil=(surface%stretch-1)/2.d0

   if (iproc.eq.0) write(8,*) 'create_surfaces: delta=',delta
   if (iproc.eq.0) write(8,*) 'create_surfaces: epsilon=',epsil

   !-----building the surface-------------
   ! type 1  : flat surface
   ! type 2  : rectangular emboss
   ! type 3  : convex spherical emboss 
   ! type 4  : concave spherical emboss
   ! type 5  : double rectangular emboss
   ! type 6  : sinus 
   ! type 7  : noisy surface 
   ! type 8  : double sinus 
   ! type 9  : cosinus sinus
   ! type 10 : slope
   ! type 11 : 2D hill

   select case(surface_type)

      case (1,2,5,6,7,8,9,10,11) 
!           if (surface%rand) then
              indix=1
              ! convex hull
              do i=1,2**levelt   +1
                 surface%x(indix)=dble(i-1)*delta
                 surface%y(indix)=0.d0
                 indix=indix+1
              end do
              do i=1,2**levelt   +1
                 surface%x(indix)=1.d0
                 surface%y(indix)=dble(i-1)*delta
                 indix=indix+1
              end do
              do i=1,2**levelt    +1
                 surface%x(indix)=1.d0-dble(i-1)*delta
                 surface%y(indix)=1.d0
                 indix=indix+1
              end do
              do i=1,2**levelt   +1
                 surface%x(indix)=0.d0
                 surface%y(indix)=1.d0-dble(i-1)*delta
                 indix=indix+1
              end do
              !interior nodes
              if(surface%rand) then
                 call random_number (surface%x(indix:surface%nsurface))
                 call random_number (surface%y(indix:surface%nsurface))
              else
                 surface%x(indix:surface%nsurface)=0.5d0
                 surface%y(indix:surface%nsurface)=0.5d0
              end if
              do j=2,2**levelt +1
                 do i=2,2**levelt  +1
                    surface%x(indix)=(((surface%x(indix)-0.5d0)*2.d0)*epsil+dble(i-1))*delta
                    surface%y(indix)=(((surface%y(indix)-0.5d0)*2.d0)*epsil+dble(j-1))*delta
                    indix=indix+1
                 enddo
              enddo
!           else
!              ! regular grid of points
!              indix=1
!              do j=1,2**levelt+1
!                 do i=1,2**levelt+1
!                    surface%x(indix)=float(i-1)*delta
!                    surface%y(indix)=float(j-1)*delta
!                    indix=indix+1 
!                 end do
!              end do
!           end if

      case (22)
           if (surface%rand) then
              indix=1
              ! convex hull
              do i=1,2**levelt
                 surface%x(indix)=(i-1)*delta
                 surface%y(indix)=0.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0
                 surface%y(indix)=(i-1)*delta
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0-(i-1)*delta
                 surface%y(indix)=1.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=0.d0
                 surface%y(indix)=1.d0-(i-1)*delta
                 indix=indix+1
              end do
              if (iproc.eq.0) write(8,*) 'nhull=',indix-1,'for surface ',is
              ! placing rows of points outside, on and outside the rectangle
              eps=0.001
              do irect=-5,3,2
                 xini=sp02+irect*eps
                 xend=sp03-irect*eps
                 yini=sp04+irect*eps
                 yend=sp05-irect*eps
                 dx=(xend-xini)/2**levelt
                 dy=(yend-yini)/2**levelt
                 do i=1,2**levelt
                    xx=xini+(i-1)*dx
                    yy=yini
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if
                 end do
                 do i=1,2**levelt
                    xx=xend
                    yy=yini+(i-1)*dy
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if
                 end do
                 do i=1,2**levelt
                    xx=xend-(i-1)*dx
                    yy=yend
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if
                 end do
                 do i=1,2**levelt
                    xx=xini
                    yy=yend-(i-1)*dy
                    if (xx.ge.0 .and. xx.le.1 .and. yy.ge.0 .and. yy.le.1) then
                       surface%x(indix)=xx
                       surface%y(indix)=yy
                       indix=indix+1
                    end if 
                 end do
              end do
              ! throwing remaining points at random
              do while (indix.le.surface%nsurface)
                 xx=ran(seed)
                 yy=ran(seed)
                 if ((xx<sp02-10*eps .or. xx>sp03+10*eps .or. yy<sp04-10*eps .or. yy>sp05+10*eps) .or. &
                     (xx>sp02+10*eps .and. xx<sp03-10*eps .and. yy>sp04+10*eps .and. yy<sp05-10*eps)) then 
                    surface%x(indix)=xx
                    surface%y(indix)=yy
                    indix=indix+1
                 end if
              end do
           else
              ! regular grid of points
              indix=1
              do j=1,2**levelt+1
                 do i=1,2**levelt+1
                 surface%x(indix)=float(i-1)*delta
                 surface%y(indix)=float(j-1)*delta
                 indix=indix+1 
                 end do
              end do
           end if   

      case (3,4) 
           if (surface%rand) then
              ! convex hull
              indix=1
              do i=1,2**levelt
                 surface%x(indix)=(i-1)*delta
                 surface%y(indix)=0.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0
                 surface%y(indix)=(i-1)*delta
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=1.d0-(i-1)*delta
                 surface%y(indix)=1.d0
                 indix=indix+1
              end do
              do i=1,2**levelt
                 surface%x(indix)=0.d0
                 surface%y(indix)=1.d0-(i-1)*delta
                 indix=indix+1
              end do
              nhull=indix-1 
              if (iproc.eq.0) write(8,*) 'nhull=',nhull,'on surface ',is
              ! placing rows of nhull points outside, in, and inside the circle
              delta_angle=2.d0*pi/nhull
              delta_radius=0.0025d0*sp04
              do icircles=-4,2
                 do i=0,nhull-1
                    angle=i*delta_angle
                    surface%x(indix)=sp02+(sp04+icircles*delta_radius)*cos(angle)
                    surface%y(indix)=sp03+(sp04+icircles*delta_radius)*sin(angle)
                    indix=indix+1
                 end do
              end do
              ! throwing the rest of the points at random provided they are not
              ! too close to the concentric circles so that it does not disturb
              ! the regular triangulation
              do while (indix.le.surface%nsurface)
                 xx=ran(seed)
                 yy=ran(seed)
                 dist=sqrt((xx-sp02)**2+(yy-sp03)**2)
                 if (dist>(sp04+10*delta_radius) .or. dist<(sp04-10*delta_radius)) then
                    surface%x(indix)=xx
                    surface%y(indix)=yy
                    indix=indix+1
                 end if
              end do
           else
              ! regular grid of points
              indix=1
              do j=1,2**levelt+1
                 do i=1,2**levelt+1
                    surface%x(indix)=float(i-1)*delta
                    surface%y(indix)=float(j-1)*delta
                    indix=indix+1 
                 end do
              end do
           end if

      case default
           call stop_run ('surface type not defined$')
   end select

   !-------------------------------------------|
   !-----building icon array-------------------|
   !-------------------------------------------|
   ! we use the NN library to compute the triangulation 
   ! of points for the surface
   ntmax=surface%nsurface*3
   nhmax=surface%nsurface
   npmax=surface%nsurface
   nnpnmax=100
   nmax=3*ntmax+npmax
   nvmax=ntmax
   allocate (points(2,surface%nsurface),stat=err) ; if (err.ne.0) call stop_run ('Error alloc points in create_surface$')
   allocate (field(surface%nsurface),stat=err)    ; if (err.ne.0) call stop_run ('Error alloc field in create_surface$')
   allocate (vertices(3,ntmax),stat=err)          ; if (err.ne.0) call stop_run ('Error alloc vertices in create_surface$')
   allocate (centres(3,ntmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc centress in create_surface$')
   allocate (neighbour(3,ntmax),stat=err)         ; if (err.ne.0) call stop_run ('Error alloc neighbour in create_surface$')
   allocate (hulltriangles(nhmax),stat=err)       ; if (err.ne.0) call stop_run ('Error alloc hulltriangles in create_surface$')
   allocate (vis_tlist(nvmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc vis_tlist in create_surface$')
   allocate (vis_elist(nvmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc vis_elist in create_surface$')
   allocate (add_tlist(nvmax),stat=err)           ; if (err.ne.0) call stop_run ('Error alloc add_tlist in icreate_surface$')
   allocate (lt_work(ntmax),stat=err)             ; if (err.ne.0) call stop_run ('Error alloc lt_work in create_surface$')
   allocate (ln_work(surface%nsurface),stat=err)  ; if (err.ne.0) call stop_run ('Error alloc ln_work in create_surface$')
   points(1,:)=surface%x
   points(2,:)=surface%y
   dmode=-2
   nmode=0
   clockwise=.true.
   eps=0.d0
   call nn2d_setup (surface%nsurface,ntmax,nhmax,npmax,nnpnmax,nmax, &
                    points,dmode,nmode,clockwise,field,nt,vertices, &
                    centres,neighbour,nh,hulltriangles,nohalt_hull, &
                    loc,nnn,nnlist,ntrilist, &
                    eps,nvmax,vis_tlist,vis_elist,add_tlist, &
                    lt_work,ln_work)
   if (iproc.eq.0) write (8,*) 'surface ',is,' has ',nt,' triangles'

   surface%icon(1:3,1:surface%nt)=vertices(1:3,1:nt)

   deallocate (points)
   deallocate (field)
   deallocate (vertices)
   deallocate (centres)
   deallocate (neighbour)
   deallocate (hulltriangles)
   deallocate (vis_tlist)
   deallocate (vis_elist)
   deallocate (add_tlist)
   deallocate (lt_work)
   deallocate (ln_work)

   !----------------------------------------------|
   !----computing z coordinate--------------------|
   !----------------------------------------------|

   call zpoints (surface%nsurface,surface%x,surface%y,surface%z,surface_type, &
                 sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10)

   !----------------------------------------------|
   !----computing the normals---------------------|  
   !----------------------------------------------|

   call compute_normals (surface%nsurface,surface%x,surface%y,surface%z, &
                         surface%nt,surface%icon,surface%xn,surface%yn,surface%zn)

   if (debug>=2) call output_surf (surface,is,'init',-1,-1) 

end if

end subroutine create_surf




!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|

subroutine zpoints (ns,x,y,z,surface_type, &
                    sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10)

!------------------------------------------------------------------------------|
!(((((((((((((((( Purpose of the routine  ))))))))))))))))))))))))))))))))))))))
!------------------------------------------------------------------------------|
! knowing the two-dimensional spatial distribution of points on the plane (x,y)
! this routine computes the z-coordinate accordingly to the surface type.
! ns is the number of points
! x,y,z are the coordinates of the points
! surface_type is the type of surface under consideration
! sp01..sp10 are the surface parameters (not all are necessarily used)
!------------------------------------------------------------------------------|
!((((((((((((((((  declaration of the subroutine arguments  ))))))))))))))))))))
!------------------------------------------------------------------------------|

use constants

implicit none

integer ns
double precision x(ns),y(ns),z(ns)
integer surface_type
double precision sp01,sp02,sp03,sp04,sp05,sp06,sp07,sp08,sp09,sp10

!------------------------------------------------------------------------------|
!(((((((((((((((( declaration of the subroutine internal variables )))))))))))))
!------------------------------------------------------------------------------|

integer i
double precision delta,dist,ran

!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|


select case(surface_type)
   case (1)
        ! a flat surface, 
        ! sp01 is the z level
        z=sp01
   case (2)
        ! rectangular emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x1,x2 
        ! sp04 and 05 are y1,y2 
        ! sp06 is the thickness 
        do i=1,ns
           if ((x(i)-sp02)*(x(i)-sp03).le.0.d0 .and. (y(i)-sp04)*(y(i)-sp05).le.0.d0) then
              z(i)=sp01-sp06
           else
              z(i)=sp01
           end if
        end do
   case (3)
        ! convex spherical emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x0,y0 
        ! sp04 is the radius 
        do i=1,ns
           dist=sqrt((x(i)-sp02)**2+(y(i)-sp03)**2)
           if (dist.le.sp04) then
              z(i)=sp01+sp04*sin(acos(dist/sp04))
           else
              z(i)=sp01
           endif
       end do
   case (4)
        ! convex spherical emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x0,y0 
        ! sp04 is the radius 
        do i=1,ns
           dist=sqrt((x(i)-sp02)**2+(y(i)-sp03)**2)
           if (dist.le.sp04) then
              z(i)=sp01-sp04*sin(acos(dist/sp04))
           else
              z(i)=sp01
           endif
       end do
   case (5)
        ! double rectangular emboss, 
        ! sp01 is the z level
        ! sp02 and 03 are x1,x2 
        ! sp04 and 05 are x3,x4 
        ! sp06 and 07 are y1,y2 
        ! sp08 and 09 are y3,y4 
        ! sp10 is the thickness 
         z=sp01
         do i=1,ns
            if ((x(i)-sp02)*(x(i)-sp03).le.0.d0 .and. (y(i)-sp06)*(y(i)-sp07).le.0.d0) then
               z(i)=z(i)-sp10
            endif
            if ((x(i)-sp04)*(x(i)-sp05).le.0.d0 .and. (y(i)-sp08)*(y(i)-sp09).le.0.d0) then
              z(i)=z(i)-sp10
           endif
        end do
   case (6)
        ! a sinus, 
        ! sp01 is the z level
        ! sp02 is the wavelength 
        ! sp03 is the amplitude 
        do i=1,ns
           z(i)=sp01+sp03*sin(x(i)*2.d0*pi/sp02)
        end do
   case (7)
        ! a noisy surface, 
        ! sp01 is the z level
        ! sp02 is the noise amplitude
        do i=1,ns
           call random_number (ran)
           z(i)=sp01+ran*sp02
        end do
   case (8)
        ! a double sinus.
        ! sp01 is the z level
        ! sp02 is the x-wavelength 
        ! sp03 is the x-amplitude 
        ! sp04 is the y-wavelength 
        ! sp05 is the y-amplitude
        do i=1,ns 
           z(i)=sp01+sp03*sin(x(i)*2.*pi/sp02)+sp05*sin(y(i)*2.*pi/sp04)
        end do
   case (9)
        ! a cosinus, 
        ! sp01 is the z level
        ! sp02 is the wavelength 
        ! sp03 is the amplitude 
        do i=1,ns
           z(i)=sp01+sp03*cos(x(i)*2.d0*pi/sp02)
        end do
   case (10)
        ! a 2D embankment
        ! sp01 is z0
        ! sp02 is y0
        ! sp03 is psi
        ! sp04 is the thickness
        do i=1,ns 
           if (y(i)<sp02) then
              z(i)=sp01
           else
              z(i)=sp01+min((y(i)-sp02)*tan(sp03*pi/180.d0),sp04)
           end if
        end do
   case (11)
       ! a 2D hill
       ! sp01 is z level
       ! sp02 height
       ! sp03 is y0
       ! sp04 is width
        do i=1,ns 
           if (y(i)>sp03-sp04/2.d0 .and. y(i)<sp03+sp04/2.d0) then
              z(i)=sp01+sp02*(sin(pi*((y(i)-sp03)/sp04+0.5d0)))**2
           else
              z(i)=sp01
           end if
        end do


   case default
         call stop_run ('surface type not defined$')
end select


return
end