Program d06abfe
! D06ABF Example Program Text
! Mark 26.2 Release. NAG Copyright 2017.
! .. Use Statements ..
Use nag_library, Only: d06abf, f06epf, nag_wp, x01aaf
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: meshout = 7, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: c, d, dnvint, r, s, theta, theta_i
Integer :: i, i1, ic, ifail, itrace, j, liwork, &
lrwork, nearest, nedge, nelt, &
nelt_near, npropa, nrae, nv, nvb, &
nvint, nvmax, nv_near
Character (1) :: pmesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: coor(:,:), rwork(:), weight(:)
Real (Kind=nag_wp) :: t(2)
Integer, Allocatable :: conn(:,:), edge(:,:), iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: cos, real, sin, tan
! .. Executable Statements ..
Write (nout,*) 'D06ABF Example Program Results'
! Skip heading in data file
Read (nin,*)
! Reading of the geometry:
! number of points on RAE aerofoil data;
! number of points on circular boundary;
! maximum number of vertices.
Read (nin,*) nrae, nvint, nvmax
Read (nin,*) pmesh
nvb = nvint + 2*nrae
nedge = nvb
lrwork = 12*nvmax + 15
liwork = 6*nedge + 32*nvmax + 2*nvb + 78
Allocate (coor(2,nvmax),rwork(lrwork),weight(nvint),conn(3,2*nvmax+5), &
edge(3,nedge),iwork(liwork))
! Circular outer boundary, radius 3 and centre (1,0)
theta = 2.0_nag_wp*x01aaf(r)/real(nvint,kind=nag_wp)
r = 3.0_nag_wp
t(1) = 1.0_nag_wp
Do i = 1, nvint
theta_i = theta*real(i-1,kind=nag_wp)
coor(1,i) = r*cos(theta_i) + t(1)
coor(2,i) = r*sin(theta_i)
End Do
! Read data for aerofoil RAE 2822
Do i = 1, nrae
Read (nin,*) i1, coor(1,nvint+i), coor(2,nvint+i)
End Do
! Transform RAE 2822 for secondary foil
theta = x01aaf(theta)/12.0_nag_wp
c = cos(theta)
s = sin(theta)
ic = nvint + nrae
! Copy and rotate coordinates by theta = pi/12
coor(1:2,ic+1:ic+nrae) = coor(1:2,nvint+1:ic)
Call f06epf(nrae,coor(1,ic+1),2,coor(2,ic+1),2,c,s)
! Reduce by 0.4 and translate to distance 0.25 from intercept at (0.75,0)
d = 0.4_nag_wp
t(1) = 0.75_nag_wp + 0.25_nag_wp*c
t(2) = -0.25_nag_wp*s
Do i = 1, nrae
coor(1:2,ic+i) = d*coor(1:2,ic+i) + t(1:2)
End Do
! Boundary edges
Do i = 1, nedge
edge(1,i) = i
edge(2,i) = i + 1
edge(3,i) = 0
End Do
! Tie up end of three boundary edges
edge(2,nvint) = 1
edge(2,nvint+nrae) = nvint + 1
edge(2,nedge) = nvint + nrae + 1
! Initialize mesh control parameters
itrace = 0
! Generation of interior vertices on the
! RAE airfoil's wake
dnvint = 2.5E0_nag_wp/real(nvint+1,kind=nag_wp)
Do i = 1, nvint
i1 = nvb + i
coor(1,i1) = 1.38E0_nag_wp + real(i,kind=nag_wp)*dnvint
coor(2,i1) = -tan(theta)*(coor(1,i1)-0.75_nag_wp)
End Do
weight(1:nvint) = 0.01E0_nag_wp
Write (nout,*)
! Loop on the propagation coef
pcoef: Do j = 1, 4
nearest = 250
Select Case (j)
Case (1)
npropa = -5
Case (2)
npropa = -1
Case (3)
npropa = 1
Case Default
npropa = 5
End Select
! 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)
Write (nout,99999) 'Mesh characteristics with NPROPA =', npropa
nv_near = ((nv+nearest/2)/nearest)*nearest
nelt_near = ((nelt+nearest/2)/nearest)*nearest
Write (nout,99998) 'NV ', nv_near, nearest
Write (nout,99998) 'NELT', nelt_near, nearest
If (pmesh=='Y') Then
! Output the mesh in a form suitable for printing
If (j==1) Then
Write (meshout,*) '# D06ABF Example Program Mesh results'
End If
Write (meshout,99999) '# Mesh line segments for NPROPA =', npropa
Do i = 1, nelt
Write (meshout,99997) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,99997) coor(1,conn(2,i)), coor(2,conn(2,i))
Write (meshout,99997) coor(1,conn(3,i)), coor(2,conn(3,i))
Write (meshout,99997) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
End If
End Do pcoef
99999 Format (1X,A,I6)
99998 Format (1X,A5,' = ',I10,' to the nearest ',I3)
99997 Format (2(2X,E13.6))
End Program d06abfe