纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之a

!*************************************************************
!*************************************************************
! 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之a
!*************************************************************
! SINOMACH
!
! 纳维斯托克斯2D-H形领域有限元法估算解决（N-S）方程(abc部)之a
!
! 周勇 20221125
!**************************************************************
!c $$$ Run, Output:
!**************************************************************
!**************************************************************
Program main

!*********************************************************************72
!
!c MAIN is the main program for HCELL.
!
! Discussion:
!
! HCELL solves the Navier-Stokes equations in an H-cell region.
!
! The fluid flow problem is formulated in terms of
! primitive variables - u,v, and p.
!
! u_t - laplacian u + (u.grad)u + grad p = f
! div u = 0
!
! Boundary conditions: (u,v)=(0,0) on top
! (u,v)=0 on left, right and bottom
!
! This version uses finite element techniques
! with piecewise linear functions on triangles to approximate
! the pressure and quadratics on triangles for the velocity
! (Taylor-Hood element), isoparametric element
!
! Input files:
!
! UP000.TXT contains the initial values of the solution coefficients.
!
! Local parameters:
!
! double precision AREA(NELEMN), the area of each element.
!
! integer BC_TYPE selects the boundary conditions, by controlling the value of ALPHA.
! Legal values are 1, for a step function, 2 for a "hat" function, 3 for a sinusoid.
!
! character ( len = * ) FILE_NAME is the name of the input file containing the
! initial values of the solution, and is also used as a template for the
! output files that store the solution at each timestep.
!
! integer INDX(MAXND,NUK), lists the indices of the U, V, and P
! variables associated with the node. It would be useful if the
! index associated with pressure also indicated whether there was
! no pressure variable associated with the node, or that it was
! prescribed. This could be done by assigning INDX(NODE,3) = 0
! for the midside nodes of the 6 node quadratic elements.
!
! integer MAXEL is an overestimate of the number of elements.
!
! integer MAXND is an overestimate of the number of nodes.
!
! integer MAXUN is an overestimate of the number of unknowns.
!
! integer MINUN is an estimate of the necessary bandwidth of the system matrix.
!
! integer MX counts the nuber of columns of nodes.
!
! integer MY counts the number of rows of nodes.
!
! integer NBAND, the bandwidth for the finite element matrix.
!
! integer NELEMN, the number of elements.
!
! integer NEQN, the total number of unknowns.
!
! integer NLBAND, NUBAND, the lower and upper half bandwidths
! for the finite element matrix.
!
! integer NNODES, the number of nodes per element.
!
! integer NODE(MAXEL,NNODES), the nodes that make up each element.
!
! integer NP, the number of nodes.
!
! integer NQUAD, the number of quadrature points used in assembly.
! (This is 3)
!
! integer NUK, the maximum number of unknowns associated with one node.
! (This is 3)
!
! integer NX counts, not quite all the elements in the X direction, but the number
! of elements plus 1.
!
! integer NX1, NX2, NX3, count the elements in the X direction in the three subregions,
!
! integer NY counts, not quite all the elements in the Y direction, but the number
! of elements plus 1.
!
! integer NY1, NY2, NY3, count the elements in the Y direction in the three subregions.
!
! integer REGION_DENSITY_X(3), REGION_DENSITY_Y(3), specifies the
! density of elements in the two coordinate directions.
!
! double precision REGION_X(4), REGION_Y(4), the coordinates of
! breakpoints that define 9 logical subrectangles.
!
! double precision XC(NP), YC(NP), the coordinates of the nodes.
!
! double precision XM(MAXEL,3), YM(MAXEL,3), the coordinates
! of quadrature points in each element.
!
Implicit Double Precision (A-H, O-Z)
!
! Set the "master parameters".
!
Parameter (nx1=45)
Parameter (nx2=15)
Parameter (nx3=45)
Parameter (ny1=5)
Parameter (ny2=1)
Parameter (ny3=5)
!
! Set parameters that depend on the master parameters.
!
Parameter (nx=nx1+nx2+nx3+1)
Parameter (ny=ny1+ny2+ny3+1)
Parameter (mx=2*nx-1)
Parameter (my=2*ny-1)
Parameter (maxel=2*(nx-1)*(ny-1))
Parameter (maxnd=mx*my)
Parameter (maxun=2*mx*my+nx*ny)
Parameter (minun=27*ny)
Parameter (nnodes=6)
Parameter (nuk=3)
Parameter (nquad=3)

Double Precision a(minun, maxun)
Double Precision area(maxel)
Integer bc_type
Double Precision f(maxun)
Double Precision g(maxun)
Integer indx(maxnd, nuk)
Integer ipivot(maxun)
Integer node(maxel, nnodes)
Character *80 node_file_name
Logical node_mask(maxnd)
Character *20 p_file_name
Integer region_density_x(3)
Integer region_density_y(3)
Double Precision region_x(4)
Double Precision region_y(4)
Double Precision reynld
Character *80 title
Double Precision uold(maxun)
Character *20 uv_file_name
Double Precision xc(maxnd)
Double Precision xm(maxel, nquad)
Double Precision yc(maxnd)
Double Precision ym(maxel, nquad)

Write (, '(a)') ' '
Write (, '(a)') 'HCELL:'
Write (, '(a)') ' FORTRAN95,2013&18 version'
Write (, '(a)') ' Solve the Navier Stokes fluid flow'
Write (, '(a)') ' equations in an H-shaped region,'
Write (, '(a)') ' using finite elements.'

bc_type = 1
p_file_name = 'p000.txt'
uv_file_name = 'uv000.txt'
!
! The density of elements in each region is assumed to be some multiple of the .
! For now, we will just use the densities.
!
region_density_x(1) = nx1
region_density_x(2) = nx2
region_density_x(3) = nx3

region_density_y(1) = ny1
region_density_y(2) = ny2
region_density_y(3) = ny3

Write (, '(a)') ' '
Write (, '(a)') ' The X direction is divided into three'
Write (, '(a)') ' regions, with element densities:'
Write (, '(a)') ' '
Write (, '(4x,3i6)')(region_density_x(i), i=1, 3)
Write (, '(a)') ' '
Write (, '(a,i6)') ' Corresponding NX = ', nx
Write (, '(a)') ' '
Write (, '(a)') ' The Y direction is divided into three'
Write (, '(a)') ' regions, with element densities:'
Write (, '(a)') ' '
Write (, '(4x,3i6)')(region_density_y(i), i=1, 3)
Write (, '(a)') ' '
Write (, '(a,i6)') ' Corresponding NY = ', ny
!
! Define the breakpoints that divide the region into 9 logical blocks.
!
region_x(1) = 0.0D+00
region_x(2) = 45.0D+00
region_x(3) = 60.0D+00
region_x(4) = 105.0D+00

region_y(1) = 0.0D+00
region_y(2) = 5.0D+00
region_y(3) = 6.0D+00
region_y(4) = 11.0D+00

Write (, '(a)') ' '
Write (, '(a)') ' The X subregions are demarked by 4 values:'
Write (, '(a)') ' '
Write (, '(4x,4f10.4)')(region_x(i), i=1, 4)
Write (, '(a)') ' '
Write (, '(a)') ' The Y subregions are demarked by 4 values:'
Write (, '(a)') ' '
Write (, '(4x,4f10.4)')(region_y(i), i=1, 4)

reynld = 1.0D0
mrow1 = minun
nsim = 3
nsteps = 10
tolns = 1.0D-6
tolopt = 1.0D-6
pi = 4.0D0*datan(1.0D0)

Write (, '(a)') ' '
Write (, '(a,i6)') ' Maximum number of nodes = ', maxnd
Write (, '(a,i6)') ' Maximum number of elements = ', maxel
Write (, '(a,i6)') ' Maximum number of unknowns = ', maxun
!
! SETGRD constructs grid, numbers unknowns, calculates areas,
! and points for midpoint quadrature rule, bandwidth and neqn
!
Call setgrd(xc, yc, area, xm, ym, region_density_x, region_density_y,&
region_x, region_y, node, indx, nlband, nuband, nband, nelemn, np,&
nnodes, nuk, nquad, neqn, maxnd, maxel)

Write (, '(a)') ' '
Write (, '(a,i6)') ' Number of nodes = ', np
Write (, '(a,i6)') ' Number of elements = ', nelemn
Write (, '(a,i6)') ' Number of unknowns = ', neqn
!
! Write the coordinates of the nodes of 6-node triangles to a file.
!
Call xy6_write('xy6.txt', np, xc, yc)
!
! Write the coordinates of the nodes of 3-node triangles to a file.
!
Call xy3_write('xy3.txt', maxnd, np, indx, xc, yc)
!
! Write the element node matrix to a file.
!
Call element_node_write('element_node.txt', maxel, nelemn, node)
!
! Make a plot of the nodes.
!
If (.True.) Then

Do i = 1, np
  If (i<=100) Then
    node_mask(i) = .True.
  Else
    node_mask(i) = .False.
  End If
End Do

title = 'H Cell Nodes'
node_file_name = 'hcell_nodes.eps'

Call node_eps(node_file_name, np, node_mask, xc, yc, title)

End If

Stop
nrow1 = nlband + nlband + nuband + 1
ncol1 = neqn

Do i = 1, neqn
f(i) = 0.0D0
End Do
!
! Read the initial value of the solution from a file.
!
If (.False.) Then

! call uv_read ( uv_file_name, neqn, uold )
! call p_read ( p_file_name, neqn, pold )

deltat = 0.0002D0
rdel = 1.0D0/deltat

Else
!
! Or assume we are doing a steady state calculation.
!
Do i = 1, neqn
uold(i) = 0.0D+00
End Do

deltat = 0.0D+00
rdel = 0.0D+00

End If
!
! Carry out the time iteration.
!
Do iter = 1, 1

If (bc_type==1) Then

  If (iter<=250) Then
    alpha = 5.0D0
  Else
    alpha = 1.0D0
  End If

Else If (bc_type==2) Then

  If (iter<=250) Then
    alpha = 80.0D0*dble(iter)*deltat + 1.0D0
  Else
    alpha = -80.0D0*dble(iter)*deltat + 9.0D0
  End If

Else If (bc_type==3) Then

  alpha = 2.0D0*sin(dble(iter)*0.01*pi)

End If

Do i = 1, neqn
  g(i) = f(i)
End Do

Do i = 1, neqn
  f(i) = 0.0D0
End Do

Call nstoke(xc, yc, area, xm, ym, a, f, g, uold, reynld, &
tolns, xlngth, ylngth, node, indx, ipivot, mrow1, nlband,&
nuband, nband, nrow1, ncol1, nelemn, np, nnodes, nuk, nquad,&
neqn, nsteps, nsim, maxnd, maxel, rdel, alpha)

!
! Save u=(gx,gy) to 'ue.dat' for 't.f'
!
Do i = 1, neqn
uold(i) = f(i)
End Do
!
! Increment the file name, and
! save u=(gx,gy) to 'up???.dat'
!
Call file_name_inc(uv_file_name)

Call uv_write(uv_file_name, maxnd, np, neqn, f, indx, xc, yc)

Call file_name_inc(p_file_name)

Call p_write(p_file_name, maxnd, np, neqn, f, indx, xc, yc)

End Do

Write (, '(a)') ' '
Write (, '(a)') 'HCELL:'
Write (*, '(a)') ' Normal end of execution.'

Stop
End Program main
Function ch_is_digit(c)

!*********************************************************************72
!
!c CH_IS_DIGIT returns .TRUE. if a character is a decimal digit.
!
Implicit None

Character c
Logical ch_is_digit

If (lge(c,'0') .And. lle(c,'9')) Then
ch_is_digit = .True.
Else
ch_is_digit = .False.
End If

Return
End Function ch_is_digit
Subroutine ch_to_digit(c, digit)

!*********************************************************************72
!
!c CH_TO_DIGIT returns the integer value of a 10 digit.
!
Implicit None

Character c
Integer digit

If (lge(c,'0') .And. lle(c,'9')) Then

digit = ichar(c) - 48

Else If (c==' ') Then

digit = 0

Else

digit = -1

End If

Return
End Subroutine ch_to_digit
Subroutine daxpy(n, da, dx, incx, dy, incy)

!*********************************************************************72
!
Double Precision dx(1), dy(1), da
Integer i, incx, incy, ix, iy, m, mp1, n

If (n<=0) Return
If (da==0.0D0) Return
If (incx==1 .And. incy==1) Goto 20
!
ix = 1
iy = 1
If (incx<0) ix = (-n+1)*incx + 1
If (incy<0) iy = (-n+1)*incy + 1
Do i = 1, n
dy(iy) = dy(iy) + da*dx(ix)
ix = ix + incx
iy = iy + incy
End Do
Return
!
20 m = mod(n, 4)
If (m==0) Goto 40
Do i = 1, m
dy(i) = dy(i) + da*dx(i)
End Do
If (n<4) Return
40 mp1 = m + 1
Do i = mp1, n, 4
dy(i) = dy(i) + da*dx(i)
dy(i+1) = dy(i+1) + da*dx(i+1)
dy(i+2) = dy(i+2) + da*dx(i+2)
dy(i+3) = dy(i+3) + da*dx(i+3)
End Do
Return
End Subroutine daxpy
Double Precision Function ddot(n, dx, incx, dy, incy)

!*********************************************************************72

Double Precision dx(1), dy(1), dtemp
Integer i, incx, incy, ix, iy, m, mp1, n

ddot = 0.0D0
dtemp = 0.0D0
If (n<=0) Return
If (incx==1 .And. incy==1) Goto 20
!
ix = 1
iy = 1
If (incx<0) ix = (-n+1)*incx + 1
If (incy<0) iy = (-n+1)*incy + 1
Do i = 1, n
dtemp = dtemp + dx(ix)*dy(iy)
ix = ix + incx
iy = iy + incy
End Do
ddot = dtemp
Return
!
! CODE FOR BOTH INCREMENTS EQUAL TO 1
!
20 m = mod(n, 5)
If (m==0) Goto 40
Do i = 1, m
dtemp = dtemp + dx(i)*dy(i)
End Do
If (n<5) Goto 60
40 mp1 = m + 1
Do i = mp1, n, 5
dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) + dx(i+2)*dy(i+2) &

• dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4)
  End Do
  60 ddot = dtemp
  Return
  End Function ddot
  Subroutine dgbfa(abd, lda, n, ml, mu, ipvt, info)

!*********************************************************************72
!
!c DGBFA factors a double precision band matrix by elimination.
!
Integer lda, n, ml, mu, ipvt(1), info
Double Precision abd(lda, 1)

Double Precision t
Integer i, idamax, i0, j, ju, jz, j0, j1, k, kp1, l, lm, m, mm, nm1
!
m = ml + mu + 1
info = 0
!
! zero initial fill-in columns
!
j0 = mu + 2
j1 = min0(n, m) - 1
If (j1<j0) Goto 30
Do jz = j0, j1
i0 = m + 1 - jz
Do i = i0, ml
abd(i, jz) = 0.0D0
End Do
End Do
30 Continue
jz = j1
ju = 0
!
! gaussian elimination with partial pivoting
!
nm1 = n - 1
If (nm1<1) Goto 130
Do k = 1, nm1
kp1 = k + 1
!
! zero next fill-in column
!
jz = jz + 1
If (jz>n) Goto 50
If (ml<1) Goto 50
Do i = 1, ml
abd(i, jz) = 0.0D0
End Do
50 Continue
!
! find l = pivot index
!
lm = min0(ml, n-k)
l = idamax(lm+1, abd(m,k), 1) + m - 1
ipvt(k) = l + k - m
!
! zero pivot implies this column already triangularized
!
If (abd(l,k)==0.0D0) Goto 100
!
! interchange if necessary
!
If (l==m) Goto 60
t = abd(l, k)
abd(l, k) = abd(m, k)
abd(m, k) = t
60 Continue
!
! compute multipliers
!
t = -1.0D0/abd(m, k)
Call dscal(lm, t, abd(m+1,k), 1)
!
! row elimination with column indexing
!
ju = min0(max0(ju,mu+ipvt(k)), n)
mm = m
If (ju<kp1) Goto 90
Do j = kp1, ju
l = l - 1
mm = mm - 1
t = abd(l, j)
If (l==mm) Goto 70
abd(l, j) = abd(mm, j)
abd(mm, j) = t
70 Continue
Call daxpy(lm, t, abd(m+1,k), 1, abd(mm+1,j), 1)
End Do
90 Continue
Goto 110
100 Continue
info = k
110 Continue
End Do
130 Continue
ipvt(n) = n
If (abd(m,n)==0.0D0) info = n
Return
End Subroutine dgbfa
Subroutine dgbsl(abd, lda, n, ml, mu, ipvt, b, job)

!*********************************************************************72
!
!c DGBSL solves a double precision band linear system.
!
Integer lda, n, ml, mu, ipvt(1), job
Double Precision abd(lda, 1), b(1)

Double Precision ddot, t
Integer k, kb, l, la, lb, lm, m, nm1

m = mu + ml + 1
nm1 = n - 1
If (job/=0) Goto 50
!
! job = 0 , solve a * x = b
! first solve l*y = b
!
If (ml==0) Goto 30
If (nm1<1) Goto 30
Do k = 1, nm1
lm = min0(ml, n-k)
l = ipvt(k)
t = b(l)
If (l==k) Goto 10
b(l) = b(k)
b(k) = t
10 Continue
Call daxpy(lm, t, abd(m+1,k), 1, b(k+1), 1)
End Do
30 Continue
!
! now solve u*x = y
!
Do kb = 1, n
k = n + 1 - kb
b(k) = b(k)/abd(m, k)
lm = min0(k, m) - 1
la = m - lm
lb = k - lm
t = -b(k)
Call daxpy(lm, t, abd(la,k), 1, b(lb), 1)
End Do
Goto 100
50 Continue
!
! job = nonzero, solve trans(a) * x = b
! first solve trans(u)*y = b
!
Do k = 1, n
lm = min0(k, m) - 1
la = m - lm
lb = k - lm
t = ddot(lm, abd(la,k), 1, b(lb), 1)
b(k) = (b(k)-t)/abd(m, k)
End Do
!
! now solve trans(l)*x = y
!
If (ml==0) Goto 90
If (nm1<1) Goto 90
Do kb = 1, nm1
k = n - kb
lm = min0(ml, n-k)
b(k) = b(k) + ddot(lm, abd(m+1,k), 1, b(k+1), 1)
l = ipvt(k)
If (l==k) Goto 70
t = b(l)
b(l) = b(k)
b(k) = t
70 Continue
End Do
90 Continue
100 Continue
Return
End Subroutine dgbsl
Subroutine dgefa(a, lda, n, ipvt, info)

!*********************************************************************72
!
!c DGEFA factors a double precision matrix by gaussian elimination.
!
Integer lda, n, ipvt(1), info
Double Precision a(lda, 1)

Double Precision t
Integer idamax, j, k, kp1, l, nm1
!
! gaussian elimination with partial pivoting
!
info = 0
nm1 = n - 1
If (nm1<1) Goto 70
Do k = 1, nm1
kp1 = k + 1
!
! find l = pivot index
!
l = idamax(n-k+1, a(k,k), 1) + k - 1
ipvt(k) = l
!
! zero pivot implies this column already triangularized
!
If (a(l,k)==0.0D0) Goto 40
!
! interchange if necessary
!
If (l==k) Goto 10
t = a(l, k)
a(l, k) = a(k, k)
a(k, k) = t
10 Continue
!
! compute multipliers
!
t = -1.0D0/a(k, k)
Call dscal(n-k, t, a(k+1,k), 1)
!
! row elimination with column indexing
!
Do j = kp1, n
t = a(l, j)
If (l==k) Goto 20
a(l, j) = a(k, j)
a(k, j) = t
20 Continue
Call daxpy(n-k, t, a(k+1,k), 1, a(k+1,j), 1)
End Do
Goto 50
40 Continue
info = k
50 Continue
End Do
70 Continue
ipvt(n) = n
If (a(n,n)==0.0D0) info = n
Return
End Subroutine dgefa
Subroutine dgesl(a, lda, n, ipvt, b, job)

!*********************************************************************72
!
!c DGESL solves the double precision system
! a * x = b or trans(a) * x = b
! using the factors computed by dgeco or dgefa.
!
Integer lda, n, ipvt(1), job
Double Precision a(lda, 1), b(1)

Double Precision ddot, t
Integer k, kb, l, nm1

nm1 = n - 1
If (job/=0) Goto 50
!
! job = 0 , solve a * x = b
! first solve ly = b
!
If (nm1<1) Goto 30
Do k = 1, nm1
l = ipvt(k)
t = b(l)
If (l==k) Goto 10
b(l) = b(k)
b(k) = t
10 Continue
Call daxpy(n-k, t, a(k+1,k), 1, b(k+1), 1)
End Do
30 Continue
!
! now solve ux = y
!
Do kb = 1, n
k = n + 1 - kb
b(k) = b(k)/a(k, k)
t = -b(k)
Call daxpy(k-1, t, a(1,k), 1, b(1), 1)
End Do
Goto 100
50 Continue
!
! job = nonzero, solve trans(a) * x = b
! first solve trans(u)*y = b
!
Do k = 1, n
t = ddot(k-1, a(1,k), 1, b(1), 1)
b(k) = (b(k)-t)/a(k, k)
End Do
!
! now solve trans(l)*x = y
!
If (nm1<1) Goto 90
Do kb = 1, nm1
k = n - kb
b(k) = b(k) + ddot(n-k, a(k+1,k), 1, b(k+1), 1)
l = ipvt(k)
If (l==k) Goto 70
t = b(l)
b(l) = b(k)
b(k) = t
70 Continue
End Do
90 Continue
100 Continue
Return
End Subroutine dgesl
Subroutine digit_inc(c)

!*********************************************************************72
!
!c DIGIT_INC increments a decimal digit.
!
Implicit None

Character c
Integer digit

Call ch_to_digit(c, digit)

If (digit==-1) Then
Return
End If

digit = digit + 1

If (digit==10) Then
digit = 0
End If

Call digit_to_ch(digit, c)

Return
End Subroutine digit_inc
Subroutine digit_to_ch(digit, c)

!*********************************************************************72
!
!c DIGIT_TO_CH returns the character representation of a decimal digit.
!
Implicit None

Character c
Integer digit

If (0<=digit .And. digit<=9) Then

c = char(digit+48)

Else

c = '*'

End If

Return
End Subroutine digit_to_ch
Subroutine dscal(n, da, dx, incx)

!*********************************************************************72
!
!c DSCAL scales a vector by a constant.
!
Double Precision da, dx(1)
Integer i, incx, m, mp1, n, nincx

If (n<=0) Return
If (incx==1) Goto 20
!
! CODE FOR INCREMENT NOT EQUAL TO 1
!
nincx = nincx
Do i = 1, nincx, incx
dx(i) = dadx(i)
End Do
Return
!
! CODE FOR INCREMENT EQUAL TO 1
!
20 m = mod(n, 5)
If (m==0) Goto 40
Do i = 1, m
dx(i) = da*dx(i)
End Do
If (n<5) Return
40 mp1 = m + 1
Do i = mp1, n, 5
dx(i) = da*dx(i)
dx(i+1) = da*dx(i+1)
dx(i+2) = da*dx(i+2)
dx(i+3) = da*dx(i+3)
dx(i+4) = da*dx(i+4)
End Do
Return
End Subroutine dscal
Subroutine element_node_bandwidth(maxnd, maxel, nnodes, nelemn,&
node, neqn, np, indx, nlband, nuband, nband)

!*********************************************************************72
!
!c ELEMENT_NODE_BANDWIDTH determines the bandwidth associated with the grid.
!
Implicit None

Integer maxnd
Integer maxel
Integer nelemn
Integer neqn
Integer nnodes
Integer np

Integer dof_max(neqn)
Integer dof_min(neqn)
Integer element
Integer i
Integer ieqn
Integer indx(maxnd, 3)
Integer l1
Integer l2
Integer n1
Integer n2
Integer nband
Integer nlband
Integer node(maxel, nnodes)
Integer nuband
Integer p1
Integer p2
Integer u1
Integer u2
Integer v1
Integer v2

Do ieqn = 1, neqn
dof_min(ieqn) = ieqn
dof_max(ieqn) = ieqn
End Do

Do element = 1, nelemn

Do l1 = 1, nnodes

  n1 = node(element, l1)

  u1 = indx(n1, 1)
  v1 = indx(n1, 2)
  p1 = indx(n1, 3)

  Do l2 = 1, nnodes

    n2 = node(element, l2)

    u2 = indx(n2, 1)
    v2 = indx(n2, 2)
    p2 = indx(n2, 3)

    If (1<=u1 .And. 1<=u2) Then
      dof_min(u1) = min(dof_min(u1), u2)
      dof_max(u1) = max(dof_max(u1), u2)
    End If

    If (1<=u1 .And. 1<=v2) Then
      dof_min(u1) = min(dof_min(u1), v2)
      dof_max(u1) = max(dof_max(u1), v2)
    End If

    If (1<=u1 .And. 1<=p2) Then
      dof_min(u1) = min(dof_min(u1), p2)
      dof_max(u1) = max(dof_max(u1), p2)
    End If

    If (1<=v1 .And. 1<=u2) Then
      dof_min(v1) = min(dof_min(v1), u2)
      dof_max(v1) = max(dof_max(v1), u2)
    End If

    If (1<=v1 .And. 1<=v2) Then
      dof_min(v1) = min(dof_min(v1), v2)
      dof_max(v1) = max(dof_max(v1), v2)
    End If

    If (1<=v1 .And. 1<=p2) Then
      dof_min(v1) = min(dof_min(v1), p2)
      dof_max(v1) = max(dof_max(v1), p2)
    End If

    If (1<=p1 .And. 1<=u2) Then
      dof_min(p1) = min(dof_min(p1), u2)
      dof_max(p1) = max(dof_max(p1), u2)
    End If

    If (1<=p1 .And. 1<=v2) Then
      dof_min(p1) = min(dof_min(p1), v2)
      dof_max(p1) = max(dof_max(p1), v2)
    End If

  End Do
End Do

End Do

nlband = 0
nuband = 0
Do ieqn = 1, neqn
nlband = max(nlband, ieqn-dof_min(ieqn))
nuband = max(nuband, dof_max(ieqn)-ieqn)
End Do

nband = nlband + nuband + 1

Write (, '(a)') ' '
Write (, '(a)') 'ELEMENT_NODE_BANDWIDTH:'
Write (, '(a,i6)') ' Lower half bandwidth = ', nlband
Write (, '(a,i6)') ' Upper half bandwidth = ', nuband
Write (*, '(a,i6)') ' Total bandwidth = ', nband

Return
End Subroutine element_node_bandwidth
Subroutine element_node_write(file_name, element_max, element_num, element_node)

!*********************************************************************72
!
!c ELEMENT_NODE_WRITE writes the element node list to a file.
!
Implicit None

Integer element_max

Integer element
Integer element_node(element_max, 6)
Integer element_num
Character (Len=*) file_name
Integer j

Open (Unit=1, File=file_name, Err=10)

Do element = 1, element_num
Write (1, *)(element_node(element,j), j=1, 6)
End Do

Close (Unit=1)

Write (, '(a)') ' '
Write (, '(a)') 'ELEMENT_NODE_WRITE:'
Write (, '(a)') ' The 3-node triangle nodal coordinates were'
Write (, '(a)') ' written to the output file:'
Write (, '(a)') ' '
Write (, '(4x,a)') file_name
Return

10 Continue

Write (, '(a)') ' '
Write (, '(a)') 'ELEMENT_NODE_WRITE - Fatal errorc'
Write (, '(a)') ' Could not open the output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') file_name

Stop
End Subroutine element_node_write
Subroutine file_name_inc(file_name)

!*********************************************************************72
!
!c FILE_NAME_INC generates the next file name in a series.
!
Implicit None

Character c
Logical ch_is_digit
Character () file_name
Integer i
Integer lens

lens = len(file_name)

Do i = lens, 1, -1

c = file_name(i:i)

If (ch_is_digit(c)) Then

  Call digit_inc(c)

  file_name(i:i) = c

  If (c/='0') Then
    Return
  End If

End If

End Do

Return
End Subroutine file_name_inc
Subroutine hcell_dof_count(region_density_x, region_density_y, dof_num)

!*********************************************************************72
!
!c HCELL_DOF_COUNT determines the number of degrees of freedom in the region.
!
Implicit None

Integer dof_num
Integer region_density_x(3)
Integer region_density_y(3)
Integer node_num_p
Integer node_num_u

node_num_u = (2region_density_x(1)+1)(2region_density_y(1)+1) + &
(2region_density_x(1)+1)(2region_density_y(3)+1) + &
(2region_density_x(2)-1)(2region_density_y(1)+1) + &
(2region_density_x(2)+1)(2region_density_y(2)-1) + &
(2region_density_x(2)-1)(2region_density_y(3)+1) + &
(2region_density_x(3)+1)(2region_density_y(1)+1) + &
(2region_density_x(3)+1)(2*region_density_y(3)+1)

node_num_p = (region_density_x(1)+1)(region_density_y(1)+1) + &
(region_density_x(1)+1)(region_density_y(3)+1) +&
(region_density_x(2)-1)(region_density_y(1)+1) +&
(region_density_x(2)+1)(region_density_y(2)-1) +&
(region_density_x(2)-1)(region_density_y(3)+1) +&
(region_density_x(3)+1)(region_density_y(1)+1) +&
(region_density_x(3)+1)*(region_density_y(3)+1)

dof_num = 2*node_num_u + node_num_p

Write (, '(a)') ' '
Write (, '(a)') 'HCELL_DOF_COUNT:'
Write (*, '(a,i6)') ' Number of degrees of freedom = ', dof_num

Return
End Subroutine hcell_dof_count