! D06ACF Example Program Text ! Mark 23 Release. NAG Copyright 2011. 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) ! .. 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) :: c, radius, x0, x1, y0, y1 ! .. Intrinsic Functions .. 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 Functions .. 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 20 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 20 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