program SubductionGeometry
implicit none
integer level,nx,ny,i,j,nx1,nxs,ny1,ny2
integer i1,i2,i3,ij,k
double precision x1,x2,x3,y1,y2,y3,z1,z2,z3
double precision xne,yne,zne,xyzn
integer nnode,nnode0,nxt,nxb,nnodemax,nelem,nelemmax
double precision xl1,yl1,yl2,xs,dip,length,thick,pi,zfree
double precision x,y,z,xsp,dist,depth
double precision, dimension(:),allocatable::xp,yp,zp,xn,yn,zn
integer,dimension(:,:),allocatable::icon
pi=atan(1.)*4.
! xl1 is where the plate starts in the x-direction;
! xs is where the subduction starts; dip is the dip of the subducting
! slab; length is it length (along the dipping plane); thick is the
! plate thickness; yl1 is where the plate starts in the y direction and
! yl2 is where it ends
xl1=.4376d0
xs=.5624d0
dip=60.d0
dip=dip*pi/180.d0
thick=.01d0
depth=.1d0
length=(depth-thick-.0001)/sin(dip)
yl1=.4376d0
yl2=.5624d0
xsp=xs-thick*tan(dip/2.d0)
zfree=1.d0
level=8
nx=2**level+1
ny=2**level+1
nnodemax=nx*ny*2
nelemmax=nnodemax*3
allocate (xp(nnodemax),yp(nnodemax),zp(nnodemax))
allocate (xn(nnodemax),yn(nnodemax),zn(nnodemax))
allocate (icon(3,nelemmax))
nx1=int(xl1*nx)
ny1=int(yl1*ny)
ny2=int(yl2*ny)
nxs=int(nx*xs)
nxb=int(length*nx)
nxt=nxb+1
nnode=0
nelem=0
do j=1,ny1
! no plate
y=float(j-1)/float(ny1-1)*yl1
do i=1,nx1
x=float(i-1)/float(nx1-1)*xl1
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny1) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
enddo
do i=nx1+1,nxs
x=xl1+(xsp-xl1)*float(i-nx1)/float(nxs-nx1)
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny1) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
enddo
do i=nxs,nx
x=xs+(1.-xs)*float(i-nxs)/float(nx-nxs)
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (i.ne.nx.and.j.ne.ny1) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
enddo
enddo
do j=ny1,ny2
! no plate
do i=1,nx1
x=float(i-1)/float(nx1-1)*xl1
y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
enddo
! horizontal plate
do i=nx1,nxs+1
x=xl1+(xsp-xl1)*float(i-nx1)/float(nxs+1-nx1)
y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
z=zfree-thick
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
if (i.ne.nxs+1.and.j.eq.ny1) call add_element (nelem,icon,nnode-2-nx,nnode-1-nx,nnode+1,nnode,nelemmax)
if (i.ne.nxs+1.and.j.eq.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+4+nx+nxb+nxt,nnode+3+nx+nxb+nxt,nelemmax)
enddo
! bottom dipping plate
do i=1,nxb+1
x=xsp+length*cos(dip)*float(i)/float(nxb+1)
y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
z=zfree-thick-length*sin(dip)*float(i)/float(nxb+1)
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
if (j.eq.ny1) call add_element (nelem,icon,nnode+nxb+nxt-2*i+3,nnode+nxb+nxt-2*i+2,nnode,nnode-1,nelemmax)
if (j.eq.ny2) call add_element (nelem,icon,nnode-1,nnode,nnode+nxb+nxt-2*i+2,nnode+nxb+nxt-2*i+3,nelemmax)
enddo
! top of dipping plate
do i=1,nxt
x=xs+length*cos(dip)-length*cos(dip)*float(i-1)/float(nxt)
y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
z=zfree-length*sin(dip)+length*sin(dip)*float(i-1)/float(nxt)
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
enddo
! no plate
do i=nxs,nx
x=xs+(1.-xs)*float(i-nxs)/float(nx-nxs)
y=yl1+float(j-ny1)/float(ny2-ny1)*(yl2-yl1)
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (i.ne.nx.and.j.ne.ny2) call add_element (nelem,icon,nnode,nnode+1,nnode+5+nx+nxb+nxt,nnode+4+nx+nxb+nxt,nelemmax)
enddo
enddo
do j=ny2,ny
! no plate
y=yl2+float(j-ny2)/float(ny-ny2)*(1.-yl2)
do i=1,nx1
x=float(i-1)/float(nx1-1)*xl1
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
enddo
do i=nx1+1,nxs
x=xl1+(xsp-xl1)*float(i-nx1)/float(nxs-nx1)
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (j.ne.ny) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
enddo
do i=nxs,nx
x=xs+(1.-xs)*float(i-nxs)/float(nx-nxs)
z=zfree
call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
if (i.ne.nx.and.j.ne.ny) call add_element (nelem,icon,nnode,nnode+1,nnode+2+nx,nnode+1+nx,nelemmax)
enddo
enddo
print*,nnode,nelem
i1=1
111 continue
do i=i1,nnode
do j=i+1,nnode
dist=(xp(i)-xp(j))**2+(yp(i)-yp(j))**2+(zp(i)-zp(j))**2
if (dist.lt.1.d-10) then
print*,'removing ',j
do ij=j+1,nnode
xp(ij-1)=xp(ij)
yp(ij-1)=yp(ij)
zp(ij-1)=zp(ij)
enddo
do ij=1,nelem
do k=1,3
if (icon(k,ij).eq.j) icon(k,ij)=i
enddo
enddo
do ij=1,nelem
do k=1,3
if (icon(k,ij).gt.j) icon(k,ij)=icon(k,ij)-1
enddo
enddo
nnode=nnode-1
i1=i
goto 111
endif
enddo
enddo
print*,nnode,nelem
print*,minval(icon(:,1:nelem)),maxval(icon(:,1:nelem))
!add normals
do i=1,nelem
i1=icon(1,i)
i2=icon(2,i)
i3=icon(3,i)
if (i1.eq.i2 .and. i2.eq.i3) print*,'element ',i,' collapsed: ',icon(:,i)
x1=xp(i1)
x2=xp(i2)
x3=xp(i3)
y1=yp(i1)
y2=yp(i2)
y3=yp(i3)
z1=zp(i1)
z2=zp(i2)
z3=zp(i3)
xne=(y2-y1)*(z3-z1)-(y3-y1)*(z2-z1)
yne=(z2-z1)*(x3-x1)-(z3-z1)*(x2-x1)
zne=(x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)
xyzn=sqrt(xne**2+yne**2+zne**2)
xne=xne/xyzn
yne=yne/xyzn
zne=zne/xyzn
do j=1,3
ij=icon(j,i)
xn(ij)=xn(ij)+xne
yn(ij)=yn(ij)+yne
zn(ij)=zn(ij)+zne
enddo
enddo
do i=1,nnode
xyzn=sqrt(xn(i)**2+yn(i)**2+zn(i)**2)
xn(i)=xn(i)/xyzn
yn(i)=yn(i)/xyzn
zn(i)=zn(i)/xyzn
enddo
open (71,file='surface-1.txt',status='unknown')
write (71,*) nnode,nelem
do i=1,nnode
write (71,*) xp(i),yp(i),zp(i),xn(i),yn(i),zn(i)
enddo
do i=1,nelem
write (71,*) (icon(j,i),j=1,3)
enddo
close (71)
open(unit=30,file='Mesh.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 ',nnode,' float'
do i=1,nnode
write(30,'(3f16.11)') xp(i),yp(i),zp(i)
enddo
write(30,'(A6, 2I10)') 'CELLS ',nelem,4*nelem
do i=1,nelem
write(30,'(9I10)')3,icon(1:3,i)-1
enddo
write(30,'(A11, I10)') 'CELL_TYPES ',nelem
do i=1,nelem
write(30,'(I2)')5 ! octree (8 nodes)
enddo
write(30,'(a11,i10)')'POINT_DATA ',nnode
write(30,'(a)')'VECTORS norm float'
do i=1,nnode
write(30,'(3f18.13)')xn(i),yn(i),zn(i)
enddo
close(30)
deallocate (xp,yp,zp,icon)
end
!--------
subroutine add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
double precision xp(nnodemax),yp(nnodemax),zp(nnodemax)
double precision x,y,z
nnode=nnode+1
if (nnode.gt.nnodemax) then
print*,nnode,nnodemax
stop 'needs to increase nnodemax'
endif
xp(nnode)=x
yp(nnode)=y
zp(nnode)=z
return
end
!----------
subroutine add_element (nelem,icon,i1,i2,i3,i4,nelemmax)
integer nelem,nelemmax,icon(3,nelemmax),i1,i2,i3,i4
nelem=nelem+1
if (nelem.gt.nelemmax) then
print*,nelem,nelemmax
stop 'nelemmax needs to be increased'
endif
icon(1,nelem)=i1
icon(2,nelem)=i2
icon(3,nelem)=i3
nelem=nelem+1
if (nelem.gt.nelemmax) then
print*,nelem,nelemmax
stop 'nelemmax needs to be increased'
endif
icon(1,nelem)=i1
icon(2,nelem)=i3
icon(3,nelem)=i4
return
end