!*************************************************************
!*************************************************************
! 纳维斯托克斯两维H形领域有限元法估算解决方程(abc部)之c
!*************************************************************
! SINOMACH
!
! 纳维斯托克斯2D-H形领域有限元法估算解决（N-S）方程(abc部)之c
!
! 周勇 20221125
!**************************************************************
!c $$$ Run, Output:
!**************************************************************
!**************************************************************
Subroutine p_write(p_file_name, node_max, node_num, neqn, c, indx, node_x, node_y)
!*********************************************************************72
!
!c P_WRITE writes the pressures values for a given timestep to a file.
!
Implicit None
Integer node_max
Integer node_num
Integer neqn
Double Precision c(neqn)
Character (Len=*) p_file_name
Integer i
Integer indx(node_max, 3)
Integer node
Double Precision node_x(node_num)
Double Precision node_y(node_num)
Double Precision p
Open (Unit=1, File=p_file_name, Err=10)
Do node = 1, node_num
i = indx(node, 3)
If (0==i) Then
Else If (0<i) Then
p = c(i)
Write (1, '(g14.6)') p
Else If (i<0) Then
p = 0.0D+00
Write (1, '(g14.6)') p
End If
End Do
Close (Unit=1)
Write (, '(a)') ' '
Write (, '(a)') 'P_WRITE:'
Write (, '(a)') ' Wrote the pressure output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') p_file_name
Return
10 Continue
Write (, '(a)') ' '
Write (, '(a)') 'P_WRITE - Fatal errorc'
Write (, '(a)') ' Could not open the output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') p_file_name
Stop
End Subroutine p_write
Subroutine quad_a_set(maxel, nelemn, nquad, area, xm, ym)
!*********************************************************************72
!
!c QUAD_A_SET sets quadrature information for the assembly routine.
!
Implicit None
Integer maxel
Integer nquad
Double Precision area(maxel)
Integer element
Integer nelemn
Double Precision xm(maxel, nquad)
Double Precision ym(maxel, nquad)
Do element = 1, nelemn
xm(element, 1) = 0.5D0
xm(element, 2) = 0.5D0
xm(element, 3) = 0.0D0
End Do
Do element = 1, nelemn
ym(element, 1) = 0.0D0
ym(element, 2) = 0.5D0
ym(element, 3) = 0.5D0
End Do
Do element = 1, nelemn
area(element) = 0.5D0
End Do
Return
End Subroutine quad_a_set
Function r8_huge()
!*********************************************************************72
!
!c R8_HUGE returns a "huge" double precision number.
!
Implicit None
Double Precision r8_huge
r8_huge = 1.0D+30
Return
End Function r8_huge
Function refbsp(x, y, iq)
!*********************************************************************72
!
!c REFBSP evaluates a linear basis functions on the reference triangle.
!
Implicit None
Integer iq
Double Precision refbsp
Double Precision x
Double Precision y
If (iq==1) Then
refbsp = 1.0D+00 - x - y
Else If (iq==2) Then
refbsp = x
Else If (iq==3) Then
refbsp = y
Else
Write (*, '(a)') ' '
Write (*, '(a)') 'REFBSP - Fatal errorc'
Write (*, '(a)') ' Illegal input value of IQ = ', iq
Stop
End If
Return
End Function refbsp
Subroutine refqbf(x, y, in, bb, bx, by)
!*********************************************************************72
!
!c REFQBF evaluates quadratic basis functions on the reference triangle.
!
Implicit None
Double Precision bb
Double Precision bx
Double Precision by
Integer in
Double Precision x
Double Precision y
If (in==1) Then
bb = (1.D0-x-y)(1.D0-2.D0x-2.D0y)
bx = -3.D0 + 4.D0x + 4.D0y
by = -3.D0 + 4.D0x + 4.D0*y
Else If (in==2) Then
bb = x*(2.D0*x-1.D0)
bx = 4.D0*x - 1.D0
by = 0.D0
Else If (in==3) Then
bb = y*(2.D0*y-1.D0)
bx = 0.D0
by = 4.D0*y - 1.D0
Else If (in==4) Then
bb = 4.D0*x*(1.D0-x-y)
bx = 4.D0*(1.D0-2.D0*x-y)
by = -4.D0*x
Else If (in==5) Then
bb = 4.D0*x*y
bx = 4.D0*y
by = 4.D0*x
Else If (in==6) Then
bb = 4.D0*y*(1.D0-x-y)
bx = -4.D0*y
by = 4.D0*(1.D0-x-2.D0*y)
Else
bb = 0.0D+00
bx = 0.0D+00
by = 0.0D+00
End If
Return
End Subroutine refqbf
Subroutine 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)
!*********************************************************************72
!
!c SETGRD sets up the grid for the problem.
!
Implicit None
Integer maxel
Integer maxnd
Integer nnodes
Integer nquad
Integer nuk
Double Precision area(maxel)
Integer indx(maxnd, nuk)
Integer nband
Integer nelemn
Integer neqn
Integer nlband
Integer node(maxel, nnodes)
Integer np
Integer nuband
Integer region_density_x(3)
Integer region_density_y(3)
Double Precision region_x(4)
Double Precision region_y(4)
Double Precision xc(maxnd)
Double Precision xlngth
Double Precision xm(maxel, nquad)
Double Precision yc(maxnd)
Double Precision ylngth
Double Precision ym(maxel, nquad)
!
! Determine NP, the number of nodes.
!
Call hcell_node_count(region_density_x, region_density_y, np)
!
! Assign coordinates to the nodes, XC and YC.
!
Call hcell_node_xy(region_density_x, region_density_y, np, region_x,&
region_y, xc, yc)
!
! Determine NELEMN, the number of elements.
!
Call hcell_element_count(region_density_x, region_density_y, nelemn)
!
! Assign nodes to elements in NODE.
!
Call hcell_element_node(region_density_x, region_density_y, maxel,&
nelemn, nnodes, node)
!
! Determine NEQN, the number of unknowns.
!
Call hcell_dof_count(region_density_x, region_density_y, neqn)
!
! Assign degrees of freedom in INDX.
!
Call hcell_dof_set(region_density_x, region_density_y, maxnd, np, indx)
!
! For the assembly routine, determine the quadrature data
! NQUAD, AREA, XM and YM.
!
Call quad_a_set(maxel, nelemn, nquad, area, xm, ym)
!
! Compute the bandwidths NLBAND, NUBAND and NBAND.
!
Call element_node_bandwidth(maxnd, maxel, nnodes, nelemn, node,&
neqn, np, indx, nlband, nuband, nband)
Return
End Subroutine setgrd
Subroutine timestamp()
!*********************************************************************72
!
!c TIMESTAMP prints out the current YMDHMS date as a timestamp.
!
Implicit None
Character *(8) ampm
Integer d
Character *(8) date
Integer h
Integer m
Integer mm
Character *(9) month(12)
Integer n
Integer s
Character *(10) time
Integer y
Save month
Data month/'January ', 'February ', 'March ', 'April ',&
'May ', 'June ', 'July ', 'August ', 'September',&
'October ', 'November ', 'December '/
Call date_and_time(date, time)
Read (date, '(i4,i2,i2)') y, m, d
Read (time, '(i2,i2,i2,1x,i3)') h, n, s, mm
If (h<12) Then
ampm = 'AM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Noon'
Else
ampm = 'PM'
End If
Else
h = h - 12
If (h<12) Then
ampm = 'PM'
Else If (h==12) Then
If (n==0 .And. s==0) Then
ampm = 'Midnight'
Else
ampm = 'AM'
End If
End If
End If
Write (*, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)') d,&
month(m), y, h, ':', n, ':', s, '.', mm, ampm
Return
End Subroutine timestamp
Subroutine trans(it, xq, yq, det, pj11, pj21, pj12, pj22, xc, yc, node, maxel)
!*********************************************************************72
!
!c TRANS transforms data between the reference and physical elements.
!
Implicit Double Precision (A-H, O-Z)
Double Precision xc()
Double Precision yc()
Integer node(maxel, *)
i1 = node(it, 1)
i2 = node(it, 2)
i3 = node(it, 3)
i4 = node(it, 4)
i5 = node(it, 5)
i6 = node(it, 6)
x1 = xc(i1)
y1 = yc(i1)
x2 = xc(i2)
y2 = yc(i2)
x3 = xc(i3)
y3 = yc(i3)
x4 = xc(i4)
y4 = yc(i4)
x5 = xc(i5)
y5 = yc(i5)
x6 = xc(i6)
y6 = yc(i6)
!
! Compute partial derivatives at point (xq,yq)
!
f1x = x1*(-3.D0+4.D0xq+4.D0yq) + x2*(4.D0xq-1.D0) +&
x44.D0*(1.D0-2.D0xq-yq) + x54.D0yq + x64.D0*(-yq)
f1y = x1*(-3.D0+4.D0xq+4.D0yq) + x3*(4.D0*yq-1.D0)&
Return
End Subroutine trans
Function ubdry(iuk, ip, xc, yc)
!*********************************************************************72
!
!c UBDRY evaluates the boundary conditions.
!
Implicit None
Integer ip
Integer iuk
Double Precision ubdry
Double Precision xc()
Double Precision yc()
ubdry = 0.0D+00
If (xc(ip)==0.0D+00) Then
If (iuk==1) Then
ubdry = 16.0D+00*(yc(ip)-0.5D+00)*(1.0D+00-yc(ip))
End If
End If
Return
End Function ubdry
Subroutine uv_read(uv_file_name, neqn, uold)
!*********************************************************************72
!
!c UV_READ returns an initial value for the solution coefficient vector.
!
Implicit None
Integer neqn
Character (Len=*) uv_file_name
Integer i
Double Precision uold(neqn)
Open (Unit=1, File=uv_file_name, Err=10)
Do i = 1, neqn
Read (1, *) uold(i)
End Do
Close (Unit=1)
Write (, '(a)') ' '
Write (, '(a)') 'UV_READ:'
Write (, '(a)') ' Read the file containing'
Write (, '(a)') ' the initial solution values:'
Write (, '(a)') ' '
Write (, '(4x,a)') uv_file_name
Return
10 Continue
Write (, '(a)') ' '
Write (, '(a)') 'UV_READ - Fatal errorc'
Write (, '(a)') ' Could not open the file containing'
Write (, '(a)') ' the initial solution values:'
Write (, '(a)') ' '
Write (, '(4x,a)') uv_file_name
Stop
End Subroutine uv_read
Subroutine uv_write(uv_file_name, node_max, node_num, neqn, c, indx,&
node_x, node_y)
!*********************************************************************72
!
!c UV_WRITE writes the velocity values for a given timestep to a file.
!
Implicit None
Integer node_max
Integer node_num
Integer neqn
Double Precision c(neqn)
Character (Len=*) uv_file_name
Integer i
Integer indx(node_max, 3)
Integer node
Double Precision node_x(node_num)
Double Precision node_y(node_num)
Double Precision u
Double Precision v
Open (Unit=1, File=uv_file_name, Err=10)
Do node = 1, node_num
i = indx(node, 1)
If (0<i) Then
u = c(i)
Else
u = 0.0D+00
End If
i = indx(node, 2)
If (0<i) Then
v = c(i)
Else
v = 0.0D+00
End If
Write (1, '(2g14.6)') u, v
End Do
Close (Unit=1)
Write (, '(a)') ' '
Write (, '(a)') 'UV_WRITE:'
Write (, '(a)') ' Wrote the output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') uv_file_name
Return
10 Continue
Write (, '(a)') ' '
Write (, '(a)') 'UV_WRITE - Fatal errorc'
Write (, '(a)') ' Could not open the output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') uv_file_name
Stop
End Subroutine uv_write
Subroutine xy3_write(file_name, maxnd, node_num, indx, node_x, node_y)
!*********************************************************************72
!
!c XY3_WRITE writes the coordinates of nodes of 3-node triangles to a file.
!
Implicit None
Integer maxnd
Integer node_num
Character (Len=*) file_name
Integer indx(maxnd, 3)
Integer node
Double Precision node_x(node_num)
Double Precision node_y(node_num)
Open (Unit=1, File=file_name, Err=10)
Do node = 1, node_num
If (0/=indx(node,3)) Then
Write (1, *) node_x(node), node_y(node)
End If
End Do
Close (Unit=1)
Write (, '(a)') ' '
Write (, '(a)') 'XY3_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)') 'XY3_WRITE - Fatal errorc'
Write (, '(a)') ' Could not open the output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') file_name
Stop
End Subroutine xy3_write
Subroutine xy6_write(file_name, node_num, node_x, node_y)
!*********************************************************************72
!
!c XY6_WRITE writes the coordinates of nodes of 6-node triangles to a file.
!
Implicit None
Integer node_num
Character (Len=*) file_name
Integer node
Double Precision node_x(node_num)
Double Precision node_y(node_num)
Open (Unit=1, File=file_name, Err=10)
Do node = 1, node_num
Write (1, *) node_x(node), node_y(node)
End Do
Close (Unit=1)
Write (, '(a)') ' '
Write (, '(a)') 'XY6_WRITE:'
Write (, '(a)') ' The 6-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)') 'XY6_WRITE - Fatal errorc'
Write (, '(a)') ' Could not open the output file:'
Write (, '(a)') ' '
Write (*, '(4x,a)') file_name
Stop
End Subroutine xy6_write
!c$$$Run, Output:
HCELL:
FORTRAN95,2013&18 version
Solve the Navier Stokes fluid flow
equations in an H-shaped region,
using finite elements.
The X direction is divided into three
regions, with element densities:
45 15 45
Corresponding NX = 106
The Y direction is divided into three
regions, with element densities:
5 1 5
Corresponding NY = 12
The X subregions are demarked by 4 values:
0.0000 45.0000 60.0000 105.0000
The Y subregions are demarked by 4 values:
0.0000 5.0000 6.0000 11.0000
Maximum number of nodes = 4853
Maximum number of elements = 2310
Maximum number of unknowns = 10978
HCELL_NODE_COUNT:
Number of nodes = 4673
HCELL_ELEMENT_COUNT:
Number of elements = 2130
HCELL_DOF_COUNT:
Number of degrees of freedom = 10618
ELEMENT_NODE_BANDWIDTH:
Lower half bandwidth = 106
Upper half bandwidth = 106
Total bandwidth = 213
Number of nodes = 4673
Number of elements = 2130
Number of unknowns = 10618
XY6_WRITE:
The 6-node triangle nodal coordinates were
written to the output file:
xy6.txt
XY3_WRITE:
The 3-node triangle nodal coordinates were
written to the output file:
xy3.txt
ELEMENT_NODE_WRITE:
The 3-node triangle nodal coordinates were
written to the output file:
element_node.txt
NODE_EPS:
An encapsulated Post file was created
containing an image of the nodes.
The file name is
hcell_nodes.eps
...Program finished with exit code 0
Press ENTER to exit console.
更多回帖
