Modified refine surface to remove points along edges that are closer together than the levelmax+4

src/refine_surface.f90

@@ -57,6 +57,7 @@ implicit none
 
 type (parameters) params
 type (sheet) surface 
+type (sheet) surface0 
 type (thread) threadinfo
 integer nadd,nrem
 integer is,istep
@@ -79,7 +80,7 @@ integer ie,j,jedge,kp,kpp,k,inode,inodep,inodepp,iadd,irem
 integer ntriangle
 integer nedge
 integer nedgen,nsurfacen
-double precision dist1,distmax,prod,distmin
+double precision dist1,distmax,prod,distmin,eps
 double precision, external :: dist
 
 
@@ -98,6 +99,7 @@ call mpi_comm_size (mpi_comm_world,nproc,ierr)
 call mpi_comm_rank (mpi_comm_world,iproc,ierr)
          
 shift=' '
+eps=1.d-8
 
 !------------------------
 ! build edge information 
@@ -108,7 +110,6 @@ ntriangle=surface%nt
 nedge=nsurfacen+ntriangle-1
 nnmax=12 ! nnmax is maximum number of neighbours in triangulation (should be around 6 on average)
 
-      
 allocate (ed(nedge),stat=threadinfo%err); call heap (threadinfo,'ed','refine_surface',size(ed),'int',+1)
 allocate (nedgepernode(nsurfacen),stat=threadinfo%err)
 call heap (threadinfo,'nedgepernode','refine_surface',size(nedgepernode),'int',+1)
@@ -119,40 +120,39 @@ call heap (threadinfo,'nodeedgenumber','refine_surface',size(nodeedgenumber),'in
 
 ed%t1=0
 ed%t2=0
-
 nedgepernode=0
-
 iedge=0
+
 do ie=1,ntriangle
-   do k=1,3
-   inode=surface%icon(k,ie)
-   kp=mod(k,3)+1
-   kpp=mod(k+1,3)+1
-   inodep=surface%icon(kp,ie)
-   inodepp=surface%icon(kpp,ie)
-   do j=1,nedgepernode(inodep)
+  do k=1,3
+    inode=surface%icon(k,ie)
+    kp=mod(k,3)+1
+    kpp=mod(k+1,3)+1
+    inodep=surface%icon(kp,ie)
+    inodepp=surface%icon(kpp,ie)
+    do j=1,nedgepernode(inodep)
       if (nodenodenumber(j,inodep).eq.inode) then
-         jedge=nodeedgenumber(j,inodep)
-         ed(jedge)%m2=inodepp
-         ed(jedge)%t2=ie
-         goto 111
+        jedge=nodeedgenumber(j,inodep)
+        ed(jedge)%m2=inodepp
+        ed(jedge)%t2=ie
+        goto 111
       endif
-   enddo
-   iedge=iedge+1
-   if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2$')
-   nedgepernode(inode)=nedgepernode(inode)+1
-   if (nedgepernode(inode).gt.nnmax) then
+    enddo
+    iedge=iedge+1
+    if (iedge.gt.nedge) call stop_run ('too many edges in delaunay2$')
+    nedgepernode(inode)=nedgepernode(inode)+1
+    if (nedgepernode(inode).gt.nnmax) then
       if (iproc.eq.0) write (8,*) 'node',inode,'has',nedgepernode(inode),'neighbours'
       call stop_run ('too many neighbours$')
-   endif
-   nodenodenumber(nedgepernode(inode),inode)=inodep
-   nodeedgenumber(nedgepernode(inode),inode)=iedge
-   ed(iedge)%n1=inode
-   ed(iedge)%n2=inodep
-   ed(iedge)%m1=inodepp
-   ed(iedge)%t1=ie
-   111 continue
-   enddo
+    endif
+    nodenodenumber(nedgepernode(inode),inode)=inodep
+    nodeedgenumber(nedgepernode(inode),inode)=iedge
+    ed(iedge)%n1=inode
+    ed(iedge)%n2=inodep
+    ed(iedge)%m1=inodepp
+    ed(iedge)%t1=ie
+    111 continue
+  enddo
 enddo
 
 if (iedge/=nedge) call stop_run ('pb_b in refine_surface$')
@@ -163,7 +163,7 @@ if (iedge/=nedge) call stop_run ('pb_b in refine_surface$')
 
 !distmax=1.d0/2.d0**surface%levelt*sqrt(2.d0)*surface%stretch
 distmax=1.d0/(2.d0**surface%levelt+1)*sqrt(2.d0)*surface%stretch   !opla
-distmin=1.d0/(2.d0**surface%levelt+4)*sqrt(2.d0)
+distmin=1.d0/(2.d0**(surface%levelt+4)+1)*sqrt(2.d0)
 
 allocate (refine(nedge),stat=err)
 allocate (refinep(nedge),stat=err)
@@ -183,13 +183,13 @@ do iedge=1+iproc,nedge,nproc
    prod=surface%xn(i1)*surface%xn(i2)+surface%yn(i1)*surface%yn(i2)+surface%zn(i1)*surface%zn(i2)
    !if (dist1.gt.distmax .or. prod.lt.prodmin) refinep(iedge)=1
    refinep(iedge)=(dist1.gt.distmax) 
-   removep(iedge)=(dist1.le.distmin)
+   removep(iedge)=(dist1.gt.eps .and. dist1.le.distmin)
+   !if (dist1.le.distmin) print *,'dist1 (',dist1,') less than distmin'
 enddo
 
 call mpi_allreduce (refinep,refine,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)
 call mpi_allreduce (removep,remove,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)
 
-
 !-----------------------
 ! compute nadd and naddp
 !-----------------------
@@ -213,10 +213,20 @@ do iedge=1,nedge
   if (remove(iedge)) then
      irem=irem+1
      if (ed(iedge)%t1.eq.0 .or. ed(iedge)%t2.eq.0) nremp=nremp+1
-     remove_list(irem)=iedge
+     !remove_list(irem)=iedge
+     remove_list(irem)=ed(iedge)%n1
   endif
 enddo
 
+! Modify remove_list to account for removal of nodes from surface
+if (nrem.gt.0) then
+  do i=1,nrem
+    do j=i,nrem
+      if (remove_list(i).lt.remove_list(j)) remove_list(j)=remove_list(j)-1
+    enddo
+  enddo
+endif
+
 if (iadd/=nadd) call stop_run ('pb_b in refine_surface$')
 if (irem/=nrem) call stop_run ('pb_b in refine_surface$')
 
@@ -308,6 +318,7 @@ deallocate (iconswap)
 nnodeint=surface%nsurface
 nelemint=surface%nt
 nedgeint=nedge
+
 do i=1,nadd
    iedge=refine_list(i)
    i1=ed(iedge)%n1
@@ -332,10 +343,24 @@ enddo
 
 if (iproc.eq.0 .and. params%debug>=1) write(*,'(a,i2,a,i4,a)') shift//'S.',is,':', nrem,' removed ptcls in refine_surface'
 
+do i=1,size(ed%n1)
+   i1=ed(i)%n1
+   i2=ed(i)%n2
+   if (i1.eq.i2) then
+     print *,'*******************************************************************'
+     print *,'Nodes for edge ',i,' have the same spatial coordinates'
+     print *,'The edges nodes are are node i1: ',i1,' and i2: ',i2
+     print *,'x1:',surface%x(i1),' y1: ',surface%y(i1),' z1: ',surface%z(i1)
+     print *,'x2:',surface%x(i2),' y2: ',surface%y(i2),' z2: ',surface%z(i2)
+     print *,'*******************************************************************'
+   endif
+enddo
+
+surface%nsurface=nnodemax
+surface%nt=nelemmax
+
 do i=1,nrem
-   iedge=remove_list(i)
-   i1=ed(iedge)%n1
-   i2=ed(iedge)%n2
+   i1=remove_list(i)
    call remove_point(i1,surface,0,surface0)
 end do
 
@@ -347,6 +372,12 @@ end if
 if (nnodeint.ne.nnodemax) call stop_run ('Error counting nodes in refine_surface$')
 if (nelemint.ne.nelemmax) call stop_run ('Error counting elements in refine_surface$')
 
+if (nrem.gt.0) then
+  nnodemax=surface%nsurface
+  nelemmax=surface%nt
+  nedgemax=nnodemax+nelemmax-1
+endif
+
 surface%nsurface=nnodemax
 surface%nt=nelemmax
 nedge=nedgemax
@@ -359,7 +390,6 @@ call heap (threadinfo,'nedgepernode','refine_surface',size(nedgepernode),'int',-
 call heap (threadinfo,'nodenodenumber','refine_surface',size(nodenodenumber),'int',-1) ; deallocate (nodenodenumber)
 call heap (threadinfo,'nodeedgenumber','refine_surface',size(nodeedgenumber),'int',-1) ; deallocate (nodeedgenumber)
 
-
 if (params%debug .ge.2) call output_surf(surface,is,'after_refine',istep,ref_count) 
 
 end subroutine