! D06DBF Example Program Text ! Mark 24 Release. NAG Copyright 2012. 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) ! .. 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 ! .. 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 Procedures .. 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 ! .. 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 Procedures .. 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