Added noise to edge length calculation for surface particle injection

src/refine_surface.f90

@@ -91,7 +91,7 @@ integer err,ierr,iproc,nproc,i,iedge,i1,i2
 integer nnodeint,nelemint,nedgeint,nelemmax,nnodemax,nedgemax
 double precision,dimension(:,:),allocatable::memswap
 integer,dimension(:,:),allocatable::iconswap
-double precision xm,ym,zm,xyzn
+double precision xm,ym,zm,xyzn,ran
 
 !------------------------------------------------------------------------------|
 !------------------------------------------------------------------------------|
@@ -189,9 +189,11 @@ do iedge=1+iproc,nedge,nproc
    !          surface%z(i1),surface%z(i2),params%distance_exponent)
    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) 
+   ! Add noise to point injection
+   call random_number(ran)
+   refinep(iedge)=(dist1.gt.distmax*(ran+2.5d0)/3.d0)
    if (params%remove_surf_pts) removep(iedge)=(dist1.gt.eps .and. dist1.le.distmin)
-   !if (dist1.le.distmin) print *,'dist1 (',dist1,') less than distmin'
+   !if (dist1.le.distmin) write(*,*) 'dist1 (',dist1,') less than distmin'
 enddo
 
 call mpi_allreduce (refinep,refine,nedge,mpi_logical,mpi_lor,mpi_comm_world,ierr)
@@ -338,6 +340,14 @@ do i=1,nadd
    iedge=refine_list(i)
    i1=ed(iedge)%n1
    i2=ed(iedge)%n2
+   if (i1 == i2 .and. iproc == 0) then
+     write(*,*) '*******************************************************************'
+     write(*,*) 'Nodes for edge ',i,' have the same spatial coordinates before refinement!'
+     write(*,*) 'The edges nodes are are node i1: ',i1,' and i2: ',i2
+     write(*,*) 'x1:',surface%x(i1),' y1: ',surface%y(i1),' z1: ',surface%z(i1)
+     write(*,*) 'x2:',surface%x(i2),' y2: ',surface%y(i2),' z2: ',surface%z(i2)
+     write(*,*) '*******************************************************************'
+   endif
    call middle (xm,ym,zm,i1,i2,surface%x,surface%y,surface%z, &
                 surface%xn,surface%yn,surface%zn,nnodeint)
    nnodeint=nnodeint+1
@@ -362,12 +372,12 @@ do i=1,size(ed%n1)
    i1=ed(i)%n1
    i2=ed(i)%n2
    if (i1 == i2 .and. iproc == 0) 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 *,'*******************************************************************'
+     write(*,*) '*******************************************************************'
+     write(*,*) 'Nodes for edge ',i,' have the same spatial coordinates'
+     write(*,*) 'The edges nodes are are node i1: ',i1,' and i2: ',i2
+     write(*,*) 'x1:',surface%x(i1),' y1: ',surface%y(i1),' z1: ',surface%z(i1)
+     write(*,*) 'x2:',surface%x(i2),' y2: ',surface%y(i2),' z2: ',surface%z(i2)
+     write(*,*) '*******************************************************************'
    endif
 enddo