module optvormod2 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,J2ctol,J3ctol,J4c(2),J4r integer, allocatable :: nvpols(:),poledges(:,:) real(kind=8), allocatable :: xpol(:,:),ypol(:,:) real(kind=8) :: Jweight(0:4) ! 0: CVT ! 1: cells with identical volume given by vol(A)/m ! 2: each cell with equally sized edges ! 3: each cell with equally sized angles ! 4: cells with volume proportional to density function f evaluated at their sites end type pdata_type public :: evalfg,approxg,drawsol,drawsolJ1,drawsolJ2 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 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 integer :: allocerr,gtmpnnz,i,ierror,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P,W ! 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)),valf(n/2) 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(4) .gt. 0.0d0 ) then do i = 1,nsites call fgcompJ4(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(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 P%n = pdata%nvpols(1) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(1,1:P%n) P%v(2,1:P%n) = pdata%ypol(1,1:P%n) P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1) 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 if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W) end if if ( pdata%Jweight(4) .gt. 0.0d0) then call compJ4(nsites,sites,i,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0) return f = f + pdata%Jweight(4) * ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(4) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz) end if if ( W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.' end if f = f + 1.0d0 else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) end if call compG(nsites,sites,i,W,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) ! call compGalternative(nsites,sites,i,W,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) 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) if ( pdata%Jweight(1) .gt. 0.0d0) then call compJ1(nsites,sites,i,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,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f = f + pdata%Jweight(2) * ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(2) * gtmpval(1:gtmpnnz) end if if ( pdata%Jweight(3) .gt. 0.0d0 ) then call compJ3(nsites,sites,i,W,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 end if call voronoi_cell_destroy(V) call polygon_destroy(W) end do f = f / nsites if ( compgrad ) g(1:n) = g(1:n) / nsites deallocate(P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine evalfg ! ****************************************************************** ! ****************************************************************** !! This function works to calculate the values of G, J1, J2, J3 subroutine evalJg(n,x,G1,j2,j3,j1,j4,ji,g,flag,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n integer, intent(out) :: flag real(kind=8), intent(out) :: j2,G1,j3,j1,j4 type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: g(n),ji(n/2) ! PARAMETERS logical, parameter :: compgrad = .false., debug = .false. integer, parameter :: gtmpnnzmax = 100000 ! LOCAL SCALARS integer :: allocerr,gtmpnnz,i,ierror,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P,W ! 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)),valf(n/2) 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 fgcompJ4(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr) end do G1 = 0.0d0 j1 = 0.0d0 j2 = 0.0d0 j3 = 0.0d0 j4 = 0.0d0 if ( compgrad ) g(1:n) = 0.0d0 maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(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 P%n = pdata%nvpols(1) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(1,1:P%n) P%v(2,1:P%n) = pdata%ypol(1,1:P%n) P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1) 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 if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W) end if call compJ4(nsites,sites,i,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return ! write(*,*) 'ftmpevalJg = ', ftmp j4 = j4 + ftmp ** 2 if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(4) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz) if ( W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.' end if ! write(*,*) 'Degenerate polygon, do not anything.' G1 = G1 + 1.0d0 j1 = j1 + 1.0d0 j2 = j2 + 1.0d0 j3 = j3 + 1.0d0 else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) end if call compG(nsites,sites,i,W,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return G1 = G1 + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + gtmpval(1:gtmpnnz) call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return j1 = j1 + ftmp ** 2 ji(i) = ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(1) * 2.0d0 * ftmp * gtmpval(1:gtmpnnz) call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return j2 = j2 + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(1) * gtmpval(1:gtmpnnz) call compJ3(nsites,sites,i,W,pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return j3 = j3 + ftmp if ( compgrad ) g(gtmpind(1:gtmpnnz)) = g(gtmpind(1:gtmpnnz)) + pdata%Jweight(2) * gtmpval(1:gtmpnnz) end if call voronoi_cell_destroy(V) call polygon_destroy(W) end do G1 = G1 / nsites j2 = j2 / nsites j3 = j3 / nsites j1 = j1 / nsites j4 = j4 / nsites if ( compgrad ) g(1:n) = g(1:n) / nsites deallocate(P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine evalJg !! This function works to calculate the volume of each cell, and the J^1_i(a) recalling that !! the function J^1(a)=1/nsites sum^{nsites}_{i=1}[J^1_i(a)]^2. subroutine evalvolJi1(n,x,vol,ji1,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(in) :: x(n) real(kind=8), intent(out) :: ji1(n/2), vol(n/2) ! PARAMETERS logical, parameter :: compgrad = .false., debug = .false. integer, parameter :: gtmpnnzmax = 100000 ! LOCAL SCALARS integer :: allocerr,i,ierror,k,maxnv,nsites,ntriag,status,ed,ednext,gtmpnnz real(kind=8) :: perim,ftmp,volW type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P,W ! 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) :: sites(2,n/2),vorxy(2,2*(n/2)),vsizeed(10),taued(2),gtmpval(gtmpnnzmax) 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 maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(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 P%n = pdata%nvpols(1) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(1,1:P%n) P%v(2,1:P%n) = pdata%ypol(1,1:P%n) P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1) 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 if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W) end if if ( W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.' end if write(*,*) 'Degenerate polygon, do not anything.' else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) end if call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),volW) vol(i) = volW call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return ji1(i) = ftmp end if call voronoi_cell_destroy(V) call polygon_destroy(W) end do deallocate(P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine evalvolJi1 ! ******************************************************************* ! ******************************************************************* !! This function works to calculate the average of each cell, the size of the edges of each cell and the J^2_i(a) recalling that !! the function J^2(a)=1/nsites sum^{nsites}_{i=1}J^2_i(a). subroutine evalsizeedJi2(n,x,avg,sizes,ji2,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(in) :: x(n) real(kind=8), intent(out) :: ji2(n/2), avg(n/2), sizes(n/2,10) ! PARAMETERS logical, parameter :: compgrad = .false., debug = .false. integer, parameter :: gtmpnnzmax = 100000 ! LOCAL SCALARS integer :: allocerr,i,ierror,k,maxnv,nsites,ntriag,status,ed,ednext,gtmpnnz real(kind=8) :: perim,ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P,W ! 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) :: sites(2,n/2),vorxy(2,2*(n/2)),vsizeed(10),taued(2),gtmpval(gtmpnnzmax) 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 maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(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 P%n = pdata%nvpols(1) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(1,1:P%n) P%v(2,1:P%n) = pdata%ypol(1,1:P%n) P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1) 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 if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W) end if if ( W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.' end if write(*,*) 'Degenerate polygon, do not anything.' else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) end if vsizeed(1:10) = 0.0d0 perim = 0.0d0 do ed = 1,W%n - 1 if ( ed .eq. W%n-1) then ednext = 1 else ednext = ed + 1 end if taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed) vsizeed(ed) = norm2( taued(1:2) ) perim = perim + vsizeed(ed) end do sizes(i,1:10)=vsizeed(1:10) avg(i) = perim / (W%n-1) call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return ji2(i) = ftmp end if call voronoi_cell_destroy(V) call polygon_destroy(W) end do deallocate(P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine evalsizeedJi2 ! ******************************************************************* ! ******************************************************************* !! This subroutine give us the minimun value of |E| and minimum value of \theta (angle between two edges) in the whole diagram !! Also gives a vector with the average size of the cell's of the edges of each Voronoi cell. subroutine eval(n,x,e,a,avgsizeed,flag,pdataptr) implicit none ! SCALAR ARGUMENTS integer, intent(in) :: n integer, intent(out) :: flag real(kind=8), intent(out) :: e,a type(c_ptr), optional, intent(in) :: pdataptr ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(n) real(kind=8), intent(out) :: avgsizeed(n/2) ! PARAMETERS logical, parameter :: debug = .false. ! LOCAL SCALARS integer :: allocerr,i,ierror,k,maxnv,nsites,ntriag,status,ed,ednext,edprev real(kind=8) :: perim type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P,W ! 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) :: fden(n/2),gden(2,n/2),sites(2,n/2),vorxy(2,2*(n/2)),valf(n/2),taued(2),tauedprev(2), & vminang(n/2),vminsizeed(n/2),vsizeed(10),vang(10) 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 maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(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 P%n = pdata%nvpols(1) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(1,1:P%n) P%v(2,1:P%n) = pdata%ypol(1,1:P%n) P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1) 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 if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W) end if if ( W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.' end if write(*,*) 'Degenerate polygon, do not anything.' else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) end if perim = 0.0d0 do ed=1,W%n-1 if ( ed .eq. W%n-1 ) then edprev = ed - 1 ednext = 1 else if ( ed .eq. 1) then edprev = W%n-1 ednext = ed + 1 else edprev = ed - 1 ednext = ed + 1 end if taued(1:2) = W%v(1:2,ednext) - W%v(1:2,ed) tauedprev(1:2)= W%v(1:2,ed) - W%v(1:2,edprev) vsizeed(ed) = norm2( taued(1:2) ) vang(ed) = acos( max( -1.0d0, min( dot_product( -tauedprev(1:2)/norm2( tauedprev(1:2) ), taued(1:2)/vsizeed(ed) ), 1.0d0) ) ) perim = perim + vsizeed(ed) end do vminang(i) = minval( vang(1:W%n-1) ) vminsizeed(i) = minval( vsizeed(1:W%n-1) ) avgsizeed(i) = perim/(W%n-1) if ( flag .ne. 0 ) return end if call voronoi_cell_destroy(V) call polygon_destroy(W) end do e = minval(vminsizeed(1:nsites)) a = minval(vminang(1:nsites)) deallocate(P%v,P%e,stat=allocerr) if ( allocerr .ne. 0 ) then write(*,*) 'In evalfg: Deallocation error.' flag = - 1 return end if end subroutine eval ! ****************************************************************** ! ****************************************************************** subroutine compG(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,ednext,nx,ny,j,k,r real(kind=8) :: vdet,volW,tmp,tinteg,tintegex ! LOCAL ARRAYS real(kind=8) :: diff1(2),diff2(2),cent(2),u(2),diff(nsites),v1(2),v2(2),v3(2),v4(3) integer :: minindex(1) ! LOCAL PARAMETERS real(kind=8), parameter :: h = 1.0d-4 if ( compgrad ) then call centroid(W,cent,volW) if ( debug ) then write(*,*) write(*,*) 'volWi= ',volW write(*,*) 'centi= ',cent(1:2) end if gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compG: 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.0 gnnz = gnnz + 2 gval(1:2) = gval(1:2) + 2.0d0 * volW * ( sites(1:2,i) - cent(1:2) ) end if f = 0.0d0 do ed = 1,W%n-1 if ( ed .eq. W%n-1 ) then ednext = 1 else ednext = ed + 1 end if 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(*,*) 'Its outer vertex, named ',ednext,'as well, corresponds to the point ',W%v(1:2,ednext) write(*,*) '====================' write(*,*) end if if ( debug ) then write(*,*) 'Next: ',ednext write(*,*) 'It corresponds to the point',W%v(1:2,ednext) end if diff1(1:2) = W%v(1:2,ed) - sites(1:2,i) diff2(1:2) = W%v(1:2,ednext) - sites(1:2,i) vdet = ( diff1(1) * diff2(2) ) - ( diff1(2) * diff2(1) ) ! vdet = ( W%v(1,ed) - sites(1,i) ) * ( W%v(2,ednext) - sites(2,i) ) - & ! ( W%v(2,ed) - sites(2,i) ) * ( W%v(1,ednext) - sites(1,i) ) !call det(W%v(1:2,ed), sites(1:2,i), W%v(1:2,ednext),vdet) if ( debug ) then write(*,*) write(*,*)'v_E - a_i = ',diff1(1:2) write(*,*)'w_E - a_i = ',diff2(1:2) write(*,*)'det = ',vdet write(*,*)'||v_E-a_i||^2= ', sum( diff1(1:2)**2 ) write(*,*)'(v_E-a_i)^T(w_E-a_i)', dot_product( diff1(1:2),diff2(1:2) ) write(*,*)'||w_E-a_i||^2= ',sum( diff2(1:2)**2 ) write(*,*) end if tmp = ( abs( vdet ) ) * ( sum( diff1(1:2)**2 ) + & dot_product( diff1(1:2),diff2(1:2) ) + sum( diff2(1:2)**2 ) ) f = f + tmp if ( debug ) then write(*,*) write(*,*)'tmp = ', tmp write(*,*)'f = ', f write(*,*) end if end do f = f /12.0d0 if ( debug ) then write(*,*)'f= ', f end if end subroutine compG ! ****************************************************************** ! ****************************************************************** subroutine normalvec(v,normv) implicit none ! ARRAY ARGUMENTS real(kind=8), intent(in) :: v(2) real(kind=8), intent(out) :: normv(2) normv(1:2) = (/ -v(2), v(1) /) end subroutine normalvec ! ****************************************************************** ! ****************************************************************** subroutine compGalternative(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 real(kind=8) :: phia,phib,phic,xi,vdet,tmp,kap,absvdet,kappe,sig ! LOCAL ARRAYS real(kind=8) :: gradphi(2),M(2,2),wrot(2),vrot(2),p(2),diff1(2),diff2(2),ai(2),q(2) xi = 1.0d0/12.0d0 f = 0.0d0 if ( debug ) then write(*,*) write(*,*) 'xi = ', xi write(*,*) end if if ( compgrad ) then gnnz = 0 if ( gnnz + 2 .gt. gnnzmax ) then write(*,*) 'In compGalternative: 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 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 ! Calculating the objective function diff1(1:2) = W%v(1:2,ed) - sites(1:2,i) diff2(1:2) = W%v(1:2,ednext) - sites(1:2,i) vdet = ( diff1(1) * diff2(2) ) - ( diff1(2) * diff2(1) ) if ( debug ) then write(*,*) write(*,*)'v_E - a_i = ',diff1(1:2) write(*,*)'w_E - a_i = ',diff2(1:2) write(*,*)'det = ',vdet write(*,*)'||v_E-a_i||^2= ', sum( diff1(1:2)**2 ) write(*,*)'(v_E-a_i)^T(w_E-a_i)', dot_product( diff1(1:2),diff2(1:2) ) write(*,*)'||w_E-a_i||^2= ',sum( diff2(1:2)**2 ) write(*,*) end if kap = sum( diff1(1:2)**2 ) + dot_product( diff1(1:2),diff2(1:2) ) + sum( diff2(1:2)**2 ) absvdet = abs(vdet) tmp = absvdet * kap f = f + tmp ! f = f + absvdet if ( debug ) then write(*,*) write(*,*)'tmp = ', tmp write(*,*)'f = ', f write(*,*) end if if ( compgrad ) then ! Processing edge (v_E,w_E) if ( debug ) then write(*,*) 'Previous: ', edprev, 'next: ',ednext write(*,*) 'They corresponds to the points ',W%v(1:2,edprev),' and ', W%v(1:2,ednext),' respectively.' end if vrot(1:2) = (/ -W%v(2,ed), W%v(1,ed) /) wrot(1:2) = (/ W%v(2,ednext), -W%v(1,ednext) /) ai(1:2) = (/ -sites(2,i), sites(1,i)/) kappe = sign(1.0d0,vdet) * kap p(1:2) = kappe*(wrot(1:2)+ai(1:2))+absvdet*(2.0d0*W%v(1:2,ed) + W%v(1:2,ednext)-3.0d0*sites(1:2,i)) q(1:2) = kappe*(vrot(1:2)-ai(1:2))+absvdet*(W%v(1:2,ed)+2.0d0*W%v(1:2,ednext)-3.0d0*sites(1:2,i)) sig = sign(1.0d0,vdet) if ( debug ) then write(*,*) write(*,*) 'vrot= ',vrot write(*,*) 'wrot= ',wrot write(*,*) 'ai*= ', ai write(*,*) 'k(E)= ',kappe write(*,*) end if gval(1:2) = gval(1:2) - xi*(3.0d0*absvdet * (W%v(1:2,ed) + W%v(1:2,ednext)-2.0d0*sites(1:2,i)) + & kappe*(wrot(1:2)+vrot(1:2))) ! if ( i .eq. 1 ) then ! write(*,*)'first term just depends of \delta ai.' ! write(*,*)'summed(1:2)= ', - xi*(absvdet * (9.0d0*W%v(1:2,ed) + 5.0d0*W%v(1:2,ednext)-14.0d0*sites(1:2,i)) + & ! kappe*(wrot(1:2)+vrot(1:2))) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if ! gval(1:2) = gval(1:2) - sig * xi * ( wrot(1:2) + vrot(1:2)) ! 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 * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) ) ! if ( i .eq. 1 ) then ! write(*,*)'Case |P(v_E)|=3.' ! write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if 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 compGalternative: 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) ), p(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(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 compGalternative: 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) ), p(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(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 * matmul( transpose( M(1:2,1:2) ), p(1:2) ) !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) ) ! if ( i .eq. 1 ) then ! write(*,*)'Case |P(v_E)|=2 (type 1).' ! write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if 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 compGalternative: 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) ), p(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(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 * matmul( transpose( M(1:2,1:2) ), p(1:2) ) !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(1:2) ) ! if ( i .eq. 1 ) then ! write(*,*)'Case |P(v_E)|=2 (type 2).' ! write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if 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 compGalternative: 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) ), p(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), wrot(1:2) + ai(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 * matmul( transpose( M(1:2,1:2) ), q(1:2) ) !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) ) ! if ( i .eq. 1 ) then ! write(*,*)'Case |P(w_E)|=3.' ! write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if 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 compGalternative: 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) ), q(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(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 compGalternative: 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) ), q(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(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 * matmul( transpose( M(1:2,1:2) ), q(1:2) ) !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) ) ! if ( i .eq. 1 ) then ! write(*,*)'Case |P(w_E)|=2 (type 1).' ! write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if 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 compGalternative: 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) ), q(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(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 * matmul( transpose( M(1:2,1:2) ), q(1:2) ) !gval(1:2) = gval(1:2) + sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) ) ! if ( i .eq. 1 ) then ! write(*,*)'Case |P(w_E)|=2 (type 2).' ! write(*,*)'summed(1:2)= ', xi * matmul( transpose( M(1:2,1:2) ), p(1:2) ) ! write(*,*)'gi(1:2)= ',gval(1:2) ! write(*,*) ! end if 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 compGalternative: 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) ), q(1:2) ) !gval(gnnz+1:gnnz+2) = sig * xi * matmul( transpose( M(1:2,1:2) ), vrot(1:2) - ai(1:2) ) gnnz = gnnz + 2 end if end if end do f = xi * f if ( debug ) then write(*,*) 'f = ',f end if end subroutine compGalternative ! ****************************************************************** ! ****************************************************************** 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:3,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 integer :: allocerr,gtmpnnz,i,ierror,k,maxnv,nsites,ntriag,status real(kind=8) :: ftmp type(pdata_type), pointer :: pdata type(VoronoiCell) :: V type(Polygon) :: P,W ! 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)) 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 fgcompJ4(sites(1:2,i),fden(i),gden(1:2,i),compgrad,pdataptr) end do maxnv = maxval( pdata%nvpols(1:pdata%npols) ) allocate(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 P%n = pdata%nvpols(1) P%deg = .false. P%v(1,1:P%n) = pdata%xpol(1,1:P%n) P%v(2,1:P%n) = pdata%ypol(1,1:P%n) P%e(1:P%n-1) = pdata%poledges(1,1:P%n-1) 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 if ( nsites .ge. 3 .and. ierror .eq. 0 ) then call VorCellInterConvPol(sites(1:2,i),P,V,W) else ! if ( ( nsites .ge. 3 .and. ierror .eq. 225 ) .or. nsites .eq. 2 ) then call VorCellInterConvPol_Collinear(i,nsites,sites,P,W) end if call compJ4(nsites,sites,i,W,pdata%volA,fden,gden,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(3,i) = ftmp if ( W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i, ' IS EMPTY.' end if f(0,i) = 1.0d0 f(1,i) = 1.0d0 f(2,i) = 1.0d0 ! f(3,i) = 1.0d0 else if ( debug ) then write(*,*) write(*,*) 'POLYGON Wi with i = ',i write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) end if call compJ1(nsites,sites,i,W,pdata%volA,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(0,i) = ftmp call compJ2(nsites,sites,i,W,pdata%J2ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(1,i) = ftmp call compJ3(nsites,sites,i,W,pdata%J3ctol,ftmp,gtmpnnzmax,gtmpnnz,gtmpind,gtmpval,flag,compgrad,debug) if ( flag .ne. 0 ) return f(2,i) = ftmp end if call voronoi_cell_destroy(V) call polygon_destroy(W) end do deallocate(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 compJ1(nsites,sites,i,W,volA,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) :: volA 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,ednext,j real(kind=8) :: mu,vol,xi ! LOCAL ARRAYS real(kind=8) :: p(2) call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),vol) 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 if ( .not. W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) 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(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W%v(1:2,ednext) write(*,*) '====================' write(*,*) end if if ( W%e(ed) .lt. 0 ) then if ( ed .eq. W%n - 1 ) then ednext = 1 else ednext = ed + 1 end if if ( debug ) then write(*,*) 'Next: ',ednext write(*,*) 'It corresponds to the point ',W%v(1:2,ednext) end if mu = norm2( W%v(1:2,ed) - W%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W%e(ed))) ) if ( debug ) then write(*,*) 'mu = ',mu end if p(1:2) = 0.5d0 * ( W%v(1:2,ed) + W%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%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W%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 Wi with i = ',i,' IS EMPTY.' write(*,*) end if end if end subroutine compJ1 ! ****************************************************************** ! ****************************************************************** subroutine compJ2(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 compJ2 ! ****************************************************************** ! ****************************************************************** subroutine compJ4(nsites,sites,i,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,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 ! 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,xi ! LOCAL ARRAYS real(kind=8) :: p(2) call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),vol) xi = nsites / volA if ( debug ) then write(*,*) write(*,*) 'xi = ',xi write(*,*) end if f = xi * vol - fden(i) if ( .not. compgrad ) then return end if if ( 2 .gt. gnnzmax ) then write(*,*) 'In compJ4: Increase gnnzmax to at least ',2 flag = - 1 return end if gind(1:2) = (/ 2*i-1, 2*i /) gval(1:2) = -gden(1:2,i) gnnz = 2 if ( .not. W%deg ) then if ( debug ) then write(*,*) write(*,*) 'POLYGON Wij with i = ',i,' and j = ',j write(*,*) 'Number of vertices: ',W%n - 1 write(*,*) 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(*,*) 'Its outer vertex, named ',ednext,' as well, corresponds to the point ',W%v(1:2,ednext) write(*,*) '====================' write(*,*) end if if ( W%e(ed) .lt. 0 ) then if ( ed .eq. W%n - 1 ) then ednext = 1 else ednext = ed + 1 end if if ( debug ) then write(*,*) 'Next: ',ednext write(*,*) 'It corresponds to the point ',W%v(1:2,ednext) end if mu = norm2( W%v(1:2,ed) - W%v(1:2,ednext) ) / norm2( sites(1:2,i) - sites(1:2,abs(W%e(ed))) ) if ( debug ) then write(*,*) 'mu = ',mu end if p(1:2) = 0.5d0 * ( W%v(1:2,ed) + W%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%e(ed))-1, 2*abs(W%e(ed)) /) gval(gnnz+1:gnnz+2) = - xi * mu * ( p(1:2) - sites(1:2,abs(W%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 Wi with i = ',i,' IS EMPTY.' write(*,*) end if end if end subroutine compJ4 ! ****************************************************************** ! ****************************************************************** 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,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 compJ3 ! ****************************************************************** ! ****************************************************************** subroutine fgcompJ4(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 real(kind=8) :: delta,param,a,b,scale,z1,z2,faux real(kind=8), parameter :: pi = 4.0d0 * atan( 1.0d0 ) ! LOCAL ARRAY real(kind=8) :: z(2) call c_f_pointer(pdataptr,pdata) !! Function psi_1 ! z(1:2) = -pdata%c(1:2)/5.0d0 + 0.40d0 * x(1:2) ! b = 0.25d0 ! a = 0.0625d0 * (2.0d0 - b - 0.5625d0) ! fden = a * ( (z(2) - ( 0.25d0 * z(1) ) ** 2 ) ** 2 + ( 0.25d0 * z(1) - 1.0d0 ) ** 2 ) + b ! if ( compgrad ) then ! gden(1) = a * ( -0.1d0 * z(1) * ( z(2) - ( 0.25d0 * z(1) ) ** 2 ) + 0.2d0 * (0.25d0 * z(1) - 1.0d0) ) ! gden(2) = a * ( 0.8d0 * ( z(2) - ( 0.25d0 * z(1) ) ** 2 ) ) ! end if !! Function psi_2 delta = 0.5d0*sqrt( 1000.d0 ) param = 2.9d0/delta**2 fden = 0.1d0 + param*( x(2) - 0.6d0*delta*sin( 2.0d0*pi*x(1)/(2.0d0*delta) ) - delta ) ** 2 if ( compgrad ) then gden(1) = 2.0d0 * param * ( x(2) - 0.6d0*delta*sin( 2.0d0*pi*x(1)/(2.0d0*delta) ) ) & * ( -0.6*2.0d0*pi*cos( 2.0d0*pi*x(1)/(2.0d0*delta) )/(2.0d0*delta) ) gden(2) = 2.0d0 * param * ( x(2) - 0.6d0*delta*sin( 2.0d0*pi*x(1)/(2.0d0*delta) ) ) end if !! Function 3 psi_3 ! fden = 0.01d0 + 20.0d0 * sum( ( x(1:2) - pdata%J4c(1:2) ) ** 2 ) / pdata%J4r ** 2 ! if ( compgrad ) gden(1:2) = 40.0d0 * ( x(1:2) - pdata%J4c(1:2) ) / pdata%J4r ** 2 end subroutine fgcompJ4 ! ****************************************************************** ! ****************************************************************** subroutine centroid(W,cent,volW) implicit none ! SCALAR ARGUMENTS type(polygon), intent(in) :: W real(kind=8), intent(out) :: volW ! ARRAY ARGUMENTS real(kind=8), intent(out) :: cent(2) ! LOCAL SCALARS integer :: j real(kind=8) :: re call polygon_area_2d(W%n-1,W%v(1:2,1:W%n-1),volW) cent = 0.0d0 do j=1,W%n-1 if ( j .eq. W%n-1 ) then re = W%v(1,j)*W%v(2,1) - W%v(1,1)*W%v(2,j) cent = cent + re*(W%v(1:2,j) + W%v(1:2,1)) else re = W%v(1,j)*W%v(2,j+1) - W%v(1,j+1)*W%v(2,j) cent = cent + re*(W%v(1:2,j) + W%v(1:2,j+1)) end if end do cent = cent/(6.0d0*volW) end subroutine centroid ! ****************************************************************** ! ****************************************************************** subroutine det(x,y,z,vdet) implicit none ! SACALAR ARGUMENTS real(kind=8), intent(out) :: vdet ! ARRAY ARGUMENTS real(kind=8), intent(in) :: x(2),y(2),z(2) vdet = ( x(1) - y(1) ) * ( z(2) - y(2) ) - ( x(2) - y(2) ) * ( z(1) - y(1) ) end subroutine det ! ****************************************************************** ! ****************************************************************** 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 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:3,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 ! write(*,*) 'ftmpscompallJ= ', f(3,1:n/2) ! write(*,*) 'J4compall = ', sum( f(3,1:n/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 !! For the paper it is not necessary to draw the 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 J4i (sum/ave/min/max) = ',sum(abs(f(3,1:n/2))),sum(abs(f(3,1:n/2)))/(n/2), & minval(abs(f(3,1:n/2))),maxval(abs(f(3,1:n/2))),' J4 itself = ',sum(f(3,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 ! ****************************************************************** ! ****************************************************************** subroutine drawsolJ1(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 :: flag,find logical :: debug real(kind=8) :: ctol integer :: cellcolorJ1(n/2) real(kind=8) :: f(0:3,n/2) character(len=80) :: paintcolor call c_f_pointer(pdataptr,pdata) ! The smaller the ctol for which the cell has positive f, the darker the color. cellcolorJ1(1:n/2) = -1 ! Means no color 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 ( abs( f(0,1:n/2) ) .gt. 1.0d-03 ) cellcolorJ1(1:n/2) = 1 end where 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 ( cellcolorJ1(i) .ne. - 1 ) then if ( cellcolorJ1(i) .eq. 1 ) then paintcolor = 'dollarbill' 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 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;') end subroutine drawsolJ1 ! ****************************************************************** ! ****************************************************************** subroutine drawsolJ2(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 :: flag,find logical :: debug real(kind=8) :: ctol, J2ctoltmp integer :: cellcolorJ2(n/2) real(kind=8) :: f(0:3,n/2) character(len=80) :: paintcolor call c_f_pointer(pdataptr,pdata) J2ctoltmp = pdata%J2ctol ! The smaller the ctol for which the cell has positive f, the darker the color. cellcolorJ2(1:n/2) = -1 ! Means no color do k = 8,1,-1 ctol = 0.1d0 * k pdata%J2ctol = 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(1,1:n/2) .gt. 1.0d-04 ) cellcolorJ2(1:n/2) = k end where end do pdata%J2ctol = J2ctoltmp 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 ( cellcolorJ2(i) .ne. - 1 ) then if ( cellcolorJ2(i) .eq. 8 ) then paintcolor = 'blueaa' else if ( cellcolorJ2(i) .eq. 7 ) then paintcolor = 'bluebb' else if ( cellcolorJ2(i) .eq. 6 ) then paintcolor = 'bluecc' else if ( cellcolorJ2(i) .eq. 5 ) then paintcolor = 'bluedd' else if ( cellcolorJ2(i) .eq. 4 ) then paintcolor = 'blueee' else if ( cellcolorJ2(i) .eq. 3 ) then paintcolor = 'bluef' else if ( cellcolorJ2(i) .eq. 2 ) then paintcolor = 'blueg' else if ( cellcolorJ2(i) .eq. 1 ) then paintcolor = 'blueh' 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 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, blueaa, bluebb, bluecc, ', & 'bluedd, blueee, bluef, blueg, blueh, bluei, 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);',/, & 'blueaa = ( 0/255,255/255,230/255);',/, & 'bluebb = ( 0/255,235/255,230/255);',/, & 'bluecc = ( 0/255,205/255,230/255);',/, & 'bluedd = ( 0/255,175/255,230/255);',/, & 'blueee = ( 0/255,145/255,230/255);',/, & 'bluef = ( 0/255,115/255,230/255);',/, & 'blueg = ( 0/255, 85/255,230/255);',/, & 'blueh = ( 0/255, 55/255,230/255);',/, & 'bluei = ( 0/255, 25/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;') end subroutine drawsolJ2 end module optvormod2