! D06DBF Example Program Text
! Mark 28.3 Release. NAG Copyright 2022.
Module d06dbfe_mod
! D06DBF 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 :: meshout = 7, nin = 5, nout = 6, &
nvb1 = 19
Logical, Parameter, Public :: pmesh = .False.
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, y0
! .. Executable Statements ..
fbnd = 0.0_nag_wp
Select Case (i)
Case (1)
! inner circle
x0 = 0.0_nag_wp
y0 = 0.0_nag_wp
radius2 = 1.0_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
Case (2)
! outer circle
x0 = 0.0_nag_wp
y0 = 0.0_nag_wp
radius2 = 5.0_nag_wp
fbnd = (x-x0)**2 + (y-y0)**2 - radius2
End Select
Return
End Function fbnd
End Module d06dbfe_mod
Program d06dbfe
! D06DBF Example Main Program
! .. Use Statements ..
Use d06dbfe_mod, Only: nout
Use nag_library, Only: nag_wp
! .. Implicit None Statement ..
Implicit None
! .. Executable Statements ..
Write (nout,*) 'D06DBF Example Program Results'
Call ex1
Call ex2
Contains
Subroutine ex1
! .. Use Statements ..
Use d06dbfe_mod, Only: meshout, nout, nvb1, pmesh
Use nag_library, Only: d06daf, d06dbf, x01aaf
! .. Local Scalars ..
Real (Kind=nag_wp) :: ddx, ddy, dx, eps, pi2, pi4, r_i, &
r_j
Integer :: i, ifail, imax, itrace, itrans, j, &
jmax, jtrans, k, ktrans, l, liwork, &
lrwork, nedge1, nedge2, nedge3, &
nelt1, nelt2, nelt3, ntrans, nv1, &
nv2, nv3
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: coor1(:,:), coor2(:,:), coor3(:,:), &
rwork(:), trans(:,:)
Integer, Allocatable :: conn1(:,:), conn2(:,:), conn3(:,:), &
edge1(:,:), edge2(:,:), edge3(:,:), &
itype(:), iwork(:), reft1(:), &
reft2(:), reft3(:)
! .. Intrinsic Procedures ..
Intrinsic :: real, sin
! .. Executable Statements ..
Write (nout,*)
Write (nout,*) 'Example 1'
Write (nout,*)
imax = nvb1 + 1
jmax = imax
nedge1 = 2*(imax-1) + 2*(jmax-1)
nedge2 = nedge1
nedge3 = nedge1 + nedge2
ntrans = 1
lrwork = 12*ntrans
! Allocate for mesh : coordinates and connectivity of the 1st domain,
! 2nd translated domain and restitched domain.
nv1 = (nvb1+1)**2
nelt1 = 2*nvb1*nvb1
nv2 = nv1
nv3 = nv1 + nv2
nelt2 = nelt1
nelt3 = nelt1 + nelt2
liwork = 2*nv1 + 3*nv2 + nelt1 + nelt2 + nedge1 + nedge2 + 1024
Allocate (coor1(2,nv1),coor2(2,nv2),coor3(2,nv3),conn1(3,nelt1), &
conn2(3,nelt2),conn3(3,nelt3),reft1(nelt1),reft2(nelt2), &
reft3(nelt3),edge1(3,nedge1),edge2(3,nedge2),edge3(3,nedge3), &
itype(ntrans),trans(6,ntrans),rwork(lrwork),iwork(liwork))
! Set up interior mesh as small perturbations on regular grid
! with regularity on the boundary of square on [0,1]x[0,1].
dx = 1.0_nag_wp/real(nvb1,kind=nag_wp)
pi2 = x01aaf(dx) + x01aaf(dx)
pi4 = pi2 + pi2
k = 0
Do j = 1, jmax
ddy = real(j-1,kind=nag_wp)*dx
Do i = 1, imax
k = k + 1
ddx = real(i-1,kind=nag_wp)*dx
coor1(1,k) = ddx + dx*0.05_nag_wp*sin(pi4*ddx)*sin(pi4*ddy)
coor1(2,k) = ddy + dx*0.05_nag_wp*sin(pi2*ddx)*sin(pi2*ddy)
End Do
End Do
! Triangulate using skew-diagonals on grid squares
k = 0
l = 0
Do i = 1, nvb1
Do j = 1, nvb1
l = l + 1
k = k + 1
conn1(1,k) = l
conn1(2,k) = l + 1
conn1(3,k) = l + nvb1 + 2
k = k + 1
conn1(1,k) = l
conn1(2,k) = l + nvb1 + 2
conn1(3,k) = l + nvb1 + 1
End Do
l = l + 1
End Do
reft1(1:nelt1) = 1
reft2(1:nelt2) = 2
! Define the edges of the boundary
Do i = 1, nvb1
edge1(1,i) = i
edge1(2,i) = i + 1
End Do
Do i = 1, nvb1
edge1(1,nvb1+i) = i*imax
edge1(2,nvb1+i) = (i+1)*imax
End Do
Do i = 1, nvb1
edge1(1,2*nvb1+i) = imax*jmax - i + 1
edge1(2,2*nvb1+i) = imax*jmax - i
End Do
Do i = 1, nvb1
edge1(1,3*nvb1+i) = (jmax-i)*imax + 1
edge1(2,3*nvb1+i) = (jmax-i-1)*imax + 1
End Do
edge1(3,1:nedge1) = 1
If (pmesh) Then
! Print interior mesh of single square
Do i = 1, nelt1
Write (meshout,99997) coor1(1,conn1(1,i)), coor1(2,conn1(1,i))
Write (meshout,99997) coor1(1,conn1(2,i)), coor1(2,conn1(2,i))
Write (meshout,99997) coor1(1,conn1(3,i)), coor1(2,conn1(3,i))
Write (meshout,99997) coor1(1,conn1(1,i)), coor1(2,conn1(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
End If
Do ktrans = 1, 2
! Translation of the 1st domain to obtain the 2nd domain
! KTRANS = 1 leading to a domains (4x2) overlapping
! KTRANS = 2 leading to a domains partition
If (ktrans==1) Then
itrans = nvb1 - 4
jtrans = nvb1 - 2
Else
itrans = nvb1
jtrans = 0
End If
itype(1:ntrans) = (/1/)
r_i = real(itrans,kind=nag_wp)/real(nvb1,kind=nag_wp)
r_j = real(jtrans,kind=nag_wp)/real(nvb1,kind=nag_wp)
trans(1,1:ntrans) = (/r_i/)
trans(2,1:ntrans) = (/r_j/)
itrace = 0
ifail = 0
Call d06daf(nv2,nedge2,nelt2,ntrans,itype,trans,coor1,edge1,conn1, &
coor2,edge2,conn2,itrace,rwork,lrwork,ifail)
edge2(3,1:nedge2) = 2
! Call to the restitching driver
itrace = 0
eps = 1.E-2_nag_wp
ifail = 0
Call d06dbf(eps,nv1,nelt1,nedge1,coor1,edge1,conn1,reft1,nv2,nelt2, &
nedge2,coor2,edge2,conn2,reft2,nv3,nelt3,nedge3,coor3,edge3,conn3, &
reft3,itrace,iwork,liwork,ifail)
If (pmesh) Then
! Output the overlapping or partitioned mesh
Write (nout,99998) nv3, nelt3, nedge3
Do i = 1, nelt3
Write (meshout,99997) coor3(1,conn3(1,i)), coor3(2,conn3(1,i))
Write (meshout,99997) coor3(1,conn3(2,i)), coor3(2,conn3(2,i))
Write (meshout,99997) coor3(1,conn3(3,i)), coor3(2,conn3(3,i))
Write (meshout,99997) coor3(1,conn3(1,i)), coor3(2,conn3(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
Do k = 1, nelt3
Write (nout,99996) conn3(1:3,k), reft3(k)
End Do
Do k = 1, nedge3
Write (nout,99998) edge3(1:3,k)
End Do
Else
If (ktrans==1) Then
Write (nout,*) &
'The restitched mesh characteristics in the overlapping case'
Else
Write (nout,*) &
'The restitched mesh characteristics in the partition case'
End If
Write (nout,99999) 'nv =', nv3
Write (nout,99999) 'nelt =', nelt3
Write (nout,99999) 'nedge =', nedge3
End If
End Do
99999 Format (1X,A,I6)
99998 Format (1X,3I10)
99997 Format (2(2X,E13.6))
99996 Format (1X,4I10)
End Subroutine ex1
Subroutine ex2
! .. Use Statements ..
Use d06dbfe_mod, Only: fbnd, meshout, nin, nout, pmesh
Use nag_library, Only: d06abf, d06acf, d06baf, d06caf, d06daf, d06dbf, &
f16dnf
! .. Local Scalars ..
Real (Kind=nag_wp) :: eps
Integer :: i, ifail, itrace, j, k, liwork, &
lrwork, maxind, maxval, ncomp, &
nedge1, nedge2, nedge3, nedge4, &
nedge5, nedmx, nelt1, nelt2, nelt3, &
nelt4, nelt5, nlines, npropa, nqint, &
ntrans, nv1, nv2, nv3, nv4, nv5, &
nvb1, nvb2, nvfix, nvint, nvmax, &
sdcrus
! .. Local Arrays ..
Real (Kind=nag_wp), Allocatable :: coor1(:,:), coor2(:,:), coor3(:,:), &
coor4(:,:), coor5(:,:), coorch(:,:), &
crus(:,:), rate(:), rwork(:), &
trans(:,:), weight(:)
Real (Kind=nag_wp) :: ruser(1)
Integer, Allocatable :: conn1(:,:), conn2(:,:), conn3(:,:), &
conn4(:,:), conn5(:,:), edge1(:,:), &
edge2(:,:), edge3(:,:), edge4(:,:), &
edge5(:,:), itype(:), iwork(:), &
lcomp(:), lined(:,:), nlcomp(:), &
numfix(:), reft1(:), reft2(:), &
reft3(:), reft4(:), reft5(:)
Integer :: iuser(1)
! .. Intrinsic Procedures ..
Intrinsic :: abs
! .. Executable Statements ..
Write (nout,*)
Write (nout,*) 'Example 2'
Write (nout,*)
! Skip headings in data file
Read (nin,*)
Read (nin,*)
! Build and mesh two domains and rotate/translate second to create
! third.
! ----------------------------------------------
! 1st domain: Annulus with straight right edge.
! First two points are end of straight edge.
! The mesh for inner circle is defined by four NWSE points with mesh
! points on the quarter-circle between each pair (fbnd, i=1).
! The mesh for the incomplete outer circle is defined by three NWS
! points and the edge points, mesh points computed using fbnd,
! i=2.
! The number of lines (1+4+4).
Read (nin,*) nlines, nvmax, nedmx
Allocate (coor1(2,nvmax),edge1(3,nedmx),lined(4,nlines),lcomp(nlines), &
coorch(2,nlines),rate(nlines))
! Characteristic points of the boundary geometry.
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)
! Allocate workspace for d06baf.
! sdcrus = 0 in this case.
sdcrus = 0
Do i = 1, nlines
If (lined(4,i)<0) Then
sdcrus = sdcrus + lined(1,i) - 2
End If
End Do
! Get max(LINED(1,:)) for computing lrwork
Call f16dnf(nlines,lined,4,maxind,maxval)
liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus
lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines
Allocate (crus(2,sdcrus),rwork(lrwork),iwork(liwork))
! The number of connected components (outer circle/edge and inner
! circle)
Read (nin,*) ncomp
Allocate (nlcomp(ncomp))
! Read the lines comprising each connected component
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
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,nvb1,coor1,nedge1,edge1,itrace,ruser,iuser,rwork, &
lrwork,iwork,liwork,ifail)
Deallocate (rwork,iwork)
! Generate mesh using Delaunay-Voronoi method
! Initialize mesh control parameters and allocate workspace.
itrace = 0
npropa = 1
nvint = 0
lrwork = 12*nvmax + 15
liwork = 6*nedge1 + 32*nvmax + 2*nvb1 + 78
Allocate (weight(nvint),rwork(lrwork),iwork(liwork), &
conn1(3,2*nvmax+5))
! Call to the 2D Delaunay-Voronoi mesh generator
ifail = 0
Call d06abf(nvb1,nvint,nvmax,nedge1,edge1,nv1,nelt1,coor1,conn1, &
weight,npropa,itrace,rwork,lrwork,iwork,liwork,ifail)
Deallocate (rwork,iwork)
! Call the smoothing routine
nvfix = 0
nqint = 10
lrwork = 2*nv1 + nelt1
liwork = 8*nelt1 + 2*nv1
Allocate (numfix(nvfix),rwork(lrwork),iwork(liwork))
ifail = 0
Call d06caf(nv1,nelt1,nedge1,coor1,edge1,conn1,nvfix,numfix,itrace, &
nqint,iwork,liwork,rwork,lrwork,ifail)
Deallocate (rwork,iwork,coorch,lined,lcomp,rate,nlcomp,crus)
! ----------------------------------------------
! 2nd domain: a rectangle (4 by 2) abutting straight edge of 1st domain.
Read (nin,*) nlines
Allocate (lined(4,nlines),lcomp(nlines),coorch(2,nlines),rate(nlines), &
coor2(2,nvmax),edge2(3,nedmx))
! Characteristic points of the boundary geometry
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 again here.
sdcrus = 0
Do i = 1, nlines
If (lined(4,i)<0) Then
sdcrus = sdcrus + lined(1,i) - 2
End If
End Do
! Generate initial mesh using Delaunay-Voronoi method
liwork = 8*nlines + nvmax + 3*nedmx + 3*sdcrus
Call f16dnf(nlines,lined,4,maxind,maxval)
lrwork = 2*(nlines+sdcrus) + 2*maxval*nlines
Allocate (crus(2,sdcrus),rwork(lrwork),iwork(liwork))
! The number of connected components (1 rectangle)
Read (nin,*) ncomp
Allocate (nlcomp(ncomp))
! Four lines of rectangle
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
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,nvb2,coor2,nedge2,edge2,itrace,ruser,iuser,rwork, &
lrwork,iwork,liwork,ifail)
Deallocate (rwork,iwork,weight)
! Remesh 2nd domain using the advancing front method
! Initialize mesh control parameters
itrace = 0
nvint = 0
lrwork = 12*nvmax + 30015
liwork = 8*nedge2 + 53*nvmax + 2*nvb2 + 10078
Allocate (weight(nvint),rwork(lrwork),iwork(liwork), &
conn2(3,2*nvmax+5))
! Call to the 2D Advancing front mesh generator
ifail = 0
Call d06acf(nvb2,nvint,nvmax,nedge2,edge2,nv2,nelt2,coor2,conn2, &
weight,itrace,rwork,lrwork,iwork,liwork,ifail)
Deallocate (rwork,iwork)
! ----------------------------------------------
! 3rd domain: rotation and translation of the 2nd domain mesh
ntrans = 1
lrwork = 12*ntrans
Allocate (rwork(lrwork),itype(ntrans),trans(6,ntrans),coor3(2,nv2), &
edge3(3,nedge2),conn3(3,nelt2))
itype(1:ntrans) = (/3/)
trans(1,1:ntrans) = (/6.0_nag_wp/)
trans(2,1:ntrans) = (/-1.0_nag_wp/)
trans(3,1:ntrans) = (/-90.0_nag_wp/)
itrace = 0
ifail = 0
Call d06daf(nv2,nedge2,nelt2,ntrans,itype,trans,coor2,edge2,conn2, &
coor3,edge3,conn3,itrace,rwork,lrwork,ifail)
Deallocate (rwork)
! ----------------------------------------------
! Combine Meshes
nv3 = nv2
nelt3 = nelt2
nedge3 = nedge2
nv4 = nv1 + nv2
nelt4 = nelt1 + nelt2
nedge4 = nedge1 + nedge2
liwork = 2*nv1 + 3*nv2 + nelt1 + nelt2 + nedge1 + nedge2 + 1024
Allocate (iwork(liwork),coor4(2,nv4),edge4(3,nedge4),conn4(3,nelt4), &
reft4(nelt4),reft1(nelt1),reft2(nelt2))
! Restitching of the mesh 1 and 2 to form the mesh 4
reft1(1:nelt1) = 1
reft2(1:nelt2) = 2
eps = 1.E-3_nag_wp
itrace = 0
ifail = 0
Call d06dbf(eps,nv1,nelt1,nedge1,coor1,edge1,conn1,reft1,nv2,nelt2, &
nedge2,coor2,edge2,conn2,reft2,nv4,nelt4,nedge4,coor4,edge4,conn4, &
reft4,itrace,iwork,liwork,ifail)
Deallocate (iwork)
nv5 = nv4 + nv3
nelt5 = nelt4 + nelt3
nedge5 = nedge4 + nedge3
liwork = 2*nv4 + 3*nv3 + nelt4 + nelt3 + nedge4 + nedge3 + 1024
Allocate (iwork(liwork),coor5(2,nv5),edge5(3,nedge5),conn5(3,nelt5), &
reft5(nelt5),reft3(nelt3))
! Restitching of the mesh 3 and 4 to form the mesh 5
reft3(1:nelt3) = 3
itrace = 0
ifail = 0
Call d06dbf(eps,nv4,nelt4,nedge4,coor4,edge4,conn4,reft4,nv3,nelt3, &
nedge3,coor3,edge3,conn3,reft3,nv5,nelt5,nedge5,coor5,edge5,conn5, &
reft5,itrace,iwork,liwork,ifail)
If (pmesh) Then
! Output the mesh
Do i = 1, nelt5
Write (meshout,99997) coor5(1,conn5(1,i)), coor5(2,conn5(1,i))
Write (meshout,99997) coor5(1,conn5(2,i)), coor5(2,conn5(2,i))
Write (meshout,99997) coor5(1,conn5(3,i)), coor5(2,conn5(3,i))
Write (meshout,99997) coor5(1,conn5(1,i)), coor5(2,conn5(1,i))
Write (meshout,*)
End Do
Write (meshout,*)
Write (nout,99998) nv5, nelt5, nedge5
Do i = 1, nv5
Write (nout,99997) coor5(1:2,i)
End Do
Do k = 1, nelt5
Write (nout,99996) conn5(1,k), conn5(2,k), conn5(3,k), reft5(k)
End Do
Do k = 1, nedge5
Write (nout,99998) edge5(1:3,k)
End Do
Else
Write (nout,*) 'The restitched mesh characteristics'
Write (nout,99999) 'nv =', nv5
Write (nout,99999) 'nelt =', nelt5
Write (nout,99999) 'nedge =', nedge5
End If
99999 Format (1X,A,I6)
99998 Format (1X,3I10)
99997 Format (2(2X,E13.6))
99996 Format (1X,4I10)
End Subroutine ex2
End Program d06dbfe