! D06ACF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
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
! .. Accessibility Statements ..
Private
Public :: fbnd
! .. Parameters ..
Integer, Parameter, Public :: 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 d06acfe_mod, Only: fbnd, nin, nout
Use nag_library, Only: d06acf, d06baf, f16dnf, nag_wp
! .. 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,*)
! Initialize 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 information
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)
! Initialize 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 (rounded to nearest 10) =', 10*((nv+5)/10)
Write (nout,99999) 'NELT (rounded to nearest 10) =', 10*((nelt+5)/10)
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