diff --git a/src/refine_surface.f90 b/src/refine_surface.f90
index d50bfe90dd100ca4ef9fe5d47e9995549424a531..53163e5f84a6138d2fd40d367c1157d38c1d9fec 100644
--- a/src/refine_surface.f90
+++ b/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