! D06BAF Example Program Text
! Mark 26.1 Release. NAG Copyright 2016.
Module d06bafe_mod
! D06BAF 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) :: radius2, x0, xa, xb, y0
! .. Executable Statements ..
xa = ruser(1)
xb = ruser(2)
x0 = ruser(3)
y0 = ruser(4)
fbnd = 0.0_nag_wp
Select Case (i)
Case (1)
! line 1,2,3, and 4: ellipse centred in (X0,Y0) with
! XA and XB as coefficients
fbnd = ((x-x0)/xa)**2 + ((y-y0)/xb)**2 - 1.0_nag_wp
Case (2)
! line 7, lower arc on letter n, is a circle centred in (X0,Y0)
! with radius SQRT(RADIUS2)
x0 = 0.5_nag_wp
y0 = 6.25_nag_wp
radius2 = 20.3125_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (3)
! line 11, upper arc on letter n, is a circle centred in (X0,Y0)
! with radius SQRT(RADIUS2)
x0 = 1.0_nag_wp
y0 = 4.0_nag_wp
radius2 = 9.0_nag_wp + (11.0_nag_wp-y0)**2
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (4)
! line 15, upper arc on letter a, is a circle centred in (X0,Y0)
! with radius SQRT(RADIUS2) touching point (5,11).
x0 = 8.5_nag_wp
y0 = 2.75_nag_wp
radius2 = (x0-5.0_nag_wp)**2 + (11.0_nag_wp-y0)**2
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (5)
! line 25, lower arc on hat of 'a', is a circle centred in (X0,Y0)
! with radius SQRT(RADIUS2) touching point (11,10).
x0 = 8.5_nag_wp
y0 = 4.0_nag_wp
radius2 = 2.5_nag_wp**2 + (10.0_nag_wp-y0)**2
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (6)
! lines 20, 21 and 22, belly of letter a, is an ellipse centered
! in (X0, Y0) with semi-axes 3.5 and 2.75.
x0 = 8.5_nag_wp
y0 = 5.75_nag_wp
fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp
Case (7)
! lines 43, 44 and 45, outer curve on bottom of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 3.5 and 2.5.
x0 = 17.5_nag_wp
y0 = 2.5_nag_wp
fbnd = ((x-x0)/3.5_nag_wp)**2 + ((y-y0)/2.5_nag_wp)**2 - 1.0_nag_wp
Case (8)
! lines 28, 29 and 30, inner curve on bottom of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 2.0 and 1.5.
x0 = 17.5_nag_wp
y0 = 2.5_nag_wp
fbnd = ((x-x0)/2.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp
Case (9)
! line 42, inner curve on lower middle of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 1.5 and 0.5.
x0 = 17.5_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/0.5_nag_wp)**2 - 1.0_nag_wp
Case (10)
! line 31, outer curve on lower middle of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 2.0 and 1.5.
x0 = 17.5_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/1.5_nag_wp)**2 - 1.0_nag_wp
Case (11)
! line 41, inner curve on upper middle of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 1.0 and 1.0.
x0 = 17.0_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/1.0_nag_wp)**2 + ((y-y0)/1.0_nag_wp)**2 - 1.0_nag_wp
Case (12)
! line 32, outer curve on upper middle of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 1.5 and 1.1573.
x0 = 16.0_nag_wp
y0 = 5.5_nag_wp
fbnd = ((x-x0)/1.5_nag_wp)**2 + ((y-y0)/1.1573_nag_wp)**2 - &
1.0_nag_wp
Case (13)
! lines 33, 33, 34, 39 and 40, upper portion of 'g', is an ellipse
! centered in (X0, Y0) with semi-axes 3.0 and 2.75.
x0 = 17.0_nag_wp
y0 = 9.25_nag_wp
fbnd = ((x-x0)/3.0_nag_wp)**2 + ((y-y0)/2.75_nag_wp)**2 - 1.0_nag_wp
End Select
Return
End Function fbnd
End Module d06bafe_mod
Program d06bafe
! D06BAF Example Main Program
! .. Use Statements ..
Use d06bafe_mod, Only: fbnd, nin, nout
Use nag_library, Only: d06abf, d06acf, d06baf, f16dnf, nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Local Scalars ..
Real (Kind=nag_wp) :: x0, xa, xb, xmax, xmin, y0, ymax, &
ymin
Integer :: i, ifail, itrace, j, k, liwork, &
lrwork, maxind, maxval, ncomp, &
nedge, nedmx, nelt, nlines, npropa, &
nv, nvb, nvint, nvmax, reftk, sdcrus
Character (1) :: pmesh
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: coor(:,:), coorch(:,:), crus(:,:), &
rate(:), rwork(:), weight(:)
Real (Kind=nag_wp) :: ruser(4)
Integer, Allocatable :: conn(:,:), edge(:,:), iwork(:), &
lcomp(:), lined(:,:), nlcomp(:)
Integer :: iuser(1)
! .. Intrinsic Procedures ..
Intrinsic :: abs
! .. Executable Statements ..
Write (nout,*) 'D06BAF Example Program Results'
Flush (nout)
! 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))
! 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
Allocate (crus(2,sdcrus),iwork(liwork),rwork(lrwork))
! The ellipse boundary which envelops the NAG Logo
! the N, the A and the G
Read (nin,*) coorch(1,1:nlines)
Read (nin,*) coorch(2,1:nlines)
Read (nin,*) crus(1,1:sdcrus)
Read (nin,*) crus(2,1:sdcrus)
! The number of connected components to the boundary
! and their information
Read (nin,*) ncomp
Allocate (nlcomp(ncomp))
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
xmin = coorch(1,4)
xmax = coorch(1,2)
ymin = coorch(2,1)
ymax = coorch(2,3)
xa = (xmax-xmin)/2.0_nag_wp
xb = (ymax-ymin)/2.0_nag_wp
x0 = (xmin+xmax)/2.0_nag_wp
y0 = (ymin+ymax)/2.0_nag_wp
ruser(1:4) = (/xa,xb,x0,y0/)
iuser(1) = 0
itrace = -1
Write (nout,*)
Flush (nout)
! Call to the 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)
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
Flush (nout)
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
npropa = 1
nvint = 0
lrwork = 12*nvmax + 15
liwork = 6*nedge + 32*nvmax + 2*nvb + 78
Allocate (weight(nvint),rwork(lrwork),iwork(liwork),conn(3,2*nvmax+5))
! 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)
Select Case (pmesh)
Case ('N')
Write (nout,*) 'Complete mesh (via the 2D Delaunay-Voronoi'
Write (nout,*) 'mesh generator) 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
Flush (nout)
End Select
Deallocate (rwork,iwork)
lrwork = 12*nvmax + 30015
liwork = 8*nedge + 53*nvmax + 2*nvb + 10078
Allocate (rwork(lrwork),iwork(liwork))
! 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 (via the 2D Advancing front mesh'
Write (nout,*) 'generator) 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 d06bafe