! D06DBF Example Program Text ! Mark 23 Release. NAG Copyright 2011. MODULE d06dbfe_mod ! D06DBF Example Program Module: ! Parameters and User-defined Routines ! .. Use Statements .. USE nag_library, ONLY : nag_wp ! .. Implicit None Statement .. IMPLICIT NONE ! .. Parameters .. INTEGER, PARAMETER :: nin = 5, nout = 6 CONTAINS FUNCTION fbnd(i,x,y,ruser,iuser) ! .. Implicit None Statement .. IMPLICIT NONE ! .. Function Return Value .. REAL (KIND=nag_wp) :: fbnd ! .. Scalar Arguments .. REAL (KIND=nag_wp), INTENT (IN) :: x, y INTEGER, INTENT (IN) :: i ! .. Array Arguments .. REAL (KIND=nag_wp), INTENT (INOUT) :: ruser(*) INTEGER, INTENT (INOUT) :: iuser(*) ! .. Local Scalars .. REAL (KIND=nag_wp) :: radius2, x0, y0 ! .. Executable Statements .. fbnd = 0.0_nag_wp SELECT CASE (i) CASE (1) ! inner circle x0 = 0.0_nag_wp y0 = 0.0_nag_wp radius2 = 1.0_nag_wp fbnd = (x-x0)**2 + (y-y0)**2 - radius2 CASE (2) ! outer circle x0 = 0.0_nag_wp y0 = 0.0_nag_wp radius2 = 5.0_nag_wp fbnd = (x-x0)**2 + (y-y0)**2 - radius2 END SELECT RETURN END FUNCTION fbnd END MODULE d06dbfe_mod PROGRAM d06dbfe ! D06DBF Example Main Program ! .. Use Statements .. USE d06dbfe_mod, ONLY : nout ! .. Implicit None Statement .. IMPLICIT NONE ! .. Executable Statements .. WRITE (nout,*) 'D06DBF Example Program Results' CALL ex1 CALL ex2 CONTAINS SUBROUTINE ex1 ! .. Use Statements .. USE nag_library, ONLY : d06daf, d06dbf, nag_wp USE d06dbfe_mod, ONLY : nin ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: eps INTEGER :: i, ifail, imax, itrace, & itrans, jmax, jtrans, k, & ktrans, liwork, lrwork, & nedge1, nedge2, nedge3, & nelt1, nelt2, nelt3, ntrans, & nv1, nv2, nv3 CHARACTER (1) :: pmesh ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: coor1(:,:), coor2(:,:), & coor3(:,:), rwork(:), & trans(:,:) INTEGER, ALLOCATABLE :: conn1(:,:), conn2(:,:), & conn3(:,:), edge1(:,:), & edge2(:,:), edge3(:,:), & itype(:), iwork(:), reft1(:), & reft2(:), reft3(:) ! .. Intrinsic Functions .. INTRINSIC real ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) 'Example 1' WRITE (nout,*) ! Skip heading in data file READ (nin,*) READ (nin,*) imax = 20 jmax = imax nedge1 = 2*(imax-1) + 2*(jmax-1) nedge2 = nedge1 nedge3 = nedge1 + nedge2 ntrans = 1 lrwork = 12*ntrans ! Read the mesh : coordinates and connectivity of the 1st domain READ (nin,*) nv1, nelt1 nv2 = nv1 nv3 = nv1 + nv2 nelt2 = nelt1 nelt3 = nelt1 + nelt2 liwork = 2*nv1 + 3*nv2 + nelt1 + nelt2 + nedge1 + nedge2 + 1024 ALLOCATE (coor1(2,nv1),coor2(2,nv2),coor3(2,nv3),conn1(3,nelt1), & conn2(3,nelt2),conn3(3,nelt3),reft1(nelt1),reft2(nelt2), & reft3(nelt3),edge1(3,nedge1),edge2(3,nedge2),edge3(3,nedge3), & itype(ntrans),trans(6,ntrans),rwork(lrwork),iwork(liwork)) DO i = 1, nv1 READ (nin,*) coor1(1:2,i) END DO DO k = 1, nelt1 READ (nin,*) conn1(1:3,k) END DO reft1(1:nelt1) = 1 reft2(1:nelt2) = 2 READ (nin,*) pmesh ! the Edges of the boundary DO i = 1, imax - 1 edge1(1,i) = i edge1(2,i) = i + 1 END DO DO i = 1, jmax - 1 edge1(1,imax-1+i) = i*imax edge1(2,imax-1+i) = (i+1)*imax END DO DO i = 1, imax - 1 edge1(1,imax-1+jmax-1+i) = imax*jmax - i + 1 edge1(2,imax-1+jmax-1+i) = imax*jmax - i END DO DO i = 1, jmax - 1 edge1(1,2*(imax-1)+jmax-1+i) = (jmax-i)*imax + 1 edge1(2,2*(imax-1)+jmax-1+i) = (jmax-i-1)*imax + 1 END DO edge1(3,1:nedge1) = 1 DO ktrans = 1, 2 ! Translation of the 1st domain to obtain the 2nd domain ! KTRANS = 1 leading to a domains overlapping ! KTRANS = 2 leading to a domains partition IF (ktrans==1) THEN itrans = imax - 5 jtrans = jmax - 3 ELSE itrans = imax - 1 jtrans = 0 END IF itype(1:ntrans) = (/ 1/) trans(1,1:ntrans) = (/ real(itrans,kind=nag_wp)/ & real(imax-1,kind=nag_wp) /) trans(2,1:ntrans) = (/ real(jtrans,kind=nag_wp)/ & real(jmax-1,kind=nag_wp) /) itrace = 0 ifail = 0 CALL d06daf(nv2,nedge2,nelt2,ntrans,itype,trans,coor1,edge1, & conn1,coor2,edge2,conn2,itrace,rwork,lrwork,ifail) edge2(3,1:nedge2) = 2 ! Call to the restitching driver itrace = 0 eps = 1.E-2_nag_wp ifail = 0 CALL d06dbf(eps,nv1,nelt1,nedge1,coor1,edge1,conn1,reft1,nv2, & nelt2,nedge2,coor2,edge2,conn2,reft2,nv3,nelt3,nedge3,coor3, & edge3,conn3,reft3,itrace,iwork,liwork,ifail) SELECT CASE (pmesh) CASE ('N') IF (ktrans==1) THEN WRITE (nout,*) 'The restitched mesh characteristics' WRITE (nout,*) 'in the overlapping case' ELSE WRITE (nout,*) 'in the partition case' END IF WRITE (nout,99999) 'NV =', nv3 WRITE (nout,99999) 'NELT =', nelt3 WRITE (nout,99999) 'NEDGE =', nedge3 CASE ('Y') ! Output the mesh WRITE (nout,99998) nv3, nelt3, nedge3 DO i = 1, nv3 WRITE (nout,99997) coor3(1:2,i) END DO DO k = 1, nelt3 WRITE (nout,99996) conn3(1:3,k), reft3(k) END DO DO k = 1, nedge3 WRITE (nout,99998) edge3(1:3,k) END DO CASE DEFAULT WRITE (nout,*) 'Problem with the printing option Y or N' END SELECT END DO 99999 FORMAT (1X,A,I6) 99998 FORMAT (1X,3I10) 99997 FORMAT (2(2X,E13.6)) 99996 FORMAT (1X,4I10) END SUBROUTINE ex1 SUBROUTINE ex2 ! .. Use Statements .. USE nag_library, ONLY : d06abf, d06acf, d06baf, d06caf, d06daf, & d06dbf, f16dnf, nag_wp USE d06dbfe_mod, ONLY : fbnd, nin ! .. Implicit None Statement .. IMPLICIT NONE ! .. Local Scalars .. REAL (KIND=nag_wp) :: eps INTEGER :: i, ifail, itrace, j, k, & liwork, lrwork, maxind, & maxval, ncomp, nedge1, & nedge2, nedge3, nedge4, & nedge5, nedmx, nelt1, nelt2, & nelt3, nelt4, nelt5, nlines, & npropa, nqint, ntrans, nv1, & nv2, nv3, nv4, nv5, nvb1, & nvb2, nvfix, nvint, nvmax, & sdcrus CHARACTER (1) :: pmesh ! .. Local Arrays .. REAL (KIND=nag_wp), ALLOCATABLE :: coor1(:,:), coor2(:,:), & coor3(:,:), coor4(:,:), & coor5(:,:), coorch(:,:), & crus(:,:), rate(:), rwork(:), & trans(:,:), weight(:) REAL (KIND=nag_wp) :: ruser(1) INTEGER, ALLOCATABLE :: conn1(:,:), conn2(:,:), & conn3(:,:), conn4(:,:), & conn5(:,:), edge1(:,:), & edge2(:,:), edge3(:,:), & edge4(:,:), edge5(:,:), & itype(:), iwork(:), lcomp(:), & lined(:,:), nlcomp(:), & numfix(:), reft1(:), & reft2(:), reft3(:), reft4(:), & reft5(:) INTEGER :: iuser(1) ! .. Intrinsic Functions .. INTRINSIC abs ! .. Executable Statements .. WRITE (nout,*) WRITE (nout,*) 'Example 2' WRITE (nout,*) ! Skip heading in data file READ (nin,*) ! Build the mesh of the 1st domain ! Initialise boundary mesh inputs: ! the number of line and of the characteristic points of ! the boundary mesh READ (nin,*) nlines, nvmax, nedmx ALLOCATE (coor1(2,nvmax),edge1(3,nedmx),lined(4,nlines), & lcomp(nlines),coorch(2,nlines),rate(nlines)) ! Characteristic points of the boundary geometry READ (nin,*) coorch(1,1:nlines) READ (nin,*) coorch(2,1:nlines) ! The Lines of the boundary mesh READ (nin,*) (lined(1:4,j),rate(j),j=1,nlines) sdcrus = 0 DO i = 1, nlines IF (lined(4,i)<0) THEN sdcrus = sdcrus + lined(1,i) - 2 END IF END DO liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus ! Get max(LINED(1,:)) for computing LRWORK CALL f16dnf(nlines,lined,4,maxind,maxval) lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines ! The number of connected components to the boundary ! and their informations READ (nin,*) ncomp ALLOCATE (nlcomp(ncomp),crus(2,sdcrus),rwork(lrwork),iwork(liwork)) j = 1 DO i = 1, ncomp READ (nin,*) nlcomp(i) k = j + abs(nlcomp(i)) - 1 READ (nin,*) lcomp(j:k) j = k + 1 END DO itrace = 0 ! Call to the 2D boundary mesh generator ifail = 0 CALL d06baf(nlines,coorch,lined,fbnd,crus,sdcrus,rate,ncomp,nlcomp, & lcomp,nvmax,nedmx,nvb1,coor1,nedge1,edge1,itrace,ruser,iuser, & rwork,lrwork,iwork,liwork,ifail) DEALLOCATE (rwork,iwork) ! mesh it using Delaunay-Voronoi method ! Initialise mesh control parameters itrace = 0 npropa = 1 nvint = 0 lrwork = 12*nvmax + 15 liwork = 6*nedge1 + 32*nvmax + 2*nvb1 + 78 ALLOCATE (weight(nvint),rwork(lrwork),iwork(liwork), & conn1(3,2*nvmax+5)) ! Call to the 2D Delaunay-Voronoi mesh generator ifail = 0 CALL d06abf(nvb1,nvint,nvmax,nedge1,edge1,nv1,nelt1,coor1,conn1, & weight,npropa,itrace,rwork,lrwork,iwork,liwork,ifail) DEALLOCATE (rwork,iwork) ! Call the smoothing routine nvfix = 0 nqint = 10 lrwork = 2*nv1 + nelt1 liwork = 8*nelt1 + 2*nv1 ALLOCATE (numfix(nvfix),rwork(lrwork),iwork(liwork)) ifail = 0 CALL d06caf(nv1,nelt1,nedge1,coor1,edge1,conn1,nvfix,numfix,itrace, & nqint,iwork,liwork,rwork,lrwork,ifail) DEALLOCATE (rwork,iwork,coorch,lined,lcomp,rate,nlcomp,crus) ! build the mesh of the 2nd domain ! Initialise boundary mesh inputs: ! the number of line and of the characteristic points of ! the boundary mesh READ (nin,*) nlines ALLOCATE (lined(4,nlines),lcomp(nlines),coorch(2,nlines), & rate(nlines),coor2(2,nvmax),edge2(3,nedmx)) ! Characteristic points of the boundary geometry READ (nin,*) coorch(1,1:nlines) READ (nin,*) coorch(2,1:nlines) ! The Lines of the boundary mesh READ (nin,*) (lined(1:4,j),rate(j),j=1,nlines) sdcrus = 0 DO i = 1, nlines IF (lined(4,i)<0) THEN sdcrus = sdcrus + lined(1,i) - 2 END IF END DO liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus ! Get max(LINED(1,:)) for computing LRWORK CALL f16dnf(nlines,lined,4,maxind,maxval) lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines ! The number of connected components to the boundary ! and their informations READ (nin,*) ncomp ALLOCATE (nlcomp(ncomp),crus(2,sdcrus),rwork(lrwork),iwork(liwork)) j = 1 DO i = 1, ncomp READ (nin,*) nlcomp(i) k = j + abs(nlcomp(i)) - 1 READ (nin,*) lcomp(j:k) j = k + 1 END DO itrace = 0 ! Call to the 2D boundary mesh generator ifail = 0 CALL d06baf(nlines,coorch,lined,fbnd,crus,sdcrus,rate,ncomp,nlcomp, & lcomp,nvmax,nedmx,nvb2,coor2,nedge2,edge2,itrace,ruser,iuser, & rwork,lrwork,iwork,liwork,ifail) DEALLOCATE (rwork,iwork,weight) ! mesh it using the advancing front method ! Initialise mesh control parameters itrace = 0 nvint = 0 lrwork = 12*nvmax + 30015 liwork = 8*nedge2 + 53*nvmax + 2*nvb2 + 10078 ALLOCATE (weight(nvint),rwork(lrwork),iwork(liwork), & conn2(3,2*nvmax+5)) ! Call to the 2D Advancing front mesh generator ifail = 0 CALL d06acf(nvb2,nvint,nvmax,nedge2,edge2,nv2,nelt2,coor2,conn2, & weight,itrace,rwork,lrwork,iwork,liwork,ifail) DEALLOCATE (rwork,iwork) ! Rotation of the 2nd domain mesh to produce ! the 3rd mesh domain ntrans = 1 lrwork = 12*ntrans ALLOCATE (rwork(lrwork),itype(ntrans),trans(6,ntrans),coor3(2,nv2), & edge3(3,nedge2),conn3(3,nelt2)) itype(1:ntrans) = (/ 3/) trans(1,1:ntrans) = (/ 6.0_nag_wp/) trans(2,1:ntrans) = (/ -1.0_nag_wp/) trans(3,1:ntrans) = (/ -90.0_nag_wp/) itrace = 0 ifail = 0 CALL d06daf(nv2,nedge2,nelt2,ntrans,itype,trans,coor2,edge2,conn2, & coor3,edge3,conn3,itrace,rwork,lrwork,ifail) DEALLOCATE (rwork) nv3 = nv2 nelt3 = nelt2 nedge3 = nedge2 nv4 = nv1 + nv2 nelt4 = nelt1 + nelt2 nedge4 = nedge1 + nedge2 liwork = 2*nv1 + 3*nv2 + nelt1 + nelt2 + nedge1 + nedge2 + 1024 ALLOCATE (iwork(liwork),coor4(2,nv4),edge4(3,nedge4),conn4(3,nelt4), & reft4(nelt4),reft1(nelt1),reft2(nelt2)) ! restitching of the mesh 1 and 2 to form the mesh 4 reft1(1:nelt1) = 1 reft2(1:nelt2) = 2 eps = 1.E-3_nag_wp itrace = 0 ifail = 0 CALL d06dbf(eps,nv1,nelt1,nedge1,coor1,edge1,conn1,reft1,nv2,nelt2, & nedge2,coor2,edge2,conn2,reft2,nv4,nelt4,nedge4,coor4,edge4, & conn4,reft4,itrace,iwork,liwork,ifail) DEALLOCATE (iwork) nv5 = nv4 + nv3 nelt5 = nelt4 + nelt3 nedge5 = nedge4 + nedge3 liwork = 2*nv4 + 3*nv3 + nelt4 + nelt3 + nedge4 + nedge3 + 1024 ALLOCATE (iwork(liwork),coor5(2,nv5),edge5(3,nedge5),conn5(3,nelt5), & reft5(nelt5),reft3(nelt3)) ! restitching of the mesh 3 and 4 to form the mesh 5 reft3(1:nelt3) = 3 itrace = 0 ifail = 0 CALL d06dbf(eps,nv4,nelt4,nedge4,coor4,edge4,conn4,reft4,nv3,nelt3, & nedge3,coor3,edge3,conn3,reft3,nv5,nelt5,nedge5,coor5,edge5, & conn5,reft5,itrace,iwork,liwork,ifail) READ (nin,*) pmesh SELECT CASE (pmesh) CASE ('N') WRITE (nout,*) 'The restitched mesh characteristics' WRITE (nout,99999) 'NV =', nv5 WRITE (nout,99999) 'NELT =', nelt5 WRITE (nout,99999) 'NEDGE =', nedge5 CASE ('Y') ! Output the mesh WRITE (nout,99998) nv5, nelt5, nedge5 DO i = 1, nv5 WRITE (nout,99997) coor5(1:2,i) END DO DO k = 1, nelt5 WRITE (nout,99996) conn5(1,k), conn5(2,k), conn5(3,k), & reft5(k) END DO DO k = 1, nedge5 WRITE (nout,99998) edge5(1:3,k) END DO CASE DEFAULT WRITE (nout,*) 'Problem with the printing option Y or N' END SELECT 99999 FORMAT (1X,A,I6) 99998 FORMAT (1X,3I10) 99997 FORMAT (2(2X,E13.6)) 99996 FORMAT (1X,4I10) END SUBROUTINE ex2 END PROGRAM d06dbfe