module optvor use iso_c_binding, only: c_ptr,c_f_pointer use VorCells_Polygons, only: Polygon,VoronoiCell,voronoi_cell,voronoi_cell_old,VorCellInterConvPol, & VorCellInterConvPol_Collinear,voronoi_cell_destroy,polygon_destroy implicit none type :: pdata_type integer :: npols real(kind=8) :: c(2),rmax,volA,xmin,xmax,ymin,ymax,drawscale,J0delta,J2c(2),J2r,J3ctol,J4ctol integer, allocatable :: nvpols(:),poledges(:,:) real(kind=8), allocatable :: xpol(:,:),ypol(:,:) real(kind=8) :: Jweight(0:5) ! 0: minimum distance between sites ! 1: cells with identical volume given by vol(A)/m ! 2: cells with volume proportional to density function f evaluated at their sites ! 3: each cell with equally sized edges ! 4: each cell with equally sized angles ! 5: each cell with edges being crossed at their middle point by their corresponding Delaunay edges end type pdata_type public :: evalfg,evalfc,evalgjac,compJ,proj,approxg,drawsol,InA,compscale contains ! ****************************************************************** ! ****************************************************************** logical function InA(p,pdataptr) implicit none ! SCALAR ARGUMENTS type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: p(2) ! LOCAL SCALARS integer :: allocerr,j,maxnv type(pdata_type), pointer :: pdata ! LOCAL ARRAYS real(kind=8), allocatable :: vertices(:,:) call c_f_pointer(pdataptr,pdata) maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(vertices(2,maxnv),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'InA: Allocation error.' stop end if j = 1 InA = .false. do while ( j .le. pdata%npols .and. .not. InA ) vertices(1,1:pdata%nvpols(j)-1) = pdata%xpol(j,1:pdata%nvpols(j)-1) vertices(2,1:pdata%nvpols(j)-1) = pdata%ypol(j,1:pdata%nvpols(j)-1) call polygon_contains_point_2d_convex(pdata%nvpols(j)-1,vertices,p,inA) j = j + 1 end do deallocate(vertices,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'InA: Deallocation error.' stop end if end function InA ! ****************************************************************** ! ****************************************************************** subroutine comparea(scale,area,pdataptr) implicit none ! SCALAR ARGUMENTS real(kind=8), intent(in) :: scale real(kind=8), intent(out) :: area type(c_ptr), optional, intent(in) :: pdataptr ! LOCAL SCALARS type(pdata_type), pointer :: pdata ! LOCAL SCALARS integer :: allocerr,j,maxnv real(kind=8) :: areaPj ! LOCAL ARRAYS real(kind=8), allocatable :: poltmp(:,:) call c_f_pointer(pdataptr,pdata) maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(poltmp(2,maxnv),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'Allocation error.' stop end if area = 0.0d0 do j = 1,pdata%npols poltmp(1,1:pdata%nvpols(j)) = scale * pdata%xpol(j,1:pdata%nvpols(j)) poltmp(2,1:pdata%nvpols(j)) = scale * pdata%ypol(j,1:pdata%nvpols(j)) call polygon_area_2d(pdata%nvpols(j)-1,poltmp(1:2,1:pdata%nvpols(j)-1),areaPj) area = area + areaPj end do deallocate(poltmp,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'Deallocation error.' stop end if end subroutine comparea ! ****************************************************************** ! ****************************************************************** subroutine compscale(nsites,scale,area,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: nsites real(kind=8), intent(out) :: area,scale type(c_ptr), optional, intent(in) :: pdataptr ! PARAMETERS real(kind=8), parameter :: eps = 0.01d0 ! LOCAL SCALARS real(kind=8) :: a,b,phia,phib,tphi,tscale,tarea type(pdata_type), pointer :: pdata call c_f_pointer(pdataptr,pdata) tscale = 1.0d0 call comparea(tscale,tarea,pdataptr) tphi = tarea / nsites - 1.0d0 a = tscale phia = tphi write(*,*) 'scale factor (a) = ',tscale,' ideal area (a) = ',tarea / nsites if ( abs( tphi ) .le. eps ) then scale = tscale area = tarea return end if if ( tphi .lt. 0.0d0 ) then do while ( tphi .lt. 0.0d0 ) tscale = 2.0d0 * tscale call comparea(tscale,tarea,pdataptr) tphi = tarea / nsites - 1.0d0 end do b = tscale phib = tphi else do while ( tphi .gt. 0.0d0 ) tscale = 0.5d0 * tscale call comparea(tscale,tarea,pdataptr) tphi = tarea / nsites - 1.0d0 end do b = tscale phib = tphi end if write(*,*) 'scale factor (b) = ',tscale,' ideal area (b) = ',tarea / nsites if ( abs( tphi ) .le. eps ) then scale = tscale area = tarea return end if 10 continue tscale = 0.5d0 * ( a + b ) call comparea(tscale,tarea,pdataptr) tphi = tarea / nsites - 1.0d0 write(*,*) 'scale factor (c) = ',tscale,' ideal area (c) = ',tarea / nsites if ( abs( tphi ) .le. eps ) then scale = tscale area = tarea return end if if ( sign(1.0d0,tphi) .eq. sign(1.0d0,phia) ) then a = tscale phia = tphi else b = tscale phib = tphi end if go to 10 end subroutine compscale ! ****************************************************************** ! ****************************************************************** subroutine proj(n,x,flag,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n integer, intent(out) :: flag type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(inout) :: x(n) ! LOCAL SCALARS integer :: i,nsites real(kind=8) :: tmp type(pdata_type), pointer :: pdata call c_f_pointer(pdataptr,pdata) flag = 0 nsites = n / 2 do i = 1,nsites tmp = norm2( x(2*i-1:2*i) - pdata%c(1:2) ) if ( tmp .gt. pdata%rmax ) then x(2*i-1:2*i) = pdata%c(1:2) + pdata%rmax * ( x(2*i-1:2*i) - pdata%c(1:2) ) / tmp end if end do end subroutine proj ! ****************************************************************** ! ****************************************************************** subroutine approxg(n,x,g,h,flag,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n integer, intent(out) :: flag real(kind=8), intent(in) :: h type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(inout) :: x(n) real(kind=8), intent(out) :: g(n) ! LOCAL SCALARS integer :: i real(kind=8) :: fplus,fminus,tmp ! LOCAL ARRAYS real(kind=8) :: gtmp(n) do i = 1,n tmp = x(i) x(i) = tmp + h call evalfg(n,x,fplus,gtmp,flag,pdataptr) x(i) = tmp - h call evalfg(n,x,fminus,gtmp,flag,pdataptr) g(i) = ( fplus - fminus ) / ( 2.0d0 * h ) x(i) = tmp end do end subroutine approxg ! ****************************************************************** ! ****************************************************************** subroutine evalfg(n,x,f,g,flag,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n integer, intent(out) :: flag real(kind=8), intent(out) :: f type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: g(n) ! PARAMETERS logical, parameter :: compgrad = .true., debug = .false. integer, parameter :: gtmpnnzmax = 100000 ! LOCAL SCALARS logical :: wempty integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P ! LOCAL ARRAYS integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(gtmpnnzmax),sites(2,n/2),vorxy(2,2*(n/2)) type(Polygon), allocatable :: W(:) flag = 0 call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .lt. 3 ) then write(*,*) 'In evalfg: nsites must be at least 3!' flag = - 1 return end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror flag = - 1 return end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In evalfg: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' flag = - 1 return end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In evalfg: triangle_circumcenter_2d returned NaN.' flag = - 1 return end if end do ! Compute density function at sites if required if ( pdata%Jweight(2) .gt. 0.0d0 ) then do i = 1,nsites call fgcompJ2(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr) end do end if f = 0.0d0 if ( compgrad ) g(1:n) = 0.0d0 maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(W(pdata%npols),P%v(2,maxnv),P%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Allocation error.' flag = - 1 return end if do i = 1,nsites if ( debug ) then write(*,*) write(*,*) 'SITE i = ',i write(*,*) 'Generator: ',sites(1:2,i) write(*,*) end if call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status) if ( status .ne. 0 ) then write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status flag = - 1 return end if do j = 1,pdata%npols P%n = pdata%nvpols(j) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(j,1:P%n) P%v(2,1:P%n) = pdata%ypol(j,1:P%n) P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W(j)) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j)) end if end do if ( pdata%Jweight(0) .gt. 0.0d0 ) then call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad) if ( flag .ne. 0 ) return f = f + pdata%Jweight(0) * ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(0) * gtmpval(1:gtmpnnz) end if if ( pdata%Jweight(1) .gt. 0.0d0 ) then call compJ1(nsites,sites,i,pdata%npols,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + pdata%Jweight(1) * ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(1) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz) end if if ( pdata%Jweight(2) .gt. 0.0d0 ) then call compJ2(nsites,sites,i,pdata%npols,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + pdata%Jweight(2) * ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(2) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz) end if wempty = .true. do j = 1,pdata%npols if ( .not. W(j)%deg ) then wempty = .false. end if end do if ( wempty ) then f = f + 1.0d0 else do j = 1,pdata%npols if ( .not. W(j)%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W(j)%n - 1 write(*,*) end if if ( pdata%Jweight(3) .gt. 0.0d0 ) then call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + pdata%Jweight(3) * ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(3) * gtmpval(1:gtmpnnz) end if if ( pdata%Jweight(4) .gt. 0.0d0 ) then call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + pdata%Jweight(4) * ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(4) * gtmpval(1:gtmpnnz) end if if ( pdata%Jweight(5) .gt. 0.0d0 ) then call compJ5(nsites,sites,i,W(j),ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + pdata%Jweight(5) * ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(5) * gtmpval(1:gtmpnnz) end if else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.' write(*,*) end if end if end do end if call voronoi_cell_destroy(V) do j = 1,pdata%npols call polygon_destroy(W(j)) end do end do f = f / nsites if ( compgrad ) g(1:n) = g(1:n) / nsites deallocate(W,P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine evalfg ! ****************************************************************** ! ****************************************************************** subroutine evalfc(n,x,f,m,c,flag,debug,pdataptr) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: debug integer, intent(in) :: m,n integer, intent(out) :: flag real(kind=8), intent(out) :: f type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: c(m) ! PARAMETER integer, parameter :: lim = 1 logical, parameter :: compgrad = .false. ! LOCAL SCALARS logical :: lmem integer :: jcnnz ! LOCAL ARRAYS integer :: jcfun(lim),jcvar(lim) real(kind=8) :: g(1),jcval(lim) call evalfcgjac(n,x,f,g,m,c,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,.false.,debug,pdataptr) end subroutine evalfc ! ****************************************************************** ! ****************************************************************** subroutine evalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,debug,pdataptr) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: debug logical, intent(out) :: lmem integer, intent(in) :: lim,m,n integer, intent(out) :: flag,jcnnz type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS integer, intent(out) :: jcfun(lim),jcvar(lim) real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: g(n),jcval(lim) ! LOCAL SCALARS real(kind=8) :: f ! LOCAL ARRAYS real(kind=8) :: c(m) call evalfcgjac(n,x,f,g,m,c,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,.true.,debug,pdataptr) end subroutine evalgjac ! ****************************************************************** ! ****************************************************************** subroutine evalfcgjac(n,x,f,g,m,c,jcfun,jcvar,jcval,jcnnz,lim,lmem,flag,compgrad,debug,pdataptr) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug logical, intent(out) :: lmem integer, intent(in) :: lim,m,n integer, intent(out) :: flag,jcnnz real(kind=8), intent(out) :: f type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS integer, intent(out) :: jcfun(lim),jcvar(lim) real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: c(m),g(n),jcval(lim) ! LOCAL SCALARS logical :: wempty integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P ! LOCAL ARRAYS integer :: gtmpind(max(n,lim)),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(max(n,lim)),sites(2,n/2),vorxy(2,2*(n/2)) type(Polygon), allocatable :: W(:) flag = 0 call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .lt. 3 ) then write(*,*) 'In evalfg: nsites must be at least 3!' flag = - 1 return end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In evalfcgjac: dtris2 returned an error different from 0 and 225: ',ierror flag = - 1 return end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In evalfcgjac: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' flag = - 1 return end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In evalfcgjac: triangle_circumcenter_2d returned NaN.' flag = - 1 return end if end do ! Compute density function at sites if required do i = 1,nsites call fgcompJ2(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr) end do f = 0.0d0 if ( compgrad ) then lmem = .false. g(1:n) = 0.0d0 jcnnz = 0 end if maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(W(pdata%npols),P%v(2,maxnv),P%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfcgjac: Allocation error.' flag = - 1 return end if do i = 1,nsites if ( debug ) then write(*,*) write(*,*) 'SITE i = ',i write(*,*) 'Generator: ',sites(1:2,i) write(*,*) end if call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status) ! call voronoi_cell_old(i,nsites,sites,ntriag,til,tnbr,vorxy,V,status) if ( status .ne. 0 ) then write(*,*) 'In evalfcgjac: Error in voronoi_cell, status = ',status flag = - 1 return end if do j = 1,pdata%npols P%n = pdata%nvpols(j) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(j,1:P%n) P%v(2,1:P%n) = pdata%ypol(j,1:P%n) P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W(j)) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j)) end if end do call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad) if ( flag .ne. 0 ) return f = f + 10.0d0 * ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + 10.0d0 * gtmpval(1:gtmpnnz) call compJ1(nsites,sites,i,pdata%npols,W,pdata%volA,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + 2.0d0 * ftmp * gtmpval(1:gtmpnnz) wempty = .true. do j = 1,pdata%npols if ( .not. W(j)%deg ) then wempty = .false. end if end do if ( wempty ) then c(i) = 1.0d0 c(nsites+i) = 1.0d0 else do j = 1,pdata%npols if ( .not. W(j)%deg ) then call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return !write(*,*) 'For i = ',i,' Constraint value = ',ftmp c(i) = ftmp if ( compgrad ) then if ( jcnnz + gtmpnnz .gt. lim ) then lmem = .true. return end if if ( any(isnan(gtmpval(1:gtmpnnz))) ) then write(*,*) 'There is a NaN here!' call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,.true.) if ( flag .ne. 0 ) return stop end if jcfun(jcnnz+1:jcnnz+gtmpnnz) = i jcvar(jcnnz+1:jcnnz+gtmpnnz) = gtmpind(1:gtmpnnz) jcval(jcnnz+1:jcnnz+gtmpnnz) = gtmpval(1:gtmpnnz) jcnnz = jcnnz + gtmpnnz end if call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return c(nsites+i) = ftmp if ( compgrad ) then if ( jcnnz + gtmpnnz .gt. lim ) then lmem = .true. return end if if ( any(isnan(gtmpval(1:gtmpnnz))) ) then write(*,*) 'There is a NaN here!' call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,max(n,lim),gtmpnnz,gtmpind,gtmpval,flag,compgrad,.true.) if ( flag .ne. 0 ) return stop end if jcfun(jcnnz+1:jcnnz+gtmpnnz) = nsites + i jcvar(jcnnz+1:jcnnz+gtmpnnz) = gtmpind(1:gtmpnnz) jcval(jcnnz+1:jcnnz+gtmpnnz) = gtmpval(1:gtmpnnz) jcnnz = jcnnz + gtmpnnz end if end if end do end if call voronoi_cell_destroy(V) do j = 1,pdata%npols call polygon_destroy(W(j)) end do end do f = f / nsites if (compgrad ) g(1:n) = g(1:n) / nsites deallocate(W,P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine evalfcgjac ! ****************************************************************** ! ****************************************************************** subroutine compJ(n,x,f,g,flag,compgrad,debug,pdataptr) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: n integer, intent(out) :: flag real(kind=8), intent(out) :: f(0:5) type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: g(n,0:5) ! PARAMETERS integer, parameter :: gtmpnnzmax = 100000 ! LOCAL SCALARS logical :: wempty integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P ! LOCAL ARRAYS integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(gtmpnnzmax),sites(2,n/2),vorxy(2,2*(n/2)) type(Polygon), allocatable :: W(:) flag = 0 call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .lt. 3 ) then write(*,*) 'In evalfg: nsites must be at least 3!' flag = - 1 return end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror flag = - 1 return end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In evalfg: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' flag = - 1 return end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In evalfg: triangle_circumcenter_2d returned NaN.' flag = - 1 return end if end do ! Compute density function at sites if required do i = 1,nsites call fgcompJ2(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr) end do f(0:5) = 0.0d0 if ( compgrad ) g(1:n,0:5) = 0.0d0 maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(W(pdata%npols),P%v(2,maxnv),P%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Allocation error.' flag = - 1 return end if do i = 1,nsites if ( debug ) then write(*,*) write(*,*) 'SITE i = ',i write(*,*) 'Generator: ',sites(1:2,i) write(*,*) end if call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status) ! call voronoi_cell_old(i,nsites,sites,ntriag,til,tnbr,vorxy,V,status) if ( status .ne. 0 ) then write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status flag = - 1 return end if do j = 1,pdata%npols P%n = pdata%nvpols(j) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(j,1:P%n) P%v(2,1:P%n) = pdata%ypol(j,1:P%n) P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W(j)) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j)) end if end do call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad) if ( flag .ne. 0 ) return f(0) = f(0) + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz),0) = g(gtmpind(1:gtmpnnz),0) + gtmpval(1:gtmpnnz) call compJ1(nsites,sites,i,pdata%npols,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(1) = f(1) + ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz),1) = g(gtmpind(1:gtmpnnz),1) + 2.0d0 * ftmp * gtmpval(1:gtmpnnz) call compJ2(nsites,sites,i,pdata%npols,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(2) = f(2) + ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz),2) = g(gtmpind(1:gtmpnnz),2) + 2.0d0 * ftmp * gtmpval(1:gtmpnnz) wempty = .true. do j = 1,pdata%npols if ( .not. W(j)%deg ) then wempty = .false. end if end do if ( wempty ) then f = f + 1.0d0 else do j = 1,pdata%npols if ( .not. W(j)%deg ) then call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(3) = f(3) + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz),3) = g(gtmpind(1:gtmpnnz),3) + gtmpval(1:gtmpnnz) call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(4) = f(4) + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz),4) = g(gtmpind(1:gtmpnnz),4) + gtmpval(1:gtmpnnz) call compJ5(nsites,sites,i,W(j),ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(5) = f(5) + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz),5) = g(gtmpind(1:gtmpnnz),5) + gtmpval(1:gtmpnnz) end if end do end if call voronoi_cell_destroy(V) do j = 1,pdata%npols call polygon_destroy(W(j)) end do end do f(0:5) = f(0:5) / nsites if ( compgrad ) g(1:n,0:5) = g(1:n,0:5) / nsites deallocate(W,P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine compJ ! ****************************************************************** ! ****************************************************************** subroutine compallJ(n,x,f,debug,flag,pdataptr) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: debug integer, intent(in) :: n integer, intent(out) :: flag real(kind=8), intent(out) :: f(0:5,n/2) type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) ! PARAMETERS logical, parameter :: compgrad = .false. integer, parameter :: gtmpnnzmax = 1 ! LOCAL SCALARS logical :: wempty integer :: allocerr,gtmpnnz,i,ierror,j,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P ! LOCAL ARRAYS integer :: gtmpind(gtmpnnzmax),indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: fden(n/2),gden(2,n/2),gtmpval(gtmpnnzmax),sites(2,n/2),vorxy(2,2*(n/2)) type(Polygon), allocatable :: W(:) flag = 0 call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .lt. 3 ) then write(*,*) 'In evalfg: nsites must be at least 3!' flag = - 1 return end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In evalfg: dtris2 returned an error different from 0 and 225: ',ierror flag = - 1 return end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In evalfg: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' flag = - 1 return end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In evalfg: triangle_circumcenter_2d returned NaN.' flag = - 1 return end if end do ! Compute density function at sites if required do i = 1,nsites call fgcompJ2(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr) end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(W(pdata%npols),P%v(2,maxnv),P%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Allocation error.' flag = - 1 return end if do i = 1,nsites if ( debug ) then write(*,*) write(*,*) 'SITE i = ',i write(*,*) 'Generator: ',sites(1:2,i) write(*,*) end if call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,V,status) ! call voronoi_cell_old(i,nsites,sites,ntriag,til,tnbr,vorxy,V,status) if ( status .ne. 0 ) then write(*,*) 'In evalfg: Error in voronoi_cell, status = ',status flag = - 1 return end if do j = 1,pdata%npols P%n = pdata%nvpols(j) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(j,1:P%n) P%v(2,1:P%n) = pdata%ypol(j,1:P%n) P%e(1:P%n-1) = pdata%poledges(j,1:P%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W(j)) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W(j)) end if end do call compJ0(nsites,sites,i,pdata%npols,W,pdata%J0delta,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad) if ( flag .ne. 0 ) return f(0,i) = ftmp call compJ1(nsites,sites,i,pdata%npols,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(1,i) = ftmp call compJ2(nsites,sites,i,pdata%npols,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(2,i) = ftmp wempty = .true. do j = 1,pdata%npols if ( .not. W(j)%deg ) then wempty = .false. end if end do if ( wempty ) then f(3,i) = 1.0d0 f(4,i) = 1.0d0 f(5,i) = 1.0d0 else do j = 1,pdata%npols if ( .not. W(j)%deg ) then call compJ3(nsites,sites,i,W(j),pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(3,i) = ftmp call compJ4(nsites,sites,i,W(j),pdata%J4ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(4,i) = ftmp call compJ5(nsites,sites,i,W(j),ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(5,i) = ftmp end if end do end if call voronoi_cell_destroy(V) do j = 1,pdata%npols call polygon_destroy(W(j)) end do end do deallocate(W,P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine compallJ ! ****************************************************************** ! ****************************************************************** subroutine compJ0(nsites,sites,i,npols,W,delta,f,gnnzmax,gnnz,gind,gval,flag,compgrad) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad integer, intent(in) :: gnnzmax,i,npols,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: delta real(kind=8), intent(out) :: f ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) type(Polygon), intent(in) :: W(npols) ! LOCAL SCALARS integer :: ed,j real(kind=8) :: dist2 ! LOCAL ARRAYS real(kind=8) :: diff(2) f = 0.0d0 if ( compgrad ) then gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ0: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gnnz = 2 gind(1:2) = (/ 2*i-1, 2*i /) gval(1:2) = 0.0d0 end if do j = 1,npols if ( .not. W(j)%deg ) then do ed = 1,W(j)%n - 1 if ( W(j)%e(ed) .lt. 0 ) then diff(1:2) = sites(1:2,i) - sites(1:2,abs(W(j)%e(ed))) dist2 = sum( diff(1:2) ** 2 ) if ( dist2 .lt. delta ** 2 ) then f = f + ( delta ** 2 - dist2 ) ** 2 if ( compgrad ) then gval(1:2) = gval(1:2) - 4.0d0 * ( delta ** 2 - dist2 ) * diff(1:2) if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ0: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = 4.0d0 * ( delta ** 2 - dist2 ) * diff(1:2) gnnz = gnnz + 2 end if end if end if end do end if end do end subroutine compJ0 ! ****************************************************************** ! ****************************************************************** subroutine compJ1alternative(nsites,sites,i,npols,W,volA,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,npols,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: volA real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W(npols) ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,edprev,ednext,j real(kind=8) :: mu,phia,phib,phic,vol,volj,xi ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),normal(2) vol = 0.0d0 do j = 1,npols if ( .not. W(j)%deg ) then !!$ volj = 0.0d0 !!$ do ed = 1,W(j)%n - 1 !!$ volj = volj + 0.5d0 * ( W(j)%v(1,ed) * W(j)%v(2,ed+1) - W(j)%v(1,ed+1) * W(j)%v(2,ed) ) !!$ end do call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj) vol = vol + volj end if end do xi = nsites / volA if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = xi * vol - 1.0d0 if ( .not. compgrad ) then return end if gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /) gval(gnnz+1:gnnz+2) = 0.0d0 gnnz = gnnz + 2 do j = 1,npols if ( .not. W(j)%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W(j)%n - 1 write(*,*) end if do ed = 1,W(j)%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W(j)%v(1:2,ed) write(*,*) '====================' write(*,*) end if ! Processing edge (v_E,w_E) if ( W(j)%e(ed) .lt. 0 ) then if ( ed .eq. 1 ) then edprev = W(j)%n - 1 ednext = ed + 1 else if ( ed .eq. W(j)%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if if ( debug ) then write(*,*) 'Previous: ',edprev,' next: ',ednext write(*,*) 'They correspond to the points ',W(j)%v(1:2,edprev),' and ',W(j)%v(1:2,ednext),' respectively.' end if ! Computing the outward unit normal vector to the edge ed call line_exp_normal_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ed),normal) if ( debug ) then write(*,*) 'Outward unit normal vector to the edge ed: ',normal(1:2) end if mu = 0.5d0 * norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) ) if ( debug ) then write(*,*) 'mu = ',mu end if ! Starting point v_E if ( debug ) then write(*,*) write(*,*) 'Processing v_E' write(*,*) end if if ( W(j)%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=3. ', & 'Neighbour cells are: ',abs(W(j)%e(edprev)),' and ',abs(W(j)%e(ed)) end if call matrixint(sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(edprev))-1, 2*abs(W(j)%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed)) end if call line_exp2imp_2d(W(j)%v(1:2,ed),W(j)%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 end if ! Ending point w_E if ( debug ) then write(*,*) write(*,*) 'Processing w_E' write(*,*) end if if ( W(j)%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W(j)%e(ed)),' and ',abs(W(j)%e(ednext)) end if call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ednext))-1, 2*abs(W(j)%e(ednext)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed)) end if call line_exp2imp_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ednext+1),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 end if else if ( debug ) then write(*,*) 'This is a fixed edge.' write(*,*) 'So, its vertices v_E and w_E are ignored.' end if end if end do else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.' write(*,*) end if end if end do end subroutine compJ1alternative ! ****************************************************************** ! ****************************************************************** subroutine compJ1(nsites,sites,i,npols,W,volA,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,npols,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: volA real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W(npols) ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,ednext,j real(kind=8) :: mu,vol,volj,xi ! LOCAL ARRAYS real(kind=8) :: p(2) vol = 0.0d0 do j = 1,npols if ( .not. W(j)%deg ) then !!$ volj = 0.0d0 !!$ do ed = 1,W(j)%n - 1 !!$ volj = volj + 0.5d0 * ( W(j)%v(1,ed) * W(j)%v(2,ed+1) - W(j)%v(1,ed+1) * W(j)%v(2,ed) ) !!$ end do call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj) vol = vol + volj end if end do xi = nsites / volA if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = xi * vol - 1.0d0 if ( .not. compgrad ) then return end if if ( 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',2 flag = - 1 return end if gind(1:2) = (/ 2*i-1, 2*i /) gval(1:2) = 0.0d0 gnnz = 2 do j = 1,npols if ( .not. W(j)%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W(j)%n - 1 write(*,*) end if do ed = 1,W(j)%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W(j)%v(1:2,ed) write(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W(j)%v(1:2,ednext) write(*,*) '====================' write(*,*) end if if ( W(j)%e(ed) .lt. 0 ) then if ( ed .eq. W(j)%n - 1 ) then ednext = 1 else ednext = ed + 1 end if if ( debug ) then write(*,*) 'Next: ',ednext write(*,*) 'It corresponds to the point ',W(j)%v(1:2,ednext) end if mu = norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W(j)%e(ed))) ) if ( debug ) then write(*,*) 'mu = ',mu end if p(1:2) = 0.5d0 * ( W(j)%v(1:2,ed) + W(j)%v(1:2,ednext) ) gval(1:2) = gval(1:2) + xi * mu * ( p(1:2) - sites(1:2,i) ) if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ1: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W(j)%e(ed))) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'This is a fixed edge.' write(*,*) 'So, its vertices v_E and w_E are ignored.' end if end if end do else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.' write(*,*) end if end if end do end subroutine compJ1 ! ****************************************************************** ! ****************************************************************** subroutine compJ2(nsites,sites,i,npols,W,volA,fden,gden,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,npols,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: volA real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W(npols) ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: fden(nsites),gden(2,nsites),sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,ednext,j real(kind=8) :: mu,vol,volj,xi ! LOCAL ARRAYS real(kind=8) :: p(2) vol = 0.0d0 do j = 1,npols if ( .not. W(j)%deg ) then !!$ volj = 0.0d0 !!$ do ed = 1,W(j)%n - 1 !!$ volj = volj + 0.5d0 * ( W(j)%v(1,ed) * W(j)%v(2,ed+1) - W(j)%v(1,ed+1) * W(j)%v(2,ed) ) !!$ end do call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj) vol = vol + volj end if end do xi = nsites / volA if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if ! f = xi * vol - 1.0d0 f = xi * vol - fden(i) if ( .not. compgrad ) then return end if if ( 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',2 flag = - 1 return end if gind(1:2) = (/ 2*i-1, 2*i /) ! gval(1:2) = 0.0d0 gval(1:2) = - gden(1:2,i) gnnz = 2 do j = 1,npols if ( .not. W(j)%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W(j)%n - 1 write(*,*) end if do ed = 1,W(j)%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W(j)%v(1:2,ed) write(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W(j)%v(1:2,ednext) write(*,*) '====================' write(*,*) end if if ( W(j)%e(ed) .lt. 0 ) then if ( ed .eq. W(j)%n - 1 ) then ednext = 1 else ednext = ed + 1 end if if ( debug ) then write(*,*) 'Next: ',ednext write(*,*) 'It corresponds to the point ',W(j)%v(1:2,ednext) end if mu = norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W(j)%e(ed))) ) if ( debug ) then write(*,*) 'mu = ',mu end if p(1:2) = 0.5d0 * ( W(j)%v(1:2,ed) + W(j)%v(1:2,ednext) ) gval(1:2) = gval(1:2) + xi * mu * ( p(1:2) - sites(1:2,i) ) if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W(j)%e(ed))) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'This is a fixed edge.' write(*,*) 'So, its vertices v_E and w_E are ignored.' end if end if end do else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.' write(*,*) end if end if end do end subroutine compJ2 ! ****************************************************************** ! ****************************************************************** subroutine compJ2alternative1(nsites,sites,i,npols,W,volA,fden,gden,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,npols,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: volA real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W(npols) ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: fden(nsites),gden(2,nsites),sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,edprev,ednext,j real(kind=8) :: mu,phia,phib,phic,vol,volj,xi ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),normal(2) vol = 0.0d0 do j = 1,npols if ( .not. W(j)%deg ) then call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj) vol = vol + volj end if end do xi = nsites / volA if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if ! f = xi * vol - 1.0d0 f = xi * vol - fden(i) if ( .not. compgrad ) then return end if gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /) ! gval(gnnz+1:gnnz+2) = 0.0d0 gval(gnnz+1:gnnz+2) = - gden(1:2,i) gnnz = gnnz + 2 do j = 1,npols if ( .not. W(j)%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W(j)%n - 1 write(*,*) end if do ed = 1,W(j)%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W(j)%v(1:2,ed) write(*,*) '====================' write(*,*) end if ! Processing edge (v_E,w_E) if ( W(j)%e(ed) .lt. 0 ) then if ( ed .eq. 1 ) then edprev = W(j)%n - 1 ednext = ed + 1 else if ( ed .eq. W(j)%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if if ( debug ) then write(*,*) 'Previous: ',edprev,' next: ',ednext write(*,*) 'They correspond to the points ',W(j)%v(1:2,edprev),' and ',W(j)%v(1:2,ednext),' respectively.' end if ! Computing the outward unit normal vector to the edge ed call line_exp_normal_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ed),normal) if ( debug ) then write(*,*) 'Outward unit normal vector to the edge ed: ',normal(1:2) end if mu = 0.5d0 * norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) ) if ( debug ) then write(*,*) 'mu = ',mu end if ! Starting point v_E if ( debug ) then write(*,*) write(*,*) 'Processing v_E' write(*,*) end if if ( W(j)%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=3. ', & 'Neighbour cells are: ',abs(W(j)%e(edprev)),' and ',abs(W(j)%e(ed)) end if call matrixint(sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(edprev))-1, 2*abs(W(j)%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed)) end if call line_exp2imp_2d(W(j)%v(1:2,ed),W(j)%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 end if ! Ending point w_E if ( debug ) then write(*,*) write(*,*) 'Processing w_E' write(*,*) end if if ( W(j)%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W(j)%e(ed)),' and ',abs(W(j)%e(ednext)) end if call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ednext))-1, 2*abs(W(j)%e(ednext)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed)) end if call line_exp2imp_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ednext+1),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W(j)%e(ed))-1, 2*abs(W(j)%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) gnnz = gnnz + 2 end if else if ( debug ) then write(*,*) 'This is a fixed edge.' write(*,*) 'So, its vertices v_E and w_E are ignored.' end if end if end do else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.' write(*,*) end if end if end do end subroutine compJ2alternative1 ! ****************************************************************** ! ****************************************************************** subroutine compJ2alternative2(nsites,sites,i,npols,W,volA,fden,gden,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,npols,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: volA real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W(npols) ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: fden(nsites),gden(2,nsites),sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,edprev,ednext,ind,j,k real(kind=8) :: mu,phia,phib,phic,sumfden,vol,volj,xi ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),normal(2) vol = 0.0d0 do j = 1,npols if ( .not. W(j)%deg ) then call polygon_area_2d(W(j)%n-1,W(j)%v(1:2,1:W(j)%n-1),volj) vol = vol + volj end if end do ! The code below corresponds to the formulation ! ! ( vol - desired_vol )^2; ! ! while the one that follows corresponds to the formulation ! ! ( vol / desired_vol - 1 )^2. !!$ sumfden = sum( fden(1:nsites) ) !!$ !!$ f = f + ( vol - fden(i) * volA / sumfden ) ** 2 !!$ !!$ xi = 2.0d0 * ( vol - fden(i) * volA / sumfden ) !!$ !!$ ind = 2 * i - 1 !!$ g(ind:ind+1) = g(ind:ind+1) - ( xi * volA / sumfden ) * gden(1:2,i) !!$ !!$ do k = 1,nsites !!$ ind = 2 * k - 1 !!$ g(ind:ind+1) = g(ind:ind+1) + ( xi * volA * fden(i) / sumfden ** 2 ) * gden(1:2,k) !!$ end do sumfden = sum( fden(1:nsites) ) xi = sumfden / ( fden(i) * volA ) if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = xi * vol - 1.0d0 if ( .not. compgrad ) then return end if if ( 2 * nsites .gt. gnnzmax ) then write(*,*) 'In compJ2: Increase gnnzmax to at least ',2 * nsites flag = - 1 return end if gnnz = 2 * nsites gind(1:gnnz) = (/ (ind,ind=1,gnnz) /) gval(1:gnnz) = 0.0d0 ind = 2 * i - 1 gval(ind:ind+1) = gval(ind:ind+1) - ( sumfden / ( fden(i) ** 2 * volA ) ) * vol * gden(1:2,i) do k = 1,nsites ind = 2 * k - 1 gval(ind:ind+1) = gval(ind:ind+1) + ( 1.0d0 / ( fden(i) * volA ) ) * vol * gden(1:2,k) end do do j = 1,npols if ( .not. W(j)%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W(j)%n - 1 write(*,*) end if do ed = 1,W(j)%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W(j)%v(1:2,ed) write(*,*) '====================' write(*,*) end if ! Processing edge (v_E,w_E) if ( W(j)%e(ed) .lt. 0 ) then if ( ed .eq. 1 ) then edprev = W(j)%n - 1 ednext = ed + 1 else if ( ed .eq. W(j)%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if if ( debug ) then write(*,*) 'Previous: ',edprev,' next: ',ednext write(*,*) 'They correspond to the points ',W(j)%v(1:2,edprev),' and ',W(j)%v(1:2,ednext),' respectively.' end if ! Computing the outward unit normal vector to the edge ed call line_exp_normal_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ed),normal) if ( debug ) then write(*,*) 'Outward unit normal vector to the edge ed: ',normal(1:2) end if mu = 0.5d0 * norm2( W(j)%v(1:2,ed) - W(j)%v(1:2,ednext) ) if ( debug ) then write(*,*) 'mu = ',mu end if ! Starting point v_E if ( debug ) then write(*,*) write(*,*) 'Processing v_E' write(*,*) end if if ( W(j)%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=3. ', & 'Neighbour cells are: ',abs(W(j)%e(edprev)),' and ',abs(W(j)%e(ed)) end if call matrixint(sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return ind = 2 * i - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return ind = 2 * abs(W(j)%e(ed)) - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),sites(1:2,abs(W(j)%e(edprev))),W(j)%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return ind = 2 * abs(W(j)%e(edprev)) - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) else if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed)) end if call line_exp2imp_2d(W(j)%v(1:2,ed),W(j)%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return ind = 2 * i - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return ind = 2 * abs(W(j)%e(ed)) - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) end if ! Ending point w_E if ( debug ) then write(*,*) write(*,*) 'Processing w_E' write(*,*) end if if ( W(j)%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W(j)%e(ed)),' and ',abs(W(j)%e(ednext)) end if call matrixint(sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return ind = 2 * i - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),sites(1:2,abs(W(j)%e(ednext))),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return ind = 2 * abs(W(j)%e(ednext)) - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixint(sites(1:2,abs(W(j)%e(ednext))),sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ed+1),M,flag) if ( flag .ne. 0 ) return ind = 2 * abs(W(j)%e(ed)) - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) else if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W(j)%e(ed)) end if call line_exp2imp_2d(W(j)%v(1:2,ednext),W(j)%v(1:2,ednext+1),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W(j)%e(ed))),sites(1:2,i),W(j)%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return ind = 2 * i - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W(j)%e(ed))),W(j)%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return ind = 2 * abs(W(j)%e(ed)) - 1 gval(ind:ind+1) = gval(ind:ind+1) + xi * mu * matmul( transpose( M(1:2,1:2) ), normal(1:2) ) end if else if ( debug ) then write(*,*) 'This is a fixed edge.' write(*,*) 'So, its vertices v_E and w_E are ignored.' end if end if end do else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j,' IS EMPTY.' write(*,*) end if end if end do end subroutine compJ2alternative2 ! ****************************************************************** ! ****************************************************************** subroutine fgcompJ2(x,fden,gden,compgrad,pdataptr) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad real(kind=8), intent(out) :: fden type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(2) real(kind=8), intent(out) :: gden(2) ! LOCAL SCALARS type(pdata_type), pointer :: pdata ! LOCAL ARRAY real(kind=8) :: z(2) call c_f_pointer(pdataptr,pdata) z(1:2) = (/ 2.0d0, 2.0d0 /) + 0.40d0 * x(1:2) if ( ( z(2) - ( 0.25d0 * z(1) ) ** 2 ) ** 2 + ( 0.25d0 * z(1) - 1.0d0 ) ** 2 - 1.0d0 .le. 0.0d0 ) then fden = 0.25d0 else fden = 1.075d0 end if if ( compgrad ) gden(1:2) = 0.0d0 !!$ fden = 2.5d0 - ( 2.5d0 - 0.5d0 ) * sum( ( x(1:2) - pdata%J2c(1:2) ) ** 2 ) / pdata%J2r ** 2 !!$ !!$ if ( compgrad ) gden(1:2) = - ( 2.5d0 - 0.5d0 ) * 2.0d0 * ( x(1:2) - pdata%J2c(1:2) ) / pdata%J2r ** 2 end subroutine fgcompJ2 ! ****************************************************************** ! ****************************************************************** subroutine compJ3(nsites,sites,i,W,ctol,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: ctol real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,edprev,ednext real(kind=8) :: glcnst,mu,perim,phia,phib,phic,sizeed,tmp,xi ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),tang(2) !write(*,*) 'Site i = ',i perim = 0.0d0 do ed = 1,W%n - 1 if ( ed .lt. W%n - 1 ) then ednext = ed + 1 else ednext = 1 end if sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) ) perim = perim + sizeed !write(*,*) 'Edge ',ed,' size = ',sizeed end do !write(*,*) 'Average size = ',perim / ( W%n - 1 ) if ( debug ) then write(*,*) write(*,*) 'Perimeter = ',perim write(*,*) end if xi = 2.0d0 / perim if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = 0.0d0 if ( compgrad ) then glcnst = 0.0d0 do ed = 1,W%n - 1 if ( ed .lt. W%n - 1 ) then ednext = ed + 1 else ednext = 1 end if sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) ) glcnst = glcnst + min( 0.0d0, ( W%n - 1 ) * sizeed / perim - ctol ) * sizeed / perim end do if ( debug ) then write(*,*) write(*,*) 'Global constant = ',glcnst write(*,*) end if gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /) gval(gnnz+1:gnnz+2) = 0.0d0 gnnz = gnnz + 2 end if do ed = 1,W%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed) write(*,*) '====================' write(*,*) end if ! Processing edge (v_E,w_E) if ( ed .eq. 1 ) then edprev = W%n - 1 ednext = ed + 1 else if ( ed .eq. W%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if if ( debug ) then write(*,*) 'Previous: ',edprev,' next: ',ednext write(*,*) 'They correspond to the points ',W%v(1:2,edprev),' and ',W%v(1:2,ednext),' respectively.' end if sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) ) if ( debug ) then write(*,*) 'Edge size = ',sizeed end if tmp = min( 0.0d0, ( W%n - 1 ) * sizeed / perim - ctol ) f = f + tmp ** 2 / ( W%n - 1 ) if ( .not. compgrad ) then cycle end if mu = glcnst - tmp if ( debug ) then write(*,*) 'mu associated wit the edge = ',mu end if if ( mu .eq. 0.0d0 ) then cycle end if tang(1:2) = ( W%v(1:2,ednext) - W%v(1:2,ed) ) / sizeed if ( debug ) then write(*,*) 'Tangent vector: ',tang(1:2) end if ! Starting point v_E if ( debug ) then write(*,*) write(*,*) 'Processing v_E' write(*,*) end if if ( W%e(ed) .lt. 0 .and. W%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(ed)) end if call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 else if ( W%e(ed) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed)) end if call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 else if ( W%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(edprev)) end if call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,ednext),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 end if ! Ending point w_E if ( debug ) then write(*,*) write(*,*) 'Processing w_E' write(*,*) end if if ( W%e(ed) .lt. 0 .and. W%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(ed)),' and ',abs(W%e(ednext)) end if call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(ednext))),sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 else if ( W%e(ed) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed)) end if call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ednext+1),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 else if ( W%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(ednext)) end if call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ed),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ3: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), tang(1:2) ) gnnz = gnnz + 2 end if end do end subroutine compJ3 ! ****************************************************************** ! ****************************************************************** subroutine compJ4(nsites,sites,i,W,ctol,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(in) :: ctol real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,edprev,edpp,ednext,nmov real(kind=8) :: ang,avgang,glcnst,mu,phia,phib,phic,sizeed,sizeedprev,tmp,xi ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),normaled(2),normaledprev(2),taued(2),tauedprev(2),zeta(2) !write(*,*) 'Site i = ',i nmov = 0 avgang = 0.0d0 do ed = 1,W%n - 1 if ( ed .eq. 1) then edprev = W%n - 1 ednext = ed + 1 else if ( ed .eq. W%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if ! If none is negative, the angle is fixed if ( W%e(ed) .lt. 0 .or. W%e(edprev) .lt. 0 ) then tauedprev(1:2) = W%v(1:2,ed) - W%v(1:2,edprev) taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed) sizeedprev = norm2( tauedprev(1:2) ) sizeed = norm2( taued(1:2) ) ang = acos( max( -1.0d0, min( dot_product( - tauedprev(1:2) / sizeedprev, taued(1:2) / sizeed ), 1.0d0 ) ) ) !write(*,*) 'vertex = ',ed,' angle = ',ang if ( debug ) then write(*,*) 'vertex = ',ed,' angle = ',ang end if nmov = nmov + 1 avgang = avgang + ang end if end do avgang = avgang / nmov !write(*,*) 'Average angle = ',avgang if ( debug ) then write(*,*) write(*,*) 'The polygon has nmov = ',nmov,' non-fixed angles.' write(*,*) 'Average angle = ',avgang write(*,*) end if xi = ( 2.0d0 / ( avgang * nmov ) ) if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = 0.0d0 if ( compgrad ) then glcnst = 0.0d0 do ed = 1,W%n - 1 if ( ed .eq. 1) then edprev = W%n - 1 ednext = ed + 1 else if ( ed .eq. W%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if ! If none is negative, the angle is fixed if ( W%e(ed) .lt. 0 .or. W%e(edprev) .lt. 0 ) then tauedprev(1:2) = W%v(1:2,ed) - W%v(1:2,edprev) taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed) sizeedprev = norm2( tauedprev(1:2) ) sizeed = norm2( taued(1:2) ) ang = acos( max( -1.0d0, min( dot_product( - tauedprev(1:2) / sizeedprev, taued(1:2) / sizeed ), 1.0d0 ) ) ) glcnst = glcnst + ( ang / avgang ) * min( 0.0d0, ang / avgang - ctol ) / nmov end if end do gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /) gval(gnnz+1:gnnz+2) = 0.0d0 gnnz = gnnz + 2 end if do ed = 1,W%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed) write(*,*) '====================' write(*,*) end if ! Processing edge (v_E,w_E) if ( ed .eq. 1) then edprev = W%n - 1 ednext = ed + 1 else if ( ed .eq. W%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if if ( debug ) then write(*,*) 'Previous: ',edprev,' next: ',ednext write(*,*) 'They correspond to the points ',W%v(1:2,edprev),' and ',W%v(1:2,ednext),' respectively.' end if ! If none is negative, the angle is fixed if ( W%e(ed) .lt. 0 .or. W%e(edprev) .lt. 0 ) then tauedprev(1:2) = W%v(1:2,ed) - W%v(1:2,edprev) taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed) sizeedprev = norm2( tauedprev(1:2) ) sizeed = norm2( taued(1:2) ) tauedprev(1:2) = tauedprev(1:2) / sizeedprev taued(1:2) = taued(1:2) / sizeed ang = acos( max( -1.0d0, min( dot_product( - tauedprev(1:2), taued(1:2) ), 1.0d0 ) ) ) tmp = min( 0.0d0, ang / avgang - ctol ) f = f + tmp ** 2 / nmov if ( .not. compgrad ) then cycle end if mu = glcnst - tmp if ( debug ) then write(*,*) 'mu = ',mu end if if ( mu .eq. 0.0d0 ) then cycle end if call line_exp_normal_2d(W%v(1:2,ednext),W%v(1:2,ed),normaled) call line_exp_normal_2d(W%v(1:2,ed),W%v(1:2,edprev),normaledprev) if ( debug ) then write(*,*) 'Normal to ed = ',normaled(1:2) write(*,*) 'Normal to edprev = ',normaledprev(1:2) end if ! Starting point v_E if ( debug ) then write(*,*) write(*,*) 'Processing v_E = w_{\hat E}' write(*,*) end if ! zeta(1:2) = ( dot_product( normaled(1:2), tauedprev(1:2) ) / sizeed ) * normaled(1:2) - & ! ( dot_product( normaledprev(1:2), taued(1:2) ) / sizeedprev ) * normaledprev(1:2) zeta(1:2) = normaled(1:2) / sizeed + normaledprev(1:2) / sizeedprev if ( debug ) then write(*,*) 'zeta = ',zeta(1:2) end if if ( W%e(ed) .lt. 0 .and. W%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(ed)) end if call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 else if ( W%e(ed) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed)) end if call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 else if ( W%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(edprev)) end if call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,ednext),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 end if ! Ending point w_E if ( debug ) then write(*,*) write(*,*) 'Processing w_E' write(*,*) end if ! zeta(1:2) = ( dot_product( normaled(1:2), tauedprev(1:2) ) / sizeed ) * normaled(1:2) zeta(1:2) = normaled(1:2) / sizeed if ( debug ) then write(*,*) 'zeta = ',zeta(1:2) end if if ( W%e(ed) .lt. 0 .and. W%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(ed)),' and ',abs(W%e(ednext)) end if call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(ednext))),sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 else if ( W%e(ed) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(ed)) end if call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ednext+1),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 else if ( W%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(ednext)) end if call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ed),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) - xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /) gval(gnnz+1:gnnz+2) = - xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 end if ! Starting point point v_{\hat E} if ( edprev .eq. 1 ) then edpp = W%n - 1 else edpp = edprev - 1 end if if ( debug ) then write(*,*) write(*,*) 'Processing v_{\hat E}' write(*,*) end if ! zeta(1:2) = ( dot_product( normaledprev(1:2), taued(1:2) ) / sizeedprev ) * normaledprev(1:2) zeta(1:2) = - normaledprev(1:2) / sizeedprev if ( debug ) then write(*,*) 'zeta = ',zeta(1:2) end if if ( W%e(edprev) .lt. 0 .and. W%e(edpp) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(edpp)) end if call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(edpp))),sites(1:2,i),W%v(1:2,edprev),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(edpp))),W%v(1:2,edprev),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edpp))-1, 2*abs(W%e(edpp)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(edpp))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,edprev),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 else if ( W%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2 (type 1). Neighbour cell is: ',abs(W%e(edprev)) end if call line_exp2imp_2d(W%v(1:2,edpp),W%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(edprev))),sites(1:2,i),W%v(1:2,edprev),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,edprev),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 else if ( W%e(edpp) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2 (type 2). Neighbour cell is: ',abs(W%e(edpp)) end if call line_exp2imp_2d(W%v(1:2,edprev),W%v(1:2,ed),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(edpp))),sites(1:2,i),W%v(1:2,edprev),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * mu * matmul( transpose( M(1:2,1:2) ),zeta(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(edpp))),W%v(1:2,edprev),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edpp))-1, 2*abs(W%e(edpp)) /) gval(gnnz+1:gnnz+2) = xi * mu * matmul( transpose( M(1:2,1:2) ), zeta(1:2) ) gnnz = gnnz + 2 end if else if ( debug ) then write(*,*) write(*,*) 'The inner vertex of this edge has a fixed angle; so it will be ignored.' write(*,*) end if end if end do end subroutine compJ4 ! ****************************************************************** ! ****************************************************************** subroutine compJ5(nsites,sites,i,W,f,gnnzmax,gnnz,gind,gval,flag,compgrad,debug) implicit none ! SCALAR ARGUMENTS logical, intent(in) :: compgrad,debug integer, intent(in) :: gnnzmax,i,nsites integer, intent(out) :: gnnz integer, intent(inout) :: flag real(kind=8), intent(out) :: f type(Polygon), intent(in) :: W ! ARRAY ARGUMENTS integer, intent(out) :: gind(gnnzmax) real(kind=8), intent(in) :: sites(2,nsites) real(kind=8), intent(out) :: gval(gnnzmax) ! LOCAL SCALARS integer :: ed,edprev,ednext,nmov real(kind=8) :: phia,phib,phic,sizeed,xi ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),mu(2),p(2),q(2),diff(2) nmov = count( W%e(1:W%n-1) .lt. 0 ) if ( debug ) then write(*,*) write(*,*) 'The polygon has nmov = ',nmov,' moving edges.' write(*,*) end if xi = 1.0d0 / nmov if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = 0.0d0 if ( compgrad ) then gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*i-1, 2*i /) gval(gnnz+1:gnnz+2) = 0.0d0 gnnz = gnnz + 2 end if do ed = 1,W%n - 1 if ( debug ) then write(*,*) write(*,*) '====================' write(*,*) 'Processing edge ed = ',ed write(*,*) 'Its inner vertex, named ',ed,' as well, corresponds to the point ',W%v(1:2,ed) write(*,*) '====================' write(*,*) end if ! Processing edge (v_E,w_E) if ( ed .eq. 1 ) then edprev = W%n - 1 ednext = ed + 1 else if ( ed .eq. W%n - 1 ) then edprev = ed - 1 ednext = 1 else edprev = ed - 1 ednext = ed + 1 end if if ( debug ) then write(*,*) 'Previous: ',edprev,' next: ',ednext write(*,*) 'They correspond to the points ',W%v(1:2,edprev),' and ',W%v(1:2,ednext),' respectively.' end if if ( W%e(ed) .lt. 0 ) then p(1:2) = 0.5d0 * ( W%v(1:2,ed) + W%v(1:2,ednext) ) q(1:2) = 0.5d0 * ( sites(1:2,i) + sites(1:2,abs(W%e(ed))) ) diff(1:2) = p(1:2) - q(1:2) sizeed = norm2( W%v(1:2,ednext) - W%v(1:2,ed) ) f = f + xi * sum( diff(1:2) ** 2 ) / sizeed ** 2 if ( .not. compgrad ) then cycle end if gval(1:2) = gval(1:2) - xi * diff(1:2) / sizeed ** 2 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * diff(1:2) / sizeed ** 2 gnnz = gnnz + 2 ! Starting point v_E if ( debug ) then write(*,*) write(*,*) 'Processing v_E' write(*,*) end if mu(1:2) = diff(1:2) / sizeed ** 2 + ( 2.0d0 * sum( diff(1:2) ** 2 ) / sizeed ** 4 ) * ( W%v(1:2,ednext) - W%v(1:2,ed) ) if ( all( mu(1:2) .eq. 0.0d0 ) ) then cycle end if if ( W%e(edprev) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(edprev)),' and ',abs(W%e(ed)) end if call matrixint(sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(edprev))),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,i),sites(1:2,abs(W%e(edprev))),W%v(1:2,ed),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(edprev))-1, 2*abs(W%e(edprev)) /) gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'We are in the case |P(v_E)|=2. Neighbour cell is: ',abs(W%e(ed)) end if call line_exp2imp_2d(W%v(1:2,ed),W%v(1:2,edprev),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ed),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) gnnz = gnnz + 2 end if ! Ending point w_E if ( debug ) then write(*,*) write(*,*) 'Processing w_E' write(*,*) end if mu(1:2) = diff(1:2) / sizeed ** 2 - ( 2.0d0 * sum( diff(1:2) ** 2 ) / sizeed ** 4 ) * ( W%v(1:2,ednext) - W%v(1:2,ed) ) if ( W%e(ednext) .lt. 0 ) then if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=3. ', & 'Neighbour cells are: ',abs(W%e(ed)),' and ',abs(W%e(ednext)) end if call matrixint(sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),sites(1:2,i),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) call matrixint(sites(1:2,i),sites(1:2,abs(W%e(ed))),sites(1:2,abs(W%e(ednext))),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ednext))-1, 2*abs(W%e(ednext)) /) gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) gnnz = gnnz + 2 call matrixint(sites(1:2,abs(W%e(ednext))),sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) gnnz = gnnz + 2 else if ( debug ) then write(*,*) 'We are in the case |P(w_E)|=2. Neighbour cell is: ',abs(W%e(ed)) end if call line_exp2imp_2d(W%v(1:2,ednext),W%v(1:2,ednext+1),phia,phib,phic) gradphi(1:2) = (/ phia, phib /) if ( gradphi(1) * sites(1,i) + gradphi(2) * sites(2,i) + phic .gt. 0.0d0 ) then gradphi(1:2) = - gradphi(1:2) end if call matrixbd(sites(1:2,abs(W%e(ed))),sites(1:2,i),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return gval(1:2) = gval(1:2) + xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) call matrixbd(sites(1:2,i),sites(1:2,abs(W%e(ed))),W%v(1:2,ednext),gradphi,M,flag) if ( flag .ne. 0 ) return if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compJ5: Increase gnnzmax to at least ',gnnz+2 flag = - 1 return end if gind(gnnz+1:gnnz+2) = (/ 2*abs(W%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = xi * matmul( transpose( M(1:2,1:2) ), mu(1:2) ) gnnz = gnnz + 2 end if end if end do end subroutine compJ5 ! ****************************************************************** ! ****************************************************************** subroutine matrixint(x,y,z,v,M,flag) ! It computes the matrix used in the calculation of the gradient ! when v belongs to T_int. implicit none ! SCALAR ARGUMENTS integer, intent(inout) :: flag ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(2),y(2),z(2),v(2) real(kind=8), intent(out) :: M(2,2) ! LOCAL SCALARS real(kind=8) :: det ! LOCAL ARRAYS real(kind=8) :: tmp(2) det = ( y(1) - x(1) ) * ( z(2) - x(2) ) - ( y(2) - x(2) ) * ( z(1) - x(1) ) if ( abs( det ) .lt. 1.0d+02 * epsilon( 1.0d0 ) ) then det = sign( 1.0d+02 * epsilon( 1.0d0 ), det ) end if ! This is to avoid a warning of the compiler that says flag is not ! being used. if ( .false. ) flag = - 1 !!$ if ( det .eq. 0.0d0 ) then !!$ write(*,*) 'In matrixint: Division by zero in the calculation of matrix M.' !!$ flag = - 1 !!$ return !!$ !!$ else if ( abs(det) .lt. 1.0d-14 ) then !!$ write(*,*) 'In matrixint: Warning: Determinant is too small.' !!$ end if tmp(1:2) = (/ y(2) - x(2), x(1) - y(1) /) M(1:2,1) = tmp(1:2) * ( v(1) - z(1) ) M(1:2,2) = tmp(1:2) * ( v(2) - z(2) ) M(1:2,1:2) = M(1:2,1:2) / det end subroutine matrixint ! ****************************************************************** ! ****************************************************************** subroutine matrixbd(y,x,v,gradphi,M,flag) ! It computes the matrix used in the calculation of the gradient ! when v belongs to T_bd. implicit none ! SCALAR ARGUMENTS integer, intent(inout) :: flag ! ARRAY ARGUMENTS real(kind=8), intent(in) :: gradphi(2),x(2),y(2),v(2) real(kind=8), intent(out) :: M(2,2) ! LOCAL SCALARS real(kind=8) :: det ! LOCAL ARRAYS real(kind=8) :: tmp(2) tmp(1:2) = (/ - gradphi(2), gradphi(1) /) det = - dot_product( tmp(1:2), y(1:2) - x(1:2) ) if ( abs( det ) .lt. 1.0d+02 * epsilon( 1.0d0 ) ) then det = sign( 1.0d+02 * epsilon( 1.0d0 ), det ) end if ! This is to avoid a warning of the compiler that says flag is not ! being used. if ( .false. ) flag = - 1 !!$ if ( det .eq. 0.0d0 ) then !!$ write(*,*) 'In matrixbd: Division by zero in the calculation of matrix M.' !!$ flag = - 1 !!$ return !!$ !!$ else if ( abs(det) .lt. 1.0d-14 ) then !!$ write(*,*) 'In matrixbd: Warning: Determinant is too small.' !!$ end if ! Outer product of (/ - gradphi(2), gradphi(1) /) with p(1:2) - x(1:2) M(1:2,1) = tmp(1:2) * ( v(1) - x(1) ) M(1:2,2) = tmp(1:2) * ( v(2) - x(2) ) M(1:2,1:2) = - M(1:2,1:2) / det end subroutine matrixbd ! ****************************************************************** ! ****************************************************************** subroutine drawsol1(n,x,filename,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) character(len=80) :: filename ! LOCAL SCALARS integer :: allocerr,i,ierror,j,k,maxnv,nsites,status,ntriag,v type(pdata_type), pointer :: pdata type(VoronoiCell) :: Vi type(Polygon) :: Wij,Pj,frame,Viframe ! LOCAL ARRAYS integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2)) call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .le. 1 ) then write(*,*) 'In drawsol: nsites must be at least 2!' stop end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror stop end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In drawsol: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' stop end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In drawsol: triangle_circumcenter_2d returned NaN.' write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k)) stop end if end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if open(unit=10,file=trim(filename)) write(10,10) 1,pdata%drawscale ! Drawing Delaunay triangulation do k = 1,ntriag do v = 1,3 write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen' end do end do ! Frame frame%n = 5 frame%deg = .false. allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /) frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /) frame%e(1:4) = (/ 0, 0, 0, 0 /) ! Drawing Voronoi diagram do i = 1,nsites call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,Vi,status) if ( status .ne. 0 ) then write(*,*) 'In drawsol: Error in voronoi_cell, status = ',status stop end if ! This is to draw the Voronoi cell restricted to the frame. if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe) end if if ( .not. Viframe%deg ) then do v = 1,Viframe%n - 1 write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil' end do end if call polygon_destroy(Viframe) ! This is to draw the Voronoi cell restricted to A. do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( Wij%e(v) .lt. 1 ) then write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire' end if end do end if call polygon_destroy(Wij) end do ! This is to draw the sites write(10,40) sites(1:2,i) call voronoi_cell_destroy(Vi) end do deallocate(frame%v,frame%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! Drawing A do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) do k = 1,Pj%n - 1 if ( Pj%e(k) .eq. 0 ) then write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black' end if end do end do ! Close figure write(10,50) close(10) deallocate(Pj%v,Pj%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! NON-EXECUTABLE STATEMENTS 10 format('prologues := 3;',/, & 'outputtemplate := "%j.mps";',/, & 'input mpcolornames;',/, & 'beginfig(',i4,');',/, & 'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', & 'orange, plum, bluea, blueb, bluec, blued, bluee, corn, deepsaffron, ferrarired;',/, & 'skyblue = (135/255,206/255,235/255);',/, & 'dollarbill = (133/255,187/255,101/255);',/, & 'teagreen = (208/255,240/255,192/255);',/, & 'desire = (234/255, 60/255, 83/255);',/, & 'jonquil = (244/255,202/255, 22/255);',/, & 'royalblue = ( 65/255,105/255,225/255);',/, & 'orange = (255/255,165/255, 0/255);',/, & 'plum = (221/255,160/255,221/255);',/, & 'bluea = ( 0/255,230/255,230/255);',/, & 'blueb = ( 0/255,191/255,230/255);',/, & 'bluec = ( 0/255,153/255,230/255);',/, & 'blued = ( 0/255,115/255,230/255);',/, & 'bluee = ( 0/255, 77/255,230/255);',/, & 'corn = (255/255,230/255,102/255);',/, & 'deepsaffron = (255/255,153/255, 51/255);',/, & 'ferrarired = (255/255, 42/255, 0/255);',/, & 'u = ',f20.10,'cm;') 20 format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';') 40 format('drawdot(',f20.10,'u,',f20.10,'u);') 50 format('endfig;',/,'end;') end subroutine drawsol1 ! ****************************************************************** ! ****************************************************************** subroutine drawsol2(n,x,filename,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) character(len=80) :: filename ! LOCAL SCALARS integer :: allocerr,i,ierror,j,k,maxnv,nsites,status,ntriag,v type(pdata_type), pointer :: pdata type(VoronoiCell) :: Vi type(Polygon) :: Wij,Pj,frame,Viframe ! LOCAL ARRAYS integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2)) integer, parameter :: nstats = 1000 integer :: flag,find logical :: debug real(kind=8) :: ctol,J3ctoltmp,J4ctoltmp integer :: cellcolorJ3(n/2) real(kind=8) :: f(0:5,n/2),stats(nstats,3) character(len=80) :: paintcolor,statsfn call c_f_pointer(pdataptr,pdata) J3ctoltmp = pdata%J3ctol J4ctoltmp = pdata%J4ctol do k = 1,nstats ctol = 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 ) pdata%J3ctol = ctol pdata%J4ctol = ctol debug = .false. call compallJ(n,x,f,debug,flag,pdataptr) if ( flag .ne. 0 ) then write(*,*) 'In drawsol, compallJ returned a non-null flag!' return end if do j = 1,3 if ( j .eq. 1 ) then find = 1 stats(k,j) = real( count( min( 0.0d0, f(find,1:n/2) + 1.0d0 - ctol ) ** 2 .le. 1.0d-04 ) ) / ( n / 2 ) else if ( j .eq. 2 ) then find = 3 stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 ) else ! if ( j .eq. 3 ) then find = 4 stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 ) end if end do end do write(statsfn,"(A6,F3.1,A1,F3.1,A4)") 'stats-',J3ctoltmp,'-',J4ctoltmp,'.dat' open(unit=10,file=trim(statsfn)) do k = 1,nstats write(10,1000) 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 ),stats(k,1:3) end do close(10) ! The smaller the ctol for which the cell has positive f, the darker the color. cellcolorJ3(1:n/2) = - 1 ! Means no color do k = 5,1,-1 ctol = 0.1d0 * k pdata%J3ctol = ctol debug = .false. call compallJ(n,x,f,debug,flag,pdataptr) if ( flag .ne. 0 ) then write(*,*) 'In drawsol, compallJ returned a non-null flag!' return end if where ( f(3,1:n/2) .gt. 1.0d-04 ) cellcolorJ3(1:n/2) = k end where end do pdata%J3ctol = J3ctoltmp pdata%J4ctol = J4ctoltmp call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .le. 1 ) then write(*,*) 'In drawsol: nsites must be at least 2!' stop end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror stop end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In drawsol: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' stop end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In drawsol: triangle_circumcenter_2d returned NaN.' write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k)) stop end if end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if open(unit=10,file=trim(filename)) write(10,10) 1,pdata%drawscale ! Drawing Delaunay triangulation do k = 1,ntriag do v = 1,3 write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen' end do end do ! Frame frame%n = 5 frame%deg = .false. allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /) frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /) frame%e(1:4) = (/ 0, 0, 0, 0 /) ! Drawing Voronoi diagram do i = 1,nsites call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,Vi,status) if ( status .ne. 0 ) then write(*,*) 'In drawsol: Error in voronoi_cell, status = ',status stop end if ! This is to draw the Voronoi cell restricted to the frame. if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe) end if if ( .not. Viframe%deg ) then do v = 1,Viframe%n - 1 write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil' end do end if call polygon_destroy(Viframe) ! This is to paint the Voronoi cell, restricted to A, if ! desired conditions does not hold. if ( cellcolorJ3(i) .ne. - 1 ) then if ( cellcolorJ3(i) .eq. 5 ) then paintcolor = 'bluea' else if ( cellcolorJ3(i) .eq. 4 ) then paintcolor = 'blueb' else if ( cellcolorJ3(i) .eq. 3 ) then paintcolor = 'bluec' else if ( cellcolorJ3(i) .eq. 2 ) then paintcolor = 'blued' else if ( cellcolorJ3(i) .eq. 1 ) then paintcolor = 'bluee' end if do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( v .eq. 1 ) then write(10,30) Wij%v(1:2,v),Wij%v(1:2,v+1) else if ( v .gt. 1 .and. v .lt. Wij%n - 1 ) then write(10,31) Wij%v(1:2,v),Wij%v(1:2,v+1) else write(10,32) Wij%v(1:2,v),trim(paintcolor) end if end do end if call polygon_destroy(Wij) end do end if ! This is to draw the Voronoi cell restricted to A. do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( Wij%e(v) .lt. 1 ) then write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire' end if end do end if call polygon_destroy(Wij) end do ! This is to draw the sites write(10,40) sites(1:2,i) call voronoi_cell_destroy(Vi) end do deallocate(frame%v,frame%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! Drawing A do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) do k = 1,Pj%n - 1 if ( Pj%e(k) .eq. 0 ) then write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black' end if end do end do ! Close figure adding a few statistics write(10,*) '% Area of A = ',pdata%volA write(10,*) '% Number of polygons that define A = ',pdata%npols write(10,*) '% Number of sites = ',nsites write(10,50) close(10) deallocate(Pj%v,Pj%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! NON-EXECUTABLE STATEMENTS 10 format('prologues := 3;',/, & 'outputtemplate := "%j.mps";',/, & 'input mpcolornames;',/, & 'beginfig(',i4,');',/, & 'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', & 'orange, plum, bluea, blueb, bluec, blued, bluee, corn, deepsaffron, ferrarired;',/, & 'skyblue = (135/255,206/255,235/255);',/, & 'dollarbill = (133/255,187/255,101/255);',/, & 'teagreen = (208/255,240/255,192/255);',/, & 'desire = (234/255, 60/255, 83/255);',/, & 'jonquil = (244/255,202/255, 22/255);',/, & 'royalblue = ( 65/255,105/255,225/255);',/, & 'orange = (255/255,165/255, 0/255);',/, & 'plum = (221/255,160/255,221/255);',/, & 'bluea = ( 0/255,230/255,230/255);',/, & 'blueb = ( 0/255,191/255,230/255);',/, & 'bluec = ( 0/255,153/255,230/255);',/, & 'blued = ( 0/255,115/255,230/255);',/, & 'bluee = ( 0/255, 77/255,230/255);',/, & 'corn = (255/255,230/255,102/255);',/, & 'deepsaffron = (255/255,153/255, 51/255);',/, & 'ferrarired = (255/255, 42/255, 0/255);',/, & 'u = ',f20.10,'cm;') 20 format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';') 30 format('fill ((',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--') 31 format(' (',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--') 32 format(' (',f20.10,'u,',f20.10,'u)--cycle) withcolor ',a30,';') 40 format('drawdot(',f20.10,'u,',f20.10,'u);') 50 format('endfig;',/,'end;') 1000 format(4(1X,F10.4)) end subroutine drawsol2 ! ****************************************************************** ! ****************************************************************** subroutine drawsol3(n,x,filename,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) character(len=80) :: filename ! LOCAL SCALARS integer :: allocerr,i,ierror,j,k,maxnv,nsites,status,ntriag,v type(pdata_type), pointer :: pdata type(VoronoiCell) :: Vi type(Polygon) :: Wij,Pj,frame,Viframe ! LOCAL ARRAYS integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: sites(2,n/2),vorxy(2,2*(n/2)) integer, parameter :: nstats = 1000 integer :: flag,find logical :: debug real(kind=8) :: ctol,J3ctoltmp,J4ctoltmp integer :: cellcolorJ4(n/2) real(kind=8) :: f(0:5,n/2),stats(nstats,3) character(len=80) :: paintcolor,statsfn call c_f_pointer(pdataptr,pdata) J3ctoltmp = pdata%J3ctol J4ctoltmp = pdata%J4ctol do k = 1,nstats ctol = 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 ) pdata%J3ctol = ctol pdata%J4ctol = ctol debug = .false. call compallJ(n,x,f,debug,flag,pdataptr) if ( flag .ne. 0 ) then write(*,*) 'In drawsol, compallJ returned a non-null flag!' return end if do j = 1,3 if ( j .eq. 1 ) then find = 1 stats(k,j) = real( count( min( 0.0d0, f(find,1:n/2) + 1.0d0 - ctol ) ** 2 .le. 1.0d-04 ) ) / ( n / 2 ) else if ( j .eq. 2 ) then find = 3 stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 ) else ! if ( j .eq. 3 ) then find = 4 stats(k,j) = real( count( f(find,1:n/2) .le. 1.0d-04 ) ) / ( n / 2 ) end if end do end do write(statsfn,"(A6,F3.1,A1,F3.1,A4)") 'stats-',J3ctoltmp,'-',J4ctoltmp,'.dat' open(unit=10,file=trim(statsfn)) do k = 1,nstats write(10,1000) 1.0d0 - ( nstats - real( k ) ) / ( nstats - 1 ),stats(k,1:3) end do close(10) ! The smaller the ctol for which the cell has positive f, the darker the color. cellcolorJ4(1:n/2) = - 1 ! Means no color do k = 7,5,-1 ctol = 0.1d0 * k pdata%J4ctol = ctol debug = .false. call compallJ(n,x,f,debug,flag,pdataptr) if ( flag .ne. 0 ) then write(*,*) 'In drawsol, compallJ returned a non-null flag!' return end if where ( f(4,1:n/2) .gt. 1.0d-04 ) cellcolorJ4(1:n/2) = k end where end do pdata%J3ctol = J3ctoltmp pdata%J4ctol = J4ctoltmp call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .le. 1 ) then write(*,*) 'In drawsol: nsites must be at least 2!' stop end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror stop end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In drawsol: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' stop end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In drawsol: triangle_circumcenter_2d returned NaN.' write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k)) stop end if end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if open(unit=10,file=trim(filename)) write(10,10) 1,pdata%drawscale ! Drawing Delaunay triangulation do k = 1,ntriag do v = 1,3 write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen' end do end do ! Frame frame%n = 5 frame%deg = .false. allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /) frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /) frame%e(1:4) = (/ 0, 0, 0, 0 /) ! Drawing Voronoi diagram do i = 1,nsites call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,Vi,status) ! call voronoi_cell_old(i,nsites,sites,ntriag,til,tnbr,vorxy,Vi,status) if ( status .ne. 0 ) then write(*,*) 'In drawsol: Error in voronoi_cell, status = ',status stop end if ! This is to draw the Voronoi cell restricted to the frame. if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe) end if if ( .not. Viframe%deg ) then do v = 1,Viframe%n - 1 write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil' end do end if call polygon_destroy(Viframe) ! This is to paint the Voronoi cell, restricted to A, if ! desired conditions does not hold. if ( cellcolorJ4(i) .ne. - 1 ) then if ( cellcolorJ4(i) .eq. 7 ) then paintcolor = 'corn' else if ( cellcolorJ4(i) .eq. 6 ) then paintcolor = 'deepsaffron' else if ( cellcolorJ4(i) .eq. 5 ) then paintcolor = 'ferrarired' end if do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( v .eq. 1 ) then write(10,30) Wij%v(1:2,v),Wij%v(1:2,v+1) else if ( v .gt. 1 .and. v .lt. Wij%n - 1 ) then write(10,31) Wij%v(1:2,v),Wij%v(1:2,v+1) else write(10,32) Wij%v(1:2,v),trim(paintcolor) end if end do end if call polygon_destroy(Wij) end do end if ! This is to draw the Voronoi cell restricted to A. do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( Wij%e(v) .lt. 1 ) then write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire' end if end do end if call polygon_destroy(Wij) end do ! This is to draw the sites write(10,40) sites(1:2,i) call voronoi_cell_destroy(Vi) end do deallocate(frame%v,frame%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! Drawing A do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) do k = 1,Pj%n - 1 if ( Pj%e(k) .eq. 0 ) then write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black' end if end do end do ! Close figure adding a few statistics write(10,*) '% Area of A = ',pdata%volA write(10,*) '% Number of polygons that define A = ',pdata%npols write(10,*) '% Number of sites = ',nsites write(10,50) close(10) deallocate(Pj%v,Pj%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! NON-EXECUTABLE STATEMENTS 10 format('prologues := 3;',/, & 'outputtemplate := "%j.mps";',/, & 'input mpcolornames;',/, & 'beginfig(',i4,');',/, & 'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', & 'orange, plum, bluea, blueb, bluec, blued, bluee, corn, deepsaffron, ferrarired;',/, & 'skyblue = (135/255,206/255,235/255);',/, & 'dollarbill = (133/255,187/255,101/255);',/, & 'teagreen = (208/255,240/255,192/255);',/, & 'desire = (234/255, 60/255, 83/255);',/, & 'jonquil = (244/255,202/255, 22/255);',/, & 'royalblue = ( 65/255,105/255,225/255);',/, & 'orange = (255/255,165/255, 0/255);',/, & 'plum = (221/255,160/255,221/255);',/, & 'bluea = ( 0/255,230/255,230/255);',/, & 'blueb = ( 0/255,191/255,230/255);',/, & 'bluec = ( 0/255,153/255,230/255);',/, & 'blued = ( 0/255,115/255,230/255);',/, & 'bluee = ( 0/255, 77/255,230/255);',/, & 'corn = (255/255,230/255,102/255);',/, & 'deepsaffron = (255/255,153/255, 51/255);',/, & 'ferrarired = (255/255, 42/255, 0/255);',/, & 'u = ',f20.10,'cm;') 20 format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';') 30 format('fill ((',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--') 31 format(' (',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u)--') 32 format(' (',f20.10,'u,',f20.10,'u)--cycle) withcolor ',a30,';') 40 format('drawdot(',f20.10,'u,',f20.10,'u);') 50 format('endfig;',/,'end;') 1000 format(4(1X,F10.4)) end subroutine drawsol3 ! ****************************************************************** ! ****************************************************************** subroutine drawsol4(n,x,filename,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) character(len=80) :: filename ! LOCAL SCALARS logical :: debug integer :: allocerr,flag,i,ierror,j,k,maxnv,nsites,status,ntriag,v type(pdata_type), pointer :: pdata type(VoronoiCell) :: Vi type(Polygon) :: Wij,Pj,frame,Viframe ! LOCAL ARRAYS integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: f(0:5,n/2),sites(2,n/2),vorxy(2,2*(n/2)) debug = .false. call compallJ(n,x,f,debug,flag,pdataptr) if ( flag .ne. 0 ) then write(*,*) 'In drawsol, compallJ returned a non-null flag!' return end if call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .le. 1 ) then write(*,*) 'In drawsol: nsites must be at least 2!' stop end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror stop end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In drawsol: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' stop end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In drawsol: triangle_circumcenter_2d returned NaN.' write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k)) stop end if end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if open(unit=10,file=trim(filename)) write(10,10) 1,pdata%drawscale ! Drawing Delaunay triangulation do k = 1,ntriag do v = 1,3 write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen' end do end do ! Frame frame%n = 5 frame%deg = .false. allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /) frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /) frame%e(1:4) = (/ 0, 0, 0, 0 /) ! Drawing Voronoi diagram do i = 1,nsites call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,Vi,status) if ( status .ne. 0 ) then write(*,*) 'In drawsol: Error in voronoi_cell, status = ',status stop end if ! This is to draw the Voronoi cell restricted to the frame. if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe) end if if ( .not. Viframe%deg ) then do v = 1,Viframe%n - 1 write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil' end do end if call polygon_destroy(Viframe) ! This is to draw the Voronoi cell restricted to A. do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( Wij%e(v) .lt. 1 ) then write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire' end if ! This is to draw the distance between the midpoint of ! Delaunay and Voronoy corresponding edges if ( Wij%e(v) .lt. 0 ) then write(10,20) 0.5d0 * ( sites(1:2,i) + sites(1:2,abs(Wij%e(v))) ), & 0.5d0 * ( Wij%v(1:2,v) + Wij%v(1:2,v+1) ),'royalblue' end if end do end if call polygon_destroy(Wij) end do ! This is to draw the sites write(10,40) sites(1:2,i) call voronoi_cell_destroy(Vi) end do deallocate(frame%v,frame%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! Drawing A do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) do k = 1,Pj%n - 1 if ( Pj%e(k) .eq. 0 ) then write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black' end if end do end do ! Close figure write(10,50) write(10,*) 'Value of J0i (sum/ave/min/max) = ',sum(abs(f(0,1:n/2))),sum(abs(f(0,1:n/2)))/(n/2), & minval(abs(f(0,1:n/2))),maxval(abs(f(0,1:n/2))) write(10,*) 'Value of J1i (sum/ave/min/max) = ',sum(abs(f(1,1:n/2))),sum(abs(f(1,1:n/2)))/(n/2), & minval(abs(f(1,1:n/2))),maxval(abs(f(1,1:n/2))),' J1 itself = ',sum(f(1,1:n/2)**2)/(n/2) write(10,*) 'Value of J5i (sum/ave/min/max) = ',sum(abs(f(5,1:n/2))),sum(abs(f(5,1:n/2)))/(n/2), & minval(abs(f(5,1:n/2))),maxval(abs(f(5,1:n/2))) close(10) deallocate(Pj%v,Pj%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! NON-EXECUTABLE STATEMENTS 10 format('prologues := 3;',/, & 'outputtemplate := "%j.mps";',/, & 'input mpcolornames;',/, & 'beginfig(',i4,');',/, & 'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', & 'orange, plum, bluea, blueb, bluec, blued, bluee, corn, deepsaffron, ferrarired;',/, & 'skyblue = (135/255,206/255,235/255);',/, & 'dollarbill = (133/255,187/255,101/255);',/, & 'teagreen = (208/255,240/255,192/255);',/, & 'desire = (234/255, 60/255, 83/255);',/, & 'jonquil = (244/255,202/255, 22/255);',/, & 'royalblue = ( 65/255,105/255,225/255);',/, & 'orange = (255/255,165/255, 0/255);',/, & 'plum = (221/255,160/255,221/255);',/, & 'bluea = ( 0/255,230/255,230/255);',/, & 'blueb = ( 0/255,191/255,230/255);',/, & 'bluec = ( 0/255,153/255,230/255);',/, & 'blued = ( 0/255,115/255,230/255);',/, & 'bluee = ( 0/255, 77/255,230/255);',/, & 'corn = (255/255,230/255,102/255);',/, & 'deepsaffron = (255/255,153/255, 51/255);',/, & 'ferrarired = (255/255, 42/255, 0/255);',/, & 'u = ',f20.10,'cm;') 20 format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';') 40 format('drawdot(',f20.10,'u,',f20.10,'u);') 50 format('endfig;',/,'end;') end subroutine drawsol4 ! ****************************************************************** ! ****************************************************************** subroutine drawsol(n,x,filename,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) character(len=80) :: filename ! LOCAL SCALARS logical :: debug integer :: allocerr,flag,i,ierror,j,k,maxnv,nsites,status,ntriag,v type(pdata_type), pointer :: pdata type(VoronoiCell) :: Vi type(Polygon) :: Wij,Pj,frame,Viframe ! LOCAL ARRAYS integer :: indices(n/2),sitetv(3*2*(n/2)),start((n/2)+1),til(3,2*(n/2)),tnbr(3,2*(n/2)) real(kind=8) :: f(0:5,n/2),sites(2,n/2),vorxy(2,2*(n/2)) debug = .false. call compallJ(n,x,f,debug,flag,pdataptr) if ( flag .ne. 0 ) then write(*,*) 'In drawsol, compallJ returned a non-null flag!' return end if call c_f_pointer(pdataptr,pdata) nsites = n / 2 indices(1:nsites) = (/ (i,i=1,nsites) /) sites(1:2,1:nsites) = reshape( x(1:n), (/ 2, nsites /) ) if ( nsites .le. 1 ) then write(*,*) 'In drawsol: nsites must be at least 2!' stop end if ! Compute Delaunay triangulation call dtris2(nsites,sites,indices,ntriag,til,tnbr,ierror) if ( ierror .ne. 0 .and. ierror .ne. 225 ) then write(*,*) 'In drawsol: dtris2 returned an error different from 0 and 225: ',ierror stop end if call voronoi_cell_patch(nsites,ntriag,til,start,sitetv) ! Compute vertices of the Voronoi diagram if ( ntriag .gt. 2 * nsites ) then write(*,*) 'In drawsol: There is something wrong. According to the dtris2, ' write(*,*) 'documentation the number of triangles in the Delaunay trangulation ' write(*,*) 'cannot be larger than twice the number of sites!' stop end if do k = 1,ntriag call triangle_circumcenter_2d(sites(1:2,til(1:3,k)),vorxy(1:2,k)) if ( isnan(vorxy(1,k)) .or. isnan(vorxy(1,k)) ) then write(*,*) 'In drawsol: triangle_circumcenter_2d returned NaN.' write(*,*) 'Triangle vertices: ',sites(1:2,til(1:3,k)) stop end if end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(Pj%v(2,maxnv),Pj%e(maxnv-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if open(unit=10,file=trim(filename)) write(10,10) 1,pdata%drawscale ! Drawing Delaunay triangulation do k = 1,ntriag do v = 1,3 write(10,20) sites(1:2,til(v,k)),sites(1:2,til(mod(v,3)+1,k)),'teagreen' end do end do ! Frame frame%n = 5 frame%deg = .false. allocate(frame%v(2,frame%n),frame%e(frame%n-1),stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Allocation error.' stop end if frame%v(1,1:5) = (/ pdata%xmin, pdata%xmax, pdata%xmax, pdata%xmin, pdata%xmin /) frame%v(2,1:5) = (/ pdata%ymin, pdata%ymin, pdata%ymax, pdata%ymax, pdata%ymin /) frame%e(1:4) = (/ 0, 0, 0, 0 /) ! Drawing Voronoi diagram do i = 1,nsites call voronoi_cell(i,nsites,sites,start,sitetv,ntriag,til,tnbr,vorxy,Vi,status) if ( status .ne. 0 ) then write(*,*) 'In drawsol: Error in voronoi_cell, status = ',status stop end if ! This is to draw the Voronoi cell restricted to the frame. if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),frame,Vi,Viframe) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,frame,Viframe) end if if ( .not. Viframe%deg ) then do v = 1,Viframe%n - 1 write(10,20) Viframe%v(1:2,v),Viframe%v(1:2,v+1),'jonquil' end do end if call polygon_destroy(Viframe) ! This is to draw the Voronoi cell restricted to A. do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),Pj,Vi,Wij) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,Pj,Wij) end if if ( .not. Wij%deg ) then do v = 1,Wij%n - 1 if ( Wij%e(v) .lt. 1 ) then write(10,20) Wij%v(1:2,v),Wij%v(1:2,v+1),'desire' end if end do end if call polygon_destroy(Wij) end do ! This is to draw the sites write(10,40) sites(1:2,i) call voronoi_cell_destroy(Vi) end do deallocate(frame%v,frame%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! Drawing A do j = 1,pdata%npols Pj%n = pdata%nvpols(j) Pj%deg = .false. Pj%v(1,1:Pj%n) = pdata%xpol(j,1:Pj%n) Pj%v(2,1:Pj%n) = pdata%ypol(j,1:Pj%n) Pj%e(1:Pj%n-1) = pdata%poledges(j,1:Pj%n-1) do k = 1,Pj%n - 1 if ( Pj%e(k) .eq. 0 ) then write(10,20) Pj%v(1:2,k),Pj%v(1:2,k+1),'black' end if end do end do ! Close figure write(10,50) write(10,*) 'Value of J2i (sum/ave/min/max) = ',sum(abs(f(2,1:n/2))),sum(abs(f(2,1:n/2)))/(n/2), & minval(abs(f(2,1:n/2))),maxval(abs(f(2,1:n/2))),' J2 itself = ',sum(f(2,1:n/2)**2)/(n/2) close(10) deallocate(Pj%v,Pj%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In drawsol: Deallocation error.' stop end if ! NON-EXECUTABLE STATEMENTS 10 format('prologues := 3;',/, & 'outputtemplate := "%j.mps";',/, & 'input mpcolornames;',/, & 'beginfig(',i4,');',/, & 'color skyblue, dollarbill, desire, jonquil, royalblue, teagreen, ', & 'orange, plum, bluea, blueb, bluec, blued, bluee, corn, deepsaffron, ferrarired;',/, & 'skyblue = (135/255,206/255,235/255);',/, & 'dollarbill = (133/255,187/255,101/255);',/, & 'teagreen = (208/255,240/255,192/255);',/, & 'desire = (234/255, 60/255, 83/255);',/, & 'jonquil = (244/255,202/255, 22/255);',/, & 'royalblue = ( 65/255,105/255,225/255);',/, & 'orange = (255/255,165/255, 0/255);',/, & 'plum = (221/255,160/255,221/255);',/, & 'bluea = ( 0/255,230/255,230/255);',/, & 'blueb = ( 0/255,191/255,230/255);',/, & 'bluec = ( 0/255,153/255,230/255);',/, & 'blued = ( 0/255,115/255,230/255);',/, & 'bluee = ( 0/255, 77/255,230/255);',/, & 'corn = (255/255,230/255,102/255);',/, & 'deepsaffron = (255/255,153/255, 51/255);',/, & 'ferrarired = (255/255, 42/255, 0/255);',/, & 'u = ',f20.10,'cm;') 20 format('draw ( ',f20.10,'u,',f20.10,'u)--(',f20.10,'u,',f20.10,'u) withcolor',a12,';') 40 format('drawdot(',f20.10,'u,',f20.10,'u);') 50 format('endfig;',/,'end;') end subroutine drawsol end module optvor