! D06BAF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d06bafe_mod ! D06BAF 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, xa, xb, y0 ! .. Executable Statements .. xa = ruser(1) xb = ruser(2) x0 = ruser(3) y0 = ruser(4) fbnd = 0.0_nag_wp Select Case (i) Case (1) ! line 1,2,3, and 4: ellipse centred in (X0,Y0) with ! XA and XB as coefficients fbnd = ((x-x0)/xa)**2 + ((y-y0)/xb)**2 - 1.0_nag_wp Case (2) ! line 7, lower arc on letter n, is a circle centred in (X0,Y0) ! with radius SQRT(RADIUS2) x0 = 0.5_nag_wp y0 = 6.25_nag_wp radius2 = 20.3125_nag_wp fbnd = (x-x0)**2 + (y-y0)**2 - radius2 Case (3) ! line 11, upper arc on letter n, is a circle centred in (X0,Y0) ! with radius SQRT(RADIUS2) x0 = 1.0_nag_wp y0 = 4.0_nag_wp radius2 = 9.0_nag_wp + (11.0_nag_wp-y0)**2 fbnd = (x-x0)**2 + (y-y0)**2 - radius2 Case (4) ! line 15, upper arc on letter a, is a circle centred in (X0,Y0) ! with radius SQRT(RADIUS2) touching point (5,11). x0 = 8.5_nag_wp y0 = 2.75_nag_wp radius2 = (x0-5.0_nag_wp)**2 + (11.0_nag_wp-y0)**2 fbnd = (x-x0)**2 + (y-y0)**2 - radius2 Case (5) ! line 25, lower arc on hat of 'a', is a circle centred in (X0,Y0) ! with radius SQRT(RADIUS2) touching point (11,10). x0 = 8.5_nag_wp y0 = 4.0_nag_wp radius2 = 2.5_nag_wp**2 + (10.0_nag_wp-y0)**2 fbnd = (x-x0)**2 + (y-y0)**2 - radius2 Case (6) ! lines 20, 21 and 22, belly of letter a, is an ellipse centered ! in (X0, Y0) with semi-axes 3.5 and 2.75. x0 = 8.5_nag_wp y0 = 5.75_nag_wp fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp Case (7) ! lines 43, 44 and 45, outer curve on bottom of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 3.5 and 2.5. x0 = 17.5_nag_wp y0 = 2.5_nag_wp fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.5_nag_wp)**2 - 1.0_nag_wp Case (8) ! lines 28, 29 and 30, inner curve on bottom of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 2.0 and 1.5. x0 = 17.5_nag_wp y0 = 2.5_nag_wp fbnd = ((x-x0)/2.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp Case (9) ! line 42, inner curve on lower middle of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 1.5 and 0.5. x0 = 17.5_nag_wp y0 = 5.5_nag_wp fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/0.5_nag_wp)**2 - 1.0_nag_wp Case (10) ! line 31, outer curve on lower middle of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 2.0 and 1.5. x0 = 17.5_nag_wp y0 = 5.5_nag_wp fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp Case (11) ! line 41, inner curve on upper middle of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 1.0 and 1.0. x0 = 17.0_nag_wp y0 = 5.5_nag_wp fbnd = ((x-x0)/1.0_nag_wp)**2 + ((y-y0)/1.0_nag_wp)**2 - 1.0_nag_wp Case (12) ! line 32, outer curve on upper middle of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 1.5 and 1.1573. x0 = 16.0_nag_wp y0 = 5.5_nag_wp fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/1.1573_nag_wp)**2 - & 1.0_nag_wp Case (13) ! lines 33, 33, 34, 39 and 40, upper portion of 'g', is an ellipse ! centered in (X0, Y0) with semi-axes 3.0 and 2.75. x0 = 17.0_nag_wp y0 = 9.25_nag_wp fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp End Select Return End Function fbnd End Module d06bafe_mod Program d06bafe ! D06BAF Example Main Program ! .. Use Statements .. Use nag_library, Only: d06abf, d06acf, d06baf, f16dnf, nag_wp Use d06bafe_mod, Only: fbnd, nin, nout ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: x0, xa, xb, xmax, xmin, y0, & ymax, ymin Integer :: i, ifail, itrace, j, k, liwork, & lrwork, maxind, maxval, ncomp, & nedge, nedmx, nelt, nlines, & npropa, nv, nvb, nvint, nvmax, & reftk, sdcrus Character (1) :: pmesh ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: coor(:,:), coorch(:,:), & crus(:,:), rate(:), rwork(:), & weight(:) Real (Kind=nag_wp) :: ruser(4) Integer, Allocatable :: conn(:,:), edge(:,:), iwork(:), & lcomp(:), lined(:,:), nlcomp(:) Integer :: iuser(1) ! .. Intrinsic Procedures .. Intrinsic :: abs ! .. Executable Statements .. Write (nout,*) 'D06BAF Example Program Results' Flush (nout) ! Skip heading in data file Read (nin,*) ! Initialise boundary mesh inputs: ! the number of line and of the characteristic points of ! the boundary mesh Read (nin,*) nlines, nvmax, nedmx Allocate (coor(2,nvmax),coorch(2,nlines),rate(nlines),edge(3,nedmx), & lcomp(nlines),lined(4,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 Allocate (crus(2,sdcrus),iwork(liwork),rwork(lrwork)) ! The ellipse boundary which envelops the NAg Logo ! the N, the A and the G Read (nin,*) coorch(1,1:nlines) Read (nin,*) coorch(2,1:nlines) Read (nin,*) crus(1,1:sdcrus) Read (nin,*) crus(2,1:sdcrus) ! The number of connected components to the boundary ! and their informations Read (nin,*) ncomp Allocate (nlcomp(ncomp)) 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 ! Data passed to the user-supplied function xmin = coorch(1,4) xmax = coorch(1,2) ymin = coorch(2,1) ymax = coorch(2,3) xa = (xmax-xmin)/2.0_nag_wp xb = (ymax-ymin)/2.0_nag_wp x0 = (xmin+xmax)/2.0_nag_wp y0 = (ymin+ymax)/2.0_nag_wp ruser(1:4) = (/xa,xb,x0,y0/) iuser(1) = 0 itrace = -1 Write (nout,*) Flush (nout) ! Call to the boundary mesh generator ifail = 0 Call d06baf(nlines,coorch,lined,fbnd,crus,sdcrus,rate,ncomp,nlcomp, & lcomp,nvmax,nedmx,nvb,coor,nedge,edge,itrace,ruser,iuser,rwork,lrwork, & iwork,liwork,ifail) Read (nin,*) pmesh Select Case (pmesh) Case ('N') Write (nout,*) 'Boundary mesh characteristics' Write (nout,99999) 'NVB =', nvb Write (nout,99999) 'NEDGE =', nedge Case ('Y') ! Output the mesh Write (nout,99998) nvb, nedge Do i = 1, nvb Write (nout,99997) i, coor(1:2,i) End Do Do i = 1, nedge Write (nout,99996) i, edge(1:3,i) End Do Flush (nout) Case Default Write (nout,*) 'Problem with the printing option Y or N' Go To 100 End Select Deallocate (rwork,iwork) ! Initialise mesh control parameters itrace = 0 npropa = 1 nvint = 0 lrwork = 12*nvmax + 15 liwork = 6*nedge + 32*nvmax + 2*nvb + 78 Allocate (weight(nvint),rwork(lrwork),iwork(liwork),conn(3,2*nvmax+5)) ! Call to the 2D Delaunay-Voronoi mesh generator ifail = 0 Call d06abf(nvb,nvint,nvmax,nedge,edge,nv,nelt,coor,conn,weight,npropa, & itrace,rwork,lrwork,iwork,liwork,ifail) Select Case (pmesh) Case ('N') Write (nout,*) 'Complete mesh (via the 2D Delaunay-Voronoi' Write (nout,*) 'mesh generator) characteristics' Write (nout,99999) 'NV =', nv Write (nout,99999) 'NELT =', nelt Case ('Y') ! Output the mesh Write (nout,99998) nv, nelt Do i = 1, nv Write (nout,99995) coor(1:2,i) End Do reftk = 0 Do k = 1, nelt Write (nout,99994) conn(1:3,k), reftk End Do Flush (nout) End Select Deallocate (rwork,iwork) lrwork = 12*nvmax + 30015 liwork = 8*nedge + 53*nvmax + 2*nvb + 10078 Allocate (rwork(lrwork),iwork(liwork)) ! Call to the 2D Advancing front mesh generator ifail = 0 Call d06acf(nvb,nvint,nvmax,nedge,edge,nv,nelt,coor,conn,weight,itrace, & rwork,lrwork,iwork,liwork,ifail) Select Case (pmesh) Case ('N') Write (nout,*) 'Complete mesh (via the 2D Advancing front mesh' Write (nout,*) 'generator) characteristics' Write (nout,99999) 'NV =', nv Write (nout,99999) 'NELT =', nelt Case ('Y') ! Output the mesh Write (nout,99998) nv, nelt Do i = 1, nv Write (nout,99995) coor(1:2,i) End Do reftk = 0 Do k = 1, nelt Write (nout,99994) conn(1:3,k), reftk End Do End Select 100 Continue 99999 Format (1X,A,I6) 99998 Format (1X,2I10) 99997 Format (2X,I4,2(2X,E13.6)) 99996 Format (1X,4I4) 99995 Format (2(2X,E13.6)) 99994 Format (1X,4I10) End Program d06bafe