Program d06aafe
! D06AAF Example Program Text
! Mark 30.1 Release. NAG Copyright 2024.
! .. Use Statements ..
Use nag_library, Only: d06aaf, nag_wp, x01aaf
! .. Implicit None Statement ..
Implicit None
! .. Parameters ..
Integer, Parameter :: meshout = 7, nin = 5, nout = 6
! .. Local Scalars ..
Real (Kind=nag_wp) :: coef, pi2, power, r, theta, theta_i, &
x0, y0
Integer :: i, ifail, itrace, liwork, lrwork, &
nedge, nelt, nv, nvb, nvb1, nvb2, &
nvb3, nvmax
Logical :: smooth
Character (1) :: pmesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: bspace(:), coor(:,:), rwork(:)
Integer, Allocatable :: conn(:,:), edge(:,:), iwork(:)
! .. Intrinsic Procedures ..
Intrinsic :: cos, max, real, sin
! .. Executable Statements ..
Write (nout,*) 'D06AAF Example Program Results'
! Skip heading in data file
Read (nin,*)
! Reading of the geometry
! Coordinates of the boundary mesh vertices and
! edges references.
Read (nin,*) nvb1, nvb2, nvb3, nvmax
nvb = nvb1 + nvb2 + nvb3
nedge = nvb
lrwork = nvmax
liwork = 16*nvmax + 2*nedge + max(4*nvmax+2,nedge-14)
Allocate (bspace(nvb),coor(2,nvmax),rwork(lrwork),conn(3,2*(nvmax- &
1)),edge(3,nedge),iwork(liwork))
! Outer circle
pi2 = 2.0_nag_wp*x01aaf(theta)
theta = pi2/real(nvb1,kind=nag_wp)
r = 1.0_nag_wp
x0 = 0.0_nag_wp
y0 = 0.0_nag_wp
Do i = 1, nvb1
theta_i = theta*real(i,kind=nag_wp)
coor(1,i) = x0 + r*cos(theta_i)
coor(2,i) = y0 + r*sin(theta_i)
End Do
! Larger inner circle
theta = pi2/real(nvb2,kind=nag_wp)
r = 0.49_nag_wp
x0 = -0.5_nag_wp
y0 = 0.0_nag_wp
Do i = 1, nvb2
theta_i = theta*real(i,kind=nag_wp)
coor(1,nvb1+i) = x0 + r*cos(theta_i)
coor(2,nvb1+i) = y0 + r*sin(theta_i)
End Do
! Smaller inner circle
theta = pi2/real(nvb3,kind=nag_wp)
r = 0.15_nag_wp
x0 = -0.5_nag_wp
y0 = 0.65_nag_wp
Do i = 1, nvb3
theta_i = theta*real(i,kind=nag_wp)
coor(1,nvb1+nvb2+i) = x0 + r*cos(theta_i)
coor(2,nvb1+nvb2+i) = y0 + r*sin(theta_i)
End Do
! Boundary edges
Do i = 1, nedge
edge(1,i) = i
edge(2,i) = i + 1
edge(3,i) = 1
End Do
edge(2,nvb1) = 1
edge(2,nvb1+nvb2) = nvb1 + 1
edge(2,nvb) = nvb1 + nvb2 + 1
! Initialize mesh control parameters
bspace(1:nvb) = 0.05E0_nag_wp
smooth = .True.
itrace = 0
coef = 0.75E0_nag_wp
power = 0.25E0_nag_wp
! Call to the mesh generator
ifail = 0
Call d06aaf(nvb,nvmax,nedge,edge,nv,nelt,coor,conn,bspace,smooth,coef, &
power,itrace,rwork,lrwork,iwork,liwork,ifail)
Write (nout,*)
Read (nin,*) pmesh
Select Case (pmesh)
Case ('N')
Write (nout,99999) 'NV =', nv
Write (nout,99999) 'NELT =', nelt
Case ('Y')
! Output the mesh in a form suitable for printing
Write (meshout,*) '# D06ABF Example Program Mesh results'
Do i = 1, nelt
Write (meshout,99998) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,99998) coor(1,conn(2,i)), coor(2,conn(2,i))
Write (meshout,99998) coor(1,conn(3,i)), coor(2,conn(3,i))
Write (meshout,99998) coor(1,conn(1,i)), coor(2,conn(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
Case Default
Write (nout,*) 'Problem with the printing option Y or N'
End Select
99999 Format (1X,A,I6)
99998 Format (2(2X,E13.6))
End Program d06aafe