! D06ACF Example Program Text ! Mark 24 Release. NAG Copyright 2012. Module d06acfe_mod ! D06ACF 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) :: c, radius, x0, x1, y0, y1 ! .. Intrinsic Procedures .. Intrinsic :: sqrt ! .. Executable Statements .. fbnd = 0.0_nag_wp Select Case (i) Case (1) ! upper NACA0012 wing beginning at the origin c = 1.008930411365_nag_wp fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*x)-0.126_nag_wp*c*x- & 0.3516_nag_wp*(c*x)**2+0.2843_nag_wp*(c*x)**3- & 0.1015_nag_wp*(c*x)**4) - c*y Case (2) ! lower NACA0012 wing beginning at the origin c = 1.008930411365_nag_wp fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*x)-0.126_nag_wp*c*x- & 0.3516_nag_wp*(c*x)**2+0.2843_nag_wp*(c*x)**3- & 0.1015_nag_wp*(c*x)**4) + c*y Case (3) x0 = ruser(1) y0 = ruser(2) radius = ruser(3) fbnd = (x-x0)**2 + (y-y0)**2 - radius**2 Case (4) ! upper NACA0012 wing beginning at (X1;Y1) c = 1.008930411365_nag_wp x1 = ruser(4) y1 = ruser(5) fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*(x- & x1))-0.126_nag_wp*c*(x-x1)-0.3516_nag_wp*(c*(x- & x1))**2+0.2843_nag_wp*(c*(x-x1))**3-0.1015_nag_wp*(c*(x-x1))**4) - & c*(y-y1) Case (5) ! lower NACA0012 wing beginning at (X1;Y1) c = 1.008930411365_nag_wp x1 = ruser(4) y1 = ruser(5) fbnd = 0.6_nag_wp*(0.2969_nag_wp*sqrt(c*(x- & x1))-0.126_nag_wp*c*(x-x1)-0.3516_nag_wp*(c*(x- & x1))**2+0.2843_nag_wp*(c*(x-x1))**3-0.1015_nag_wp*(c*(x-x1))**4) + & c*(y-y1) End Select Return End Function fbnd End Module d06acfe_mod Program d06acfe ! D06ACF Example Main Program ! .. Use Statements .. Use nag_library, Only: d06acf, d06baf, f16dnf, nag_wp Use d06acfe_mod, Only: fbnd, nin, nout ! .. Implicit None Statement .. Implicit None ! .. Local Scalars .. Real (Kind=nag_wp) :: dnvint, radius, x0, x1, y0, y1 Integer :: i, ifail, itrace, j, k, liwork, & lrwork, maxind, maxval, ncomp, & nedge, nedmx, nelt, nlines, nv, & nvb, nvint, nvint2, nvmax, & reftk, sdcrus Character (1) :: pmesh ! .. Local Arrays .. Real (Kind=nag_wp), Allocatable :: coor(:,:), coorch(:,:), & crus(:,:), rate(:), rwork(:), & weight(:) Real (Kind=nag_wp) :: ruser(5) Integer, Allocatable :: conn(:,:), edge(:,:), iwork(:), & lcomp(:), lined(:,:), nlcomp(:) Integer :: iuser(1) ! .. Intrinsic Procedures .. Intrinsic :: abs, real ! .. Executable Statements .. Write (nout,*) 'D06ACF Example Program Results' ! 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)) 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 (crus(2,sdcrus),nlcomp(ncomp),iwork(liwork),rwork(lrwork)) 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 x0 = 1.5_nag_wp y0 = 0.0_nag_wp radius = 4.5_nag_wp x1 = 0.8_nag_wp y1 = -0.3_nag_wp ruser(1:5) = (/x0,y0,radius,x1,y1/) iuser(1) = 0 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,nvb,coor,nedge,edge,itrace,ruser,iuser,rwork,lrwork, & iwork,liwork,ifail) Write (nout,*) 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 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 ! Generation of interior vertices ! for the wake of the first NACA nvint = 40 lrwork = 12*nvmax + 30015 liwork = 8*nedge + 53*nvmax + 2*nvb + 10078 Allocate (weight(nvint),rwork(lrwork),conn(3,2*nvmax+5),iwork(liwork)) nvint2 = 20 dnvint = 5.0_nag_wp/real(nvint2+1,kind=nag_wp) Do i = 1, nvint2 reftk = nvb + i coor(1,reftk) = 1.0_nag_wp + real(i,kind=nag_wp)*dnvint coor(2,reftk) = 0.0_nag_wp End Do weight(1:nvint2) = 0.05_nag_wp ! for the wake of the second one dnvint = 4.19_nag_wp/real(nvint2+1,kind=nag_wp) Do i = nvint2 + 1, nvint reftk = nvb + i coor(1,reftk) = 1.8_nag_wp + real(i-nvint2,kind=nag_wp)*dnvint coor(2,reftk) = -0.3_nag_wp End Do weight((nvint2+1):nvint) = 0.05_nag_wp ! 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 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 d06acfe