!*************************************************************
!*************************************************************
! 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之b
!*************************************************************
! SINOMACH
!
! 纳维斯托克斯2D-H形领域有限元法估算解决（N-S）方程(abc部)之b
!
! 周勇 20221125
!**************************************************************
!c $$$ Run, Output:
!**************************************************************
!**************************************************************
Subroutine hcell_dof_set(region_density_x, region_density_y, maxnd,&
node_num, node_dof_index)
!*********************************************************************72
!
!c HCELL_DOF_SET assigns degrees of freedom to each node.
!
Implicit None
Integer maxnd
Integer col
Integer dof
Logical found
Integer i
Integer j
Integer node
Integer node_dof_index(maxnd, 3)
Integer node_num
Integer region_density_x(3)
Integer region_density_y(3)
Integer row
node = 0
dof = 0
j = 0
!
! Working in column 1, and NOT handling nodes in extreme right.
!
Do col = 1, 2*region_density_x(1)
!
! Working in row 1.
!
j = j + 1
i = 0
Do row = 1, 2*region_density_y(1) + 1
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
!
! Working in row 3.
!
i = i + 2*region_density_y(2) - 1
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
End Do
!
! Working in column 2, including nodes on extreme right and left.
!
Do col = 1, 2*region_density_x(2) + 1
!
! Working in row 1.
!
j = j + 1
i = 0
Do row = 1, 2*region_density_y(1) + 1
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
!
! Working in row 2.
!
Do row = 2, 2*region_density_y(2)
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
!
! Working in row 3.
!
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
End Do
!
! Working in column 3, and NOT handling nodes in extreme left.
!
Do col = 2, 2*region_density_x(3) + 1
!
! Working in row 1.
!
j = j + 1
i = 0
Do row = 1, 2*region_density_y(1) + 1
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
!
! Working in row 3.
!
i = i + 2*region_density_y(2) - 1
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1
dof = dof + 1
node_dof_index(node, 1) = dof
dof = dof + 1
node_dof_index(node, 2) = dof
If (mod(j,2)==1 .And. mod(i,2)==1) Then
dof = dof + 1
node_dof_index(node, 3) = dof
Else
node_dof_index(node, 3) = 0
End If
End Do
End Do
!
! Replace one continuity equation by a pressure = 0 equation.
!
found = .False.
Do node = 1, node_num
If (0<node_dof_index(node,3)) Then
node_dof_index(node, 3) = -1
found = .True.
Exit
End If
End Do
If (.Not. found) Then
Write (, '(a)') ' '
Write (, '(a)') 'HCELL_DOF_SET - Fatal errorc'
Write (, '(a)') ' Could not find a degree of freedom'
Write (, '(a)') ' associated with the continuity equation.'
Stop
End If
Return
End Subroutine hcell_dof_set
Subroutine hcell_element_count(region_density_x, region_density_y, element_num)
!*********************************************************************72
!
!c HCELL_ELEMENT_COUNT determines the number of elements in the region.
!
Implicit None
Integer element_num
Integer region_density_x(3)
Integer region_density_y(3)
element_num = 2*region_density_x(1)region_density_y(1) + &
2region_density_x(1)region_density_y(3) + &
2region_density_x(2)region_density_y(1) +&
2region_density_x(2)region_density_y(2) + &
2region_density_x(2)region_density_y(3) +&
2region_density_x(3)region_density_y(1) + &
2region_density_x(3)*region_density_y(3)
Write (, '(a)') ' '
Write (, '(a)') 'HCELL_ELEMENT_COUNT:'
Write (*, '(a,i6)') ' Number of elements = ', element_num
Return
End Subroutine hcell_element_count
Subroutine hcell_element_node(region_density_x, region_density_y,&
maxel, element_num, nnodes, element_node)
!*********************************************************************72
!
!c HCELL_ELEMENT_NODE determines the nodes that make up each element.
!
Implicit None
Integer maxel
Integer nnodes
Integer col
Integer col2
Integer element
Integer element_node(maxel, nnodes)
Integer element_num
Integer inc1
Integer inc2
Integer node_sw
Integer region_density_x(3)
Integer region_density_y(3)
Integer row
Integer row2
element = 0
Do col = 1, 3
Do col2 = 1, region_density_x(col)
Do row = 1, 3
If (row/=2 .Or. col==2) Then
If (col==1) Then
If (col2<region_density_x(1)) Then
If (row==1) Then
If (col2==1) Then
node_sw = 1
Else
node_sw = node_sw + inc1 + 1
End If
Else
node_sw = node_sw + 1
End If
inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
Else If (row==1) Then
node_sw = node_sw + inc1 + 1
inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
Else If (row==3) Then
node_sw = node_sw + 1
inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(2)-1)&
+ (2*region_density_y(3)+1)
End If
Else If (col==2) Then
If (row==1) Then
node_sw = node_sw + inc1 + 1
End If
inc1 = (2*region_density_y(1)+1) + (2*region_density_y(2)-1)&
+ (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(2)-1)&
+ (2*region_density_y(3)+1)
Else If (col==3) Then
If (1<col2) Then
If (row==1) Then
node_sw = node_sw + inc1 + 1
Else
node_sw = node_sw + 1
End If
inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
Else If (row==1) Then
node_sw = node_sw + inc1 + 1
inc1 = (2*region_density_y(1)+1) +&
(2*region_density_y(2)-1) + (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
Else If (row==3) Then
node_sw = node_sw + (2*region_density_y(2)-1) + 1
inc1 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
inc2 = (2*region_density_y(1)+1) + (2*region_density_y(3)+1)
End If
End If
Do row2 = 1, region_density_y(row)
element = element + 1
element_node(element, 1) = node_sw
element_node(element, 2) = node_sw + inc1 + inc2
element_node(element, 3) = node_sw + 2
element_node(element, 4) = node_sw + inc1
element_node(element, 5) = node_sw + inc1 + 1
element_node(element, 6) = node_sw + 1
element = element + 1
element_node(element, 1) = node_sw + inc1 + inc2 + 2
element_node(element, 2) = node_sw + 2
element_node(element, 3) = node_sw + inc1 + inc2
element_node(element, 4) = node_sw + inc1 + 2
element_node(element, 5) = node_sw + inc1 + 1
element_node(element, 6) = node_sw + inc1 + inc2 + 1
node_sw = node_sw + 2
End Do
End If
End Do
End Do
End Do
Return
End Subroutine hcell_element_node
Subroutine hcell_node_count(region_density_x, region_density_y, node_num)
!*********************************************************************72
!
!c HCELL_NODE_NUM determines the number of nodes in the region.
!
Implicit None
Integer node_num
Integer region_density_x(3)
Integer region_density_y(3)
!
! Count the nodes.
!
node_num = (2region_density_x(1)+1)(2*region_density_y(1)+1)&
Write (, '(a)') ' '
Write (, '(a)') 'HCELL_NODE_COUNT:'
Write (*, '(a,i6)') ' Number of nodes = ', node_num
Return
End Subroutine hcell_node_count
Subroutine hcell_node_xy(region_density_x, region_density_y,&
node_num, region_x, region_y, node_x, node_y)
!*********************************************************************72
!
!c HCELL_NODE_XY assigns coordinates to each node.
!
Implicit None
Integer node_num
Integer col
Integer i
Integer j
Integer node
Double Precision node_x(node_num)
Double Precision node_y(node_num)
Integer region_density_x(3)
Integer region_density_y(3)
Double Precision region_x(4)
Double Precision region_y(4)
Integer row
node = 0
j = 0
!
! Working in column 1, except for the extreme right.
!
Do col = 1, 2*region_density_x(1)
j = j + 1
i = 0
!
! Working in row 1.
!
Do row = 1, 2*region_density_y(1) + 1
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(1)+1-col)*region_x(1)+&
dble(-1+col)*region_x(2))/dble(2*region_density_x(1))
node_y(node) = (dble(2*region_density_y(1)+1-row)*region_y(1)+&
dble(-1+row)*region_y(2))/dble(2*region_density_y(1))
End Do
!
! Working in row 3.
!
i = i + 2*region_density_y(2) - 1
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(1)+1-col)*region_x(1)+&
dble(-1+col)*region_x(2))/dble(2*region_density_x(1))
node_y(node) = (dble(2*region_density_y(3)+1-row)*region_y(3)+&
dble(-1+row)*region_y(4))/dble(2*region_density_y(3))
End Do
End Do
!
! Working in column 2, including extreme left and right.
!
Do col = 1, 2*region_density_x(2) + 1
!
! Working in row 1.
!
j = j + 1
i = 0
Do row = 1, 2*region_density_y(1) + 1
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(2)+1-col)*region_x(2)+&
dble(-1+col)*region_x(3))/dble(2*region_density_x(2))
node_y(node) = (dble(2*region_density_y(1)+1-row)*region_y(1)+&
dble(-1+row)*region_y(2))/dble(2*region_density_y(1))
End Do
!
! Working in row 2.
!
Do row = 2, 2*region_density_y(2)
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(2)+1-col)*region_x(2)+&
dble(-1+col)*region_x(3))/dble(2*region_density_x(2))
node_y(node) = (dble(2*region_density_y(2)+1-row)*region_y(2)+&
dble(-1+row)*region_y(3))/dble(2*region_density_y(2))
End Do
!
! Working in row 3.
!
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(2)+1-col)*region_x(2)+&
dble(-1+col)*region_x(3))/dble(2*region_density_x(2))
node_y(node) = (dble(2*region_density_y(3)+1-row)*region_y(3)+&
dble(-1+row)*region_y(4))/dble(2*region_density_y(3))
End Do
End Do
!
! Working in column 3, except for extreme left.
!
Do col = 2, 2*region_density_x(3) + 1
!
! Working in row 1.
!
j = j + 1
i = 0
Do row = 1, 2*region_density_y(1) + 1
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(3)+1-col)*region_x(3)+&
dble(-1+col)*region_x(4))/dble(2*region_density_x(3))
node_y(node) = (dble(2*region_density_y(1)+1-row)*region_y(1)+&
dble(-1+row)*region_y(2))/dble(2*region_density_y(1))
End Do
i = i + 2*region_density_y(2) - 1
!
! Working in row 3.
!
Do row = 1, 2*region_density_y(3) + 1
i = i + 1
node = node + 1
node_x(node) = (dble(2*region_density_x(3)+1-col)&
*region_x(3)+dble(-1+col)*region_x(4))/dble(2*region_density_x(3))
node_y(node) = (dble(2*region_density_y(3)+1-row)&
*region_y(3)+dble(-1+row)*region_y(4))/dble(2*region_density_y(3))
End Do
End Do
Return
End Subroutine hcell_node_xy
Integer Function idamax(n, dx, incx)
!*********************************************************************72
!
!c IDAMAX finds the vector element of largest magnitude.
!
Double Precision dx(1), dmax
Integer i, incx, ix, n
idamax = 0
If (n<1) Return
idamax = 1
If (n==1) Return
If (incx==1) Goto 20
!
! CODE FOR INCREMENT NOT EQUAL TO 1
!
ix = 1
dmax = dabs(dx(1))
ix = ix + incx
Do i = 2, n
If (dabs(dx(ix))<=dmax) Goto 5
idamax = i
dmax = dabs(dx(ix))
5 ix = ix + incx
End Do
Return
!
! CODE FOR INCREMENT EQUAL TO 1
!
20 dmax = dabs(dx(1))
Do i = 2, n
If (dabs(dx(i))<=dmax) Goto 30
idamax = i
dmax = dabs(dx(i))
30 End Do
Return
End Function idamax
Subroutine node_eps(node_file_name, np, node_mask, xc, yc, title)
!*********************************************************************72
!
!c NODE_EPS creates an EPS file containing an image of the nodes.
!
Implicit None
Integer np
Double Precision ave_x
Double Precision ave_y
Integer circle_size
Double Precision dif
Integer eps_unit
Integer eps_x
Integer eps_y
Character (Len=*) node_file_name
Integer i
Integer ios
Integer ip1
Integer j
Integer k
Integer node
Logical node_mask(np)
Integer node_unmasked_num
Double Precision node_x_max
Double Precision node_x_min
Double Precision node_y_max
Double Precision node_y_min
Double Precision r8_huge
Double Precision scale
Character (Len=40) string
Character (Len=*) title
Double Precision xc(np)
Double Precision yc(np)
!
! Determine the range and number of the unmasked nodes.
!
node_x_min = r8_huge()
node_x_max = -r8_huge()
node_y_min = r8_huge()
node_y_max = -r8_huge()
node_unmasked_num = 0
Do node = 1, np
If (node_mask(node)) Then
node_unmasked_num = node_unmasked_num + 1
node_x_min = min(node_x_min, xc(node))
node_x_max = max(node_x_max, xc(node))
node_y_min = min(node_y_min, yc(node))
node_y_max = max(node_y_max, yc(node))
End If
End Do
If (node_unmasked_num<=200) Then
circle_size = 3
Else If (node_unmasked_num<=800) Then
circle_size = 2
Else
circle_size = 1
End If
If (node_y_max-node_y_min<node_x_max-node_x_min) Then
scale = node_x_max - node_x_min
dif = (node_x_max-node_x_min) - (node_y_max-node_y_min)
node_y_max = node_y_max + 0.5dif
node_y_min = node_y_min - 0.5dif
Else
scale = node_y_max - node_y_min
dif = (node_y_max-node_y_min) - (node_x_max-node_x_min)
node_x_max = node_x_max + 0.5dif
node_x_min = node_x_min - 0.5dif
End If
eps_unit = 1
Open (Unit=eps_unit, File=node_file_name, Status='unknown')
Write (eps_unit, '(a)') '%cPS-Adobe-3.0 EPSF-3.0'
Write (eps_unit, '(a)') '%%Creator: node_eps(hcell.f)'
Write (eps_unit, '(a,a)') '%% : ', node_file_name
Write (eps_unit, '(a)') '%%Pages: 1'
Write (eps_unit, '(a)') '%%BoundingBox: 36 36 576 756'
Write (eps_unit, '(a)') '%%Document-Fonts: Times-Roman'
Write (eps_unit, '(a)') '%%LanguageLevel: 1'
Write (eps_unit, '(a)') '%%EndComments'
Write (eps_unit, '(a)') '%%BeginProlog'
Write (eps_unit, '(a)') '/inch {72 mul} def'
Write (eps_unit, '(a)') '%%EndProlog'
Write (eps_unit, '(a)') '%%Page: 1 1'
Write (eps_unit, '(a)') 'save'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Set RGB line color.'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.9000 0.9000 0.9000 setrgbcolor'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Draw a gray border around the page.'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') 'newpath'
Write (eps_unit, '(a)') ' 36 126 moveto'
Write (eps_unit, '(a)') ' 576 126 lineto'
Write (eps_unit, '(a)') ' 576 666 lineto'
Write (eps_unit, '(a)') ' 36 666 lineto'
Write (eps_unit, '(a)') ' 36 126 lineto'
Write (eps_unit, '(a)') 'stroke'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Set RGB line color.'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Label the plot:'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor'
Write (eps_unit, '(a)') '/Times-Roman findfont 0.50 inch scalefont setfont'
Write (eps_unit, '(a)') ' 36 666 moveto'
Write (eps_unit, '(a,a,a)') '(',title, ') show'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Define a clipping polygon'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 36 126 moveto'
Write (eps_unit, '(a)') ' 576 126 lineto'
Write (eps_unit, '(a)') ' 576 666 lineto'
Write (eps_unit, '(a)') ' 36 666 lineto'
Write (eps_unit, '(a)') ' 36 126 lineto'
Write (eps_unit, '(a)') 'clip newpath'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Draw filled dots at each node:'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 1.0000 setrgbcolor'
Do node = 1, np
If (node_mask(node)) Then
eps_x = int((node_x_max-xc(node))*61.0+(+xc(node)-node_x_min)*551.0)/scale
eps_y = int((node_y_max-yc(node))*151.0+(yc(node)-node_y_min)*641.0)/scale
Write (eps_unit, '(a,i4,2x,i4,2x,i4,a)') 'newpath ',&
eps_x, eps_y, circle_size, ' 0 360 arc closepath fill'
End If
End Do
!
! Label each node, but only if there aren't too many of themc
!
If (node_unmasked_num<=200) Then
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% Label the nodes:'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') ' 0.0000 0.0000 0.0000 setrgbcolor'
Write (eps_unit, '(a)') '/Times-Roman findfont 0.20 inch scalefont setfont'
Do node = 1, np
If (node_mask(node)) Then
eps_x = int((node_x_max-xc(node))*61.0+(+xc(node)-node_x_min)*551.0)/scale
eps_y = int((node_y_max-yc(node))*151.0+(yc(node)-node_y_min)*641.0)/scale
Write (eps_unit, '(i4,2x,i4,a,i4,a)') eps_x, eps_y + 5,&
' moveto (', node, ') show'
End If
End Do
End If
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') 'restore showpage'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '% End of page'
Write (eps_unit, '(a)') '%'
Write (eps_unit, '(a)') '%%Trailer'
Write (eps_unit, '(a)') '%%EOF'
Close (Unit=eps_unit)
Write (, '(a)') ' '
Write (, '(a)') 'NODE_EPS:'
Write (, '(a)') ' An encapsulated Post file was created'
Write (, '(a)') ' containing an image of the nodes.'
Write (, '(a)') ' The file name is'
Write (, '(4x,a)') node_file_name
Return
End Subroutine node_eps
Subroutine 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)
!*********************************************************************72
!
!c NSTOKE solves the Navier-Stokes equations using Taylor-Hood elements.
!
Implicit Double Precision (A-H, O-Z)
Double Precision a(mrow1, )
Double Precision area()
Double Precision f()
Double Precision g()
Integer indx(maxnd, 3)
Integer ipivot()
Integer it
Integer iter
Integer node(maxel, )
Double Precision reynld
Double Precision un(2)
Double Precision unx(2)
Double Precision uny(2)
Double Precision uold()
Double Precision uold_qp
Double Precision visc
Double Precision vold_qp
Double Precision x
Double Precision xc()
Double Precision xm(maxel, )
Double Precision y
Double Precision yc()
Double Precision ym(maxel, *)
!
! zero arrays
! g array contains the previous iterate
! f array contains the right hand side initially and then the current
! iterate is overwritten on f
!
visc = 1.D0/reynld
!
! matrix assembly triangle by triangle
! nsim is the number of simple iterations performed
!
csim = 0.D0
Do iter = 1, nsteps
niter = iter
If (iter>nsim) csim = 1.D0
Do i = 1, nrow1
Do j = 1, ncol1
a(i, j) = 0.D0
End Do
End Do
Do it = 1, nelemn
arr = area(it)/3.D0
Do iquad = 1, nquad
y = ym(it, iquad)
x = xm(it, iquad)
Call trans(it, x, y, det, xix, xiy, etax, etay, xc, yc, node, maxel)
ar = arr*det
Do kk = 1, 2
un(kk) = 0.D0
uny(kk) = 0.D0
unx(kk) = 0.D0
End Do
uold_qp = 0.0D0
vold_qp = 0.0D0
Do iq = 1, nnodes
Call refqbf(x, y, iq, bb, tbx, tby)
bx = tbx*xix + tby*etax
by = tbx*xiy + tby*etay
ip = node(it, iq)
Do iuk = 1, 2
iun = indx(ip, iuk)
If (0<iun) Then
un(iuk) = un(iuk) + bb*g(iun)
unx(iuk) = unx(iuk) + bx*g(iun)
uny(iuk) = uny(iuk) + by*g(iun)
If (iuk==1) uold_qp = uold_qp + bb*uold(iun)
If (iuk==2) vold_qp = vold_qp + bb*uold(iun)
Else If (iun<0) Then
ubc = alpha*ubdry(iuk, ip, xc, yc)
un(iuk) = un(iuk) + bb*ubc
unx(iuk) = unx(iuk) + bx*ubc
uny(iuk) = uny(iuk) + by*ubc
If (iuk==1) uold_qp = uold_qp + bb*ubc
If (iuk==2) vold_qp = vold_qp + bb*ubc
End If
End Do
End Do
Do iq = 1, nnodes
ip = node(it, iq)
Call refqbf(x, y, iq, bb, tbx, tby)
bx = tbx*xix + tby*etax
by = tbx*xiy + tby*etay
bbl = refbsp(x, y, iq)
Do iuk = 1, nuk
i = indx(ip, iuk)
If (i<=0) Goto 210
If (iuk==1) f(i) = f(i) + csim*((un(1)*unx(1)+un(2)*uny(1))&
*bb)*ar + rdel*uold_qp*bb*ar
If (iuk==2) f(i) = f(i) + csim*((un(1)*unx(2)+un(2)*uny(2))&
*bb)*ar + rdel*vold_qp*bb*ar
Do iqq = 1, nnodes
ipp = node(it, iqq)
Call refqbf(x, y, iqq, bbb, tbbx, tbby)
bbx = tbbx*xix + tbby*etax
bby = tbbx*xiy + tbby*etay
bbbl = refbsp(x, y, iqq)
Do iukk = 1, nuk
j = indx(ipp, iukk)
If (j==0) Goto 190
aij = 0.D0
If (i==neqn) Goto 190
If (iuk==1) Then
If (iukk==1) aij = visc*(by*bby+bx*bbx)&
+ (bbb*unx(1)*bb)*csim + bb*bbx*un(1) +&
bb*bby*un(2) + rdel*(bb*bbb)
If (iukk==2) aij = csim*(bb*bbb*uny(1))
If (iukk==3) aij = -bx*bbbl
Else If (iuk==2) Then
If (iukk==2) aij = (visc*(by*bby+bx*bbx)+&
(bb*bbb*uny(2))*csim+bb*bby*un(2)+bb*bbx*un(1)) + rdel*(bb*bbb)
If (iukk==1) aij = csim*(bb*bbb*unx(2))
If (iukk==3) aij = -by*bbbl
Else
If (iukk==1) aij = bbx*bbl
If (iukk==2) aij = bby*bbl
End If
If (j<0) Goto 185
iuse = i - j + nband
a(iuse, j) = a(iuse, j) + aij*ar
Goto 190
!
! add terms to rhs for inhomogeneous boundary condition
!
185 Continue
f(i) = f(i) - aralphaubdry(iukk, ipp, xc, yc)*aij
190 End Do
End Do
210 End Do
End Do
End Do
End Do
f(neqn) = 0.D0
Do j = neqn - nlband, neqn - 1
i = neqn - j + nband
a(i, j) = 0.D0
End Do
a(nband, neqn) = 1.D0
!
! solve system
!
job = 0
Call dgbfa(a, mrow1, neqn, nlband, nuband, ipivot, info)
Call dgbsl(a, mrow1, neqn, nlband, nuband, ipivot, f, job)
!
! check for convergence
!
diff = 0.D0
Do i = 1, neqn
diff = diff + (g(i)-f(i))**2
End Do
diff = sqrt(diff)
Write (*, *) ' Iteration ', iter, ' Difference is ', diff
If (diff<=tolns) Then
Return
End If
Do i = 1, neqn
g(i) = f(i)
f(i) = 0.D0
End Do
End Do
Return
End Subroutine nstoke
更多回帖
