BallGeometry.f90

      program BallGeometry

      implicit none

      integer level,i,j,nrad,nz,nradp,ip
      integer i1,i2,i3,i4,ij,k,np
      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,nnodep
      double precision xc,yc,zc,rad,radp,phi,theta
      double precision x,y,z,xsp,dist,pi
      double precision, dimension(:),allocatable::xp,yp,zp,xn,yn,zn
      integer,dimension(:,:),allocatable::icon

      pi=atan(1.)*4.

! xc yc zc is the centre of the ball
! rad is radius

      xc=.501d0
      yc=.501d0
      zc=.495d0
      xc=.5d0
      yc=.5d0
!      zc=.9d0
      rad=.025d0
! reference size
      rad=.05d0

      level=9

      nz=int(2.d0**level*rad*pi)+1
      if ((nz/2)*2.eq.nz) nz=nz+1
      nnodemax=nz*nz*2
      nelemmax=nnodemax*3

      allocate (xp(nnodemax),yp(nnodemax),zp(nnodemax))
      allocate (xn(nnodemax),yn(nnodemax),zn(nnodemax))
      allocate (icon(3,nelemmax))

      nnode=0
      nelem=0

      ip=1
      do j=1,nz
      theta=-pi/2.+pi*float(j-1)/float(nz-1)
      radp=rad*cos(theta)
      z=rad*sin(theta)
      nrad=1+int(radp*pi*2.d0*2.d0**level)
      print*,nz,nrad
      do i=1,nrad
      phi=2.d0*pi*float(i-1)/float(nrad)
      x=radp*sin(phi)
      y=radp*cos(phi)
      call add_point (nnode,x,y,z,xp,yp,zp,nnodemax)
! bottom half
        if (j.ne.1.and.j.le.nz/2+1) then
        i1=nnode
        i2=nnode+1
        if (i.eq.nrad) i2=nnodep+nradp+1
        i3=nnodep+int((i*nradp)/nrad)+1
        ip=nnodep+int(((i-1)*nradp)/nrad)+1
        if (i3.ge.nnodep+nradp+1) i3=nnodep+1
        call add_element (nelem,icon,i1,i2,i3,nelemmax)
        if (i3.ne.ip) call add_element (nelem,icon,i1,i3,ip,nelemmax)
        ip=i3
! top half (minus very top)
        elseif (j.ne.1.and.j.ne.nz) then
        i1=nnode
        i2=nnode+1
        if (i.eq.nrad) i2=nnode-nrad+1
        i3=nnodep+int((i*nradp)/nrad)+1
        if (i.eq.1) i3=nnodep+2
        np=1
        if (i.ne.1) np=i3-ip
        do k=1,np
        call add_element (nelem,icon,i3-k,i1,i3-k+1,nelemmax)
        enddo
        if (i.eq.nrad) call add_element (nelem,icon,i2,nnodep+1,nnode-nrad,nelemmax)
        call add_element (nelem,icon,i1,i2,i3,nelemmax)
        ip=i3
! very top
        elseif (j.eq.nz) then
        do k=1,nradp
        i1=nnodep+k
        i2=i1+1
        if (i2.eq.nnode) i2=nnodep+1
        i3=nnode
        call add_element (nelem,icon,i1,i3,i2,nelemmax)
        enddo
        endif
      enddo
      nradp=nrad
      nnodep=nnode-nrad
      enddo

      print*,nnode,nelem

!add normals

112   continue
      do i=1,nelem
      i1=icon(1,i)
      i2=icon(2,i)
      i3=icon(3,i)
      if (i1.eq.i2 .or. i2.eq.i3 .or. i3.eq.i1)then
      print*,'element ',i,' collapsed: ',icon(:,i)
        do j=i+1,nelem
        icon(:,j-1)=icon(:,j)
        enddo
      nelem=nelem-1
      goto 112
      endif
      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

      xp=xp+xc
      yp=yp+yc
      zp=zp+zc

      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,nelemmax)

      integer nelem,nelemmax,icon(3,nelemmax),i1,i2,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)=i2
      icon(3,nelem)=i3

      return
      end