-
Dave Whipp authoredDave Whipp authored
create_surfaces.f90 23.24 KiB
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! ||===\\ |
! || \\ |
! || || //==\\ || || //==|| ||/==\\ |
! || || || || || || || || || || |
! || // || || || || || || || |
! ||===// \\==// \\==\\ \\==\\ || |
! |
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
! |
! 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
integer :: dmode,nt,nbmax,it,nhull,icircles,irect
integer,dimension(:), allocatable :: hulltriangles
integer,dimension(:), allocatable :: vis_tlist,vis_elist,add_tlist,nb
integer,dimension(:,:),allocatable :: vertices,neighbour
integer,dimension(:,:),allocatable :: nn
logical*1,dimension(:),allocatable :: lt_work,ln_work
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
integer nnn(1),nnlist(1),ntrilist(1)
!------------------------------------------------------------------------------|
!------------------------------------------------------------------------------|
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
loc=1
nohalt_hull=0
clockwise=.true.
eps=1.d-8
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