program main !*****************************************************************************80 ! !! bump() solves a steady flow problem in a channel with a bump obstruction. ! ! Discussion: ! ! The Navier Stokes equations: ! ! The primitive variable formulation involves horizontal velocity U, ! vertical velocity V, and pressure P. The equations are: ! ! U dUdx + V dUdy + dPdx - mu*(ddU/dxdx + ddU/dydy) = F1 ! ! U dVdx + V dVdy + dPdy - mu*(ddV/dxdx + ddV/dydy) = F2 ! ! dUdx + dVdy = 0 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! John Burkardt, Max Gunzburger, Janet Peterson, ! Discretization of Cost and Sensitivities in Shape Optimization, ! in Computation and Control IV, ! edited by Bowers and Lund, ! Proceedings of the Fourth Bozeman Conference on Computation ! and Control, ! Birkhaeuser, 1995. ! ! Local Parameters: ! ! Local, real ( kind = 8 ) A(MAXROW,NEQN), contains the matrix ! in LINPACK general band storage mode. ! ! Local, real ( kind = 8 ) APROF, the value of the parameter for ! which the target profile was generated. ! ! Local, real ( kind = 8 ) AREA(NELEMN), the area of the elements. ! ! Local, real ( kind = 8 ) DCDA(MY), the sensitivities. ! ! Local, real ( kind = 8 ) F(NEQN), is used as a right hand side ! vector in LINSYS, and is overwritten there by the ! new solution (which is then copied into g). ! ! Local, real ( kind = 8 ) G(NEQN), the current solution vector, ! in which are stored pressures and velocities. ! ! Local, integer IBUMP, determines where isoparametric elements ! will be used. ! 0, no isoparametric elements will be used. ! 1, isoparametric elements will be used only for the ! elements which directly impinge on the bump. ! 2, isoparametric elements will be used for all elements which ! are above the bump. ! 3, isoparametric elements will be used for all elements. ! ! Local, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Local, integer INSC(NP), contains, for each node I, the ! index of the pressure P at that node, or 0. ! ! Local, integer ISOTRI(NELEMN), contains, for each element, ! a 1 if the element is isoparametric, or 0 otherwise. ! ! Local, integer IWRITE, controls the amount of output printed. ! ! Local, logical LONG. ! .TRUE. if region is "long and thin", and ! .FALSE. if region is "tall and skinny". ! ! Local, integer MAXELM, the maximum number of elements. ! This is actually simply the first dimension of the array NODE. ! ! Local, integer MAXEQN, the maximum number of equations, ! and functions. ! ! Local, integer MAXNEW, the maximum number of Newton steps ! to take in one iteration. ! ! Local, integer MAXNP, the maximum number of nodes. ! This is actually simply the first dimension of the array INDX. ! ! Local, integer MAXROW, the first dimension of the matrix A. ! ! Local, integer MAXSEC, the maximum number of secant steps ! to take. ! ! Local, integer MX, the number of nodes in the X direction. ! ! Local, integer MY, the number of nodes in the Y direction. ! ! Local, integer NBAND, the bandwidth of the linear system. ! ! Local, integer NBLEFT, the column at which the left corner ! of the bump lies. ! ! Local, integer NBRITE, the column at which the right corner ! of the bump lies. ! ! Local, integer NELEMN, the number of elements. ! ! Local, integer NEQN, the number of equations, and functions. ! ! Local, integer NLBAND, the lower bandwidth of the matrix. ! ! Local, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! ! Local, integer NODE(MAXELM,NNODES), contains, for each ! element, the global node indices of the element nodes. ! ! Local, integer NODEX0, the node whose Y coordinate is zero, ! and whose X coordinate is XPROF. This is the first ! node along the line where the velocity profile is ! measured. ! ! Local, integer NP, the number of nodes. ! ! Local, integer NQUAD, the number of quadrature points per ! element. ! ! Local, integer NROW, the number of rows need to store the ! matrix A. ! ! Local, integer NUMNEW, total number of Newton steps taken. ! ! Local, integer NX, controls the spacing of nodes and elements ! in the X direction. There are 2*NX+1 nodes in the X ! direction. ! ! Local, integer NY, controls the spacing of nodes and elements ! in the Y direction. There are 2*NY+1 nodes in the Y ! direction. ! ! Local, real ( kind = 8 ) PHI(NELEMN,NQUAD,NNODES,3). Each entry of PHI ! contains the value of a quadratic basis function or its derivative, ! evaluated at a quadrature point. ! In particular, PHI(I,J,K,1) is the value of the quadratic basis ! function associated with local node K in element I, evaluatated ! at quadrature point J. ! PHI(I,J,K,2) is the X derivative of that same basis function, ! PHI(I,J,K,3) is the Y derivative of that same basis function. ! ! Local, real ( kind = 8 ) PSI(NELEMN,NQUAD,NNODES). Each entry of PSI ! contains the value of a linear basis function evaluated at a ! quadrature point. ! PSI(I,J,K) is the value of the linear basis function associated ! with local node K in element I, evaluated at quadrature point J. ! ! Local, real ( kind = 8 ) RES(NEQN), holds the residual. ! ! Local, real ( kind = 8 ) REYNLD, the value of the Reynolds ! number. ! ! Local, real ( kind = 8 ) SENS(NEQN), the sensitivities. ! ! Local, real ( kind = 8 ) TOLNEW, convergence tolerance for Newton ! iteration. ! ! Local, real ( kind = 8 ) TOLSEC, convergence tolerance for secant ! iteration. ! ! Local, real ( kind = 8 ) UPROF(MY), the horizontal velocity ! along the profile line. ! ! Local, real ( kind = 8 ) XBLEFT, the left X coordinate of the bump. ! ! Local, real ( kind = 8 ) XBRITE, the right X coordinate of the bump. ! ! Local, real ( kind = 8 ) XC(NP), the X coordinates of the nodes. ! ! Local, real ( kind = 8 ) XLNGTH, the length of the region. ! ! Local, real ( kind = 8 ) XM(NELEMN,NQUAD), contains the X coordinates ! of the quadrature points for each element. ! ! Local, real ( kind = 8 ) XPROF, the X coordinate at which the ! profile is measured. This value should be a grid value! ! ! Local, integer XY_UNIT, the unit to which XY graphics ! is written. ! ! Local, real ( kind = 8 ) YC(NP), the Y coordinates of the nodes. ! ! Local, real ( kind = 8 ) YLNGTH, the height of the region. ! ! Local, real ( kind = 8 ) YM(NELEMN,NQUAD), contains the Y coordinates ! of the quadrature points for each element. ! implicit none integer, parameter :: maxnew = 4 integer, parameter :: maxsec = 10 integer, parameter :: nx = 21 integer, parameter :: ny = 7 ! ! This assignment should really read (maxrow = 27*min(nx,ny)) ! integer, parameter :: maxrow = 27*ny integer, parameter :: nelemn = 2*(nx-1)*(ny-1) integer, parameter :: mx = 2*nx-1 integer, parameter :: my = 2*ny-1 integer, parameter :: maxeqn = 2*mx*my+nx*ny integer, parameter :: np = mx*my integer, parameter :: nnodes = 6 integer, parameter :: nquad = 3 real ( kind = 8 ) a(maxrow,maxeqn) real ( kind = 8 ) anew real ( kind = 8 ) anext real ( kind = 8 ) aold real ( kind = 8 ) aprof real ( kind = 8 ) area(nelemn) real ( kind = 8 ) cpu1 real ( kind = 8 ) cpu2 real ( kind = 8 ) dcda(my) real ( kind = 8 ), external :: ddot real ( kind = 8 ) f(maxeqn) real ( kind = 8 ) g(maxeqn) real ( kind = 8 ) gr(my,my) integer i integer ibump integer iline(my) integer indx(np,2) integer insc(np) integer isotri(nelemn) integer iter integer itype ! integer iukk integer iwrite integer j logical long integer nband integer neqn integer nlband integer node(nelemn,nnodes) ! integer npara integer nrow integer numnew integer numsec ! real ( kind = 8 ) para real ( kind = 8 ) phi(nelemn,nquad,nnodes,3) real ( kind = 8 ) psi(nelemn,nquad,nnodes) real ( kind = 8 ) r(my) real ( kind = 8 ) res(maxeqn) real ( kind = 8 ) reynld real ( kind = 8 ) rjpnew real ( kind = 8 ) rjpold ! real ( kind = 8 ) rtemp(my) real ( kind = 8 ) sens(maxeqn) ! real tarray(2) real ( kind = 8 ) temp real ( kind = 8 ) test real ( kind = 8 ) tolnew real ( kind = 8 ) tolsec real ( kind = 8 ) uprof(my) character ( len = 255 ) uv_file integer uv_unit real ( kind = 8 ) xbleft real ( kind = 8 ) xbrite real ( kind = 8 ) xc(np) real ( kind = 8 ) xlngth real ( kind = 8 ) xm(nelemn,nquad) real ( kind = 8 ) xprof character ( len = 255 ) xy_file integer xy_unit real ( kind = 8 ) yc(np) real ( kind = 8 ) ylngth real ( kind = 8 ) ym(nelemn,nquad) real ( kind = 8 ) ypert call cpu_time ( cpu1 ) call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'bump():' write ( *, '(a)' ) ' Fortran90 version' write ( *, '(a)' ) ' Control problem for channel flow over a bump.' ! ! Set the input data. ! anew = 0.0D+00 anext = 0.3 aold = 0.0D+00 aprof = 0.25D+00 ibump = 2 iwrite = 1 numnew = 0 numsec = 0 reynld = 1.0D+00 rjpnew = 0.0D+00 rjpold = 0.0D+00 tolnew = 0.0001 tolsec = 0.0001 uv_file = 'uv_000.txt' xbleft = 1.0D+00 xbrite = 3.0D+00 xlngth = 10.0D+00 xprof = 4.0D+00 xy_file = 'xy_000.txt' ylngth = 3.0D+00 write ( *, * ) ' ' write ( *, * ) ' The bump will be generated with a height of ', aprof write ( *, * ) ' ' write ( *, * ) ' NX = ', nx write ( *, * ) ' NY = ', ny write ( *, * ) ' Number of elements = ', nelemn write ( *, * ) ' Reynolds number = ', reynld write ( *, * ) ' Secant tolerance = ', tolsec write ( *, * ) ' Newton tolerance = ', tolnew ! ! SETGRD constructs grid, numbers unknowns, calculates areas, ! and points for midpoint quadrature rule. ! call setgrd ( ibump, indx, insc, isotri, iwrite, long, maxeqn, mx, my, & nelemn, neqn, nnodes, node, np, nx, ny, xbleft, xbrite, xlngth ) ! ! Set the X Y coordinates of grid points. ! ypert = aprof call setxy ( iwrite, long, mx, my, np, nx, ny, xc, xlngth, yc, ylngth, ypert ) ! ! Set the quadrature points. ! call setqud ( area, isotri, iwrite, nelemn, nnodes, node, np, nquad, & xc, xm, yc, ym ) ! ! Set basis functions values at grid points ! call setbas ( isotri, nelemn, nnodes, node, np, nquad, phi, psi, xc, xm, & yc, ym ) ! ! Find points on velocity profile sampling line ! call setlin ( iline, indx, iwrite, long, mx, my, np, nx, ny, xlngth, xprof ) ! ! Compute the bandwidth. ! call setban ( indx, insc, maxrow, nband, nelemn, nlband, nnodes, node, & np, nrow ) ! ! Generate the solution of the Navier Stokes equations, using ! an initial estimate of G = 0. ! g(1:neqn) = 0.0D+00 call nstoke ( a, area, f, g, indx, insc, isotri, maxnew, & maxrow, nband, nelemn, neqn, nlband, nnodes, node, np, nquad, & nrow, numnew, phi, psi, reynld, tolnew, xc, xm, yc, ym ) ! ! Compute the residual for the solution G. ! call resid ( area, g, indx, insc, isotri, iwrite, nelemn, neqn, nnodes, & node, np, nquad, phi, psi, res, reynld, xc, xm, yc, ym ) ! ! Copy the flow along the profile line. ! call getg ( g, iline, my, neqn, uprof ) if ( 1 <= iwrite ) then write(*,*)' ' write(*,*)'Velocity profile:' write(*,*)' ' write(*,'(5g14.6)') (uprof(i),i = 1,my) end if ! ! Calculate Gram matrix GR and vector R = line integral of UPROF*phi ! call gram ( gr, iline, indx, iwrite, my, nelemn, nnodes, node, & np, r, uprof, xc, xprof, yc ) ! ! XY_WRITE creates a TABLE "XY" file ! call file_name_inc ( xy_file ) call get_unit ( xy_unit ) open ( unit = xy_unit, file = xy_file, status = 'replace' ) call xy_write ( xy_unit, np, xc, yc ) close ( unit = xy_unit ) call file_name_inc ( uv_file ) call get_unit ( uv_unit ) open ( unit = uv_unit, file = uv_file, status = 'replace' ) call uv_write ( f, indx, uv_unit, neqn, np, yc ) close ( unit = uv_unit ) ! ! Destroy information about true solution before beginning ! secant iteration. ! g(1:neqn) = 0.0D+00 ! ! Secant iteration loop: ! do iter = 1, maxsec write(*,*) ' ' write(*,*) 'Secant iteration ',iter numsec = numsec + 1 ! ! Update the grid. ! ypert = anew call setxy ( iwrite, long, mx, my, np, nx, ny, xc, xlngth, yc, & ylngth, ypert ) ! ! Set the quadrature points. ! call setqud ( area, isotri, iwrite, nelemn, nnodes, node, np, nquad, & xc, xm, yc, ym ) ! ! Set basis functions values at grid points. ! call setbas(isotri,nelemn,nnodes,node,np,nquad,phi,psi,xc,xm,yc,ym) ! ! Solve for the flow at the new value of the parameter. ! call nstoke(a,area,f,g,indx,insc,isotri,maxnew, & maxrow,nband,nelemn,neqn,nlband,nnodes,node,np,nquad, & nrow,numnew,phi,psi,reynld,tolnew,xc,xm,yc,ym) ! ! Get the velocity profile UPROF along sampling line. ! call getg ( g, iline, my, neqn, uprof ) if ( 1 <= iwrite ) then write ( *, * ) ' ' write ( *, * ) 'Velocity profile:' write ( *, * ) ' ' write ( *, '(5g14.6)' ) uprof(1:my) end if ! ! Solve linear system for sensitivities. ! itype = -2 call linsys ( a, area, sens, g, indx, insc, isotri, itype, maxrow, & nband, nelemn, neqn, nlband, nnodes, node, np, nquad, nrow, & phi, psi, reynld, xc, xm, yc, ym ) ! ! Get the sensitivities DCDA along sampling line. ! call getg ( sens, iline, my, neqn, dcda ) if ( 2 <= iwrite ) then write(*,*)' ' write(*,*)'Sensitivities:' write(*,*)' ' write(*,'(5g14.6)') dcda(1:my) end if ! ! Evaluate J prime at current value of parameter where J is ! the functional to be minimized. ! ! JPRIME = 2.0D+00 * DCDA(I) * ( GR(I,J) * UPROF(J) - R(I) ) ! rjpnew = 0.0D+00 do i = 1, my temp = -r(i) do j = 1, my temp = temp + gr(i,j) * uprof(j) end do rjpnew = rjpnew + 2.0D+00 * dcda(i) * temp end do ! ! Write data to graphics files. ! call file_name_inc ( xy_file ) call get_unit ( xy_unit ) open ( unit = xy_unit, file = xy_file, status = 'replace' ) call xy_write ( xy_unit, np, xc, yc ) close ( unit = xy_unit ) call file_name_inc ( uv_file ) open ( unit = uv_unit, file = uv_file, status = 'replace' ) call uv_write ( f, indx, uv_unit, neqn, np, yc ) close ( unit = uv_unit ) write ( *, * ) ' ' write ( *, * ) ' Parameter = ', anew, ' J prime=', rjpnew ! ! Update the estimate of the parameter using the secant step. ! if ( 1 < iter ) then anext = aold - rjpold * ( anew - aold ) / ( rjpnew - rjpold ) end if aold = anew anew = anext rjpold = rjpnew if ( anew /= 0.0D+00 ) then test = abs ( anew - aold ) / anew else test = 0.0D+00 end if write ( *, * ) ' New value of parameter = ', anew write ( *, * ) ' Convergence test = ', test if ( abs ( ( anew - aold ) ) <= abs ( anew ) * tolsec .and. 1 < iter ) then write ( *, * ) 'Secant iteration converged.' go to 40 end if end do write ( *, * ) ' Secant iteration failed to converge.' 40 continue ! ! Produce total CPU time used. ! call cpu_time ( cpu2 ) write ( *, '(a)' ) ' ' write ( *, * ) ' Total execution time = ', cpu2 - cpu1, ' seconds.' write ( *, * ) ' Number of secant steps = ', numsec write ( *, * ) ' Number of Newton steps = ', numnew ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'bump():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end function bsp ( it, iq, id, nelemn, nnodes, node, np, xc, xq, yc, yq ) !*****************************************************************************80 ! !! BSP evaluates the linear basis function associated with pressure. ! ! Discussion: ! ! The local reference element: ! ! ^ ! | ! 1 3 ! | |\ ! | | \ ! | | \ ! | | \ ! 0 1----2 ! | ! +--0----1----> ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 01 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IT, the element in which the basis function ! is defined. ! ! Input, integer IQ, the index (1, 2 or 3) of the local node at ! which the basis function is to be evaluated. ! ! Input, integer ID, the index (1, 2 or 3) of the local node ! associated with the basis function. ! implicit none integer nelemn integer nnodes integer np real ( kind = 8 ) bsp real ( kind = 8 ) d integer g1 integer g2 integer g3 integer i4_wrap integer id integer iq integer it integer l1 integer l2 integer l3 integer node(nelemn,nnodes) real ( kind = 8 ) xc(np) real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) yq ! ! L1, L2, L3 are the appropriate local node indices. ! l1 = iq l2 = i4_wrap ( iq + 1, 1, 3 ) l3 = i4_wrap ( iq + 2, 1, 3 ) ! ! G1, G2, G3 are the global node indices. ! g1 = node(it,l1) g2 = node(it,l2) g3 = node(it,l3) d = (xc(g2)-xc(g1))*(yc(g3)-yc(g1))-(xc(g3)-xc(g1))*(yc(g2)-yc(g1)) if ( id == 1 ) then bsp = 1.0D+00 +((yc(g2)-yc(g3))*(xq-xc(g1))+(xc(g3)-xc(g2))*(yq-yc(g1)))/d else if ( id == 2 ) then bsp = (yc(g2)-yc(g3))/d else if ( id == 3 ) then bsp = (xc(g3)-xc(g2))/d else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BSP - Fatal error!' write ( *, '(a)' ) ' Illegal local index value for linear basis.' write ( *, '(a)' ) ' Legal values are 1, 2 or 3.' write ( *, '(a,i6)' ) ' The input value was ID = ', id stop end if return end subroutine daxpy ( n, da, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DAXPY computes constant times a vector plus a vector. ! ! Discussion: ! ! Uses unrolled loops for increments equal to one. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! Jack Dongarra ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of elements in DX and DY. ! ! Input, real ( kind = 8 ) DA, the multiplier of DX. ! ! Input, real ( kind = 8 ) DX(*), the first vector. ! ! Input, integer INCX, the increment between successive ! entries of DX. ! ! Input/output, real ( kind = 8 ) DY(*), the second vector. ! On output, DY(*) has been replaced by DY(*) + DA * DX(*). ! ! Input, integer INCY, the increment between successive ! entries of DY. ! implicit none real ( kind = 8 ) da real ( kind = 8 ) dx(*) real ( kind = 8 ) dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n if ( n <= 0 ) then return end if if ( da == 0.0D+00 ) then return end if ! ! Code for unequal increments or equal increments ! not equal to 1. ! if ( incx /= 1 .or. incy /= 1 ) then if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n dy(iy) = dy(iy) + da * dx(ix) ix = ix + incx iy = iy + incy end do ! ! Code for both increments equal to 1. ! else m = mod ( n, 4 ) do i = 1, m dy(i) = dy(i) + da * dx(i) end do do i = m+1, 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 end if return end function ddot ( n, dx, incx, dy, incy ) !*****************************************************************************80 ! !! DDOT forms the dot product of two vectors. ! ! Discussion: ! ! This routine uses unrolled loops for increments equal to one. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! Jack Dongarra ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vectors. ! ! Input, real ( kind = 8 ) DX(*), the first vector. ! ! Input, integer INCX, the increment between successive ! entries in DX. ! ! Input, real ( kind = 8 ) DY(*), the second vector. ! ! Input, integer INCY, the increment between successive ! entries in DY. ! ! Output, real ( kind = 8 ) DDOT, the sum of the product of the ! corresponding entries of DX and DY. ! implicit none real ( kind = 8 ) ddot real ( kind = 8 ) dtemp real ( kind = 8 ) dx(*) real ( kind = 8 ) dy(*) integer i integer incx integer incy integer ix integer iy integer m integer n ddot = 0.0D+00 dtemp = 0.0D+00 if ( n <= 0 ) then return end if ! ! Code for unequal increments or equal increments ! not equal to 1. ! if ( incx /= 1 .or. incy /= 1 ) then if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1, n dtemp = dtemp + dx(ix) * dy(iy) ix = ix + incx iy = iy + incy end do ! ! Code for both increments equal to 1. ! else m = mod ( n, 5 ) do i = 1, m dtemp = dtemp + dx(i) * dy(i) end do do i = m+1, 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 end if ddot = dtemp return end subroutine dgbfa ( abd, lda, n, ml, mu, ipvt, info ) !*****************************************************************************80 ! !! DGBFA factors a real band matrix by elimination. ! ! Discussion: ! ! DGBFA is usually called by DGBCO, but it can be called ! directly with a saving in time if RCOND is not needed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! FORTRAN90 translation by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input/output, real ( kind = 8 ) ABD(LDA,N). On input, the matrix in band ! storage. The columns of the matrix are stored in the columns of ABD ! and the diagonals of the matrix are stored in rows ML+1 through ! 2*ML+MU+1 of ABD. On output, an upper triangular matrix in band storage ! and the multipliers which were used to obtain it. The factorization ! can be written A = L*U where L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! Input, integer LDA, the leading dimension of the array ABD. ! 2*ML + MU + 1 <= LDA is required. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, MU, the number of diagonals below and ! above the main diagonal. 0 <= ML < N, 0 <= MU < N. ! ! Output, integer IPVT(N), the pivot indices. ! ! Output, integer INFO, error flag. ! 0, normal value. ! K, if U(K,K) == 0.0D+00. This is not an error condition for this ! subroutine, but it does indicate that DGBSL will divide by zero if ! called. Use RCOND in DGBCO for a reliable indication of singularity. ! implicit none integer lda integer n real ( kind = 8 ) abd(lda,n) integer i integer i0 integer info integer ipvt(n) integer idamax integer j integer j0 integer j1 integer ju integer jz integer k integer l integer lm integer m integer ml integer mm integer mu real ( kind = 8 ) t m = ml + mu + 1 info = 0 ! ! Zero initial fill-in columns. ! j0 = mu + 2 j1 = min ( n, m ) - 1 do jz = j0, j1 i0 = m + 1 - jz do i = i0, ml abd(i,jz) = 0.0D+00 end do end do jz = j1 ju = 0 ! ! Gaussian elimination with partial pivoting. ! do k = 1, n-1 ! ! Zero out the next fill-in column. ! jz = jz + 1 if ( jz <= n ) then abd(1:ml,jz) = 0.0D+00 end if ! ! Find L = pivot index. ! lm = min ( 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.0D+00 ) then info = k ! ! Interchange if necessary. ! else if ( l /= m ) then t = abd(l,k) abd(l,k) = abd(m,k) abd(m,k) = t end if ! ! Compute multipliers. ! t = -1.0D+00 / abd(m,k) call dscal ( lm, t, abd(m+1,k), 1 ) ! ! Row elimination with column indexing. ! ju = min ( max ( ju, mu+ipvt(k) ), n ) mm = m do j = k+1, ju l = l - 1 mm = mm - 1 t = abd(l,j) if ( l /= mm ) then abd(l,j) = abd(mm,j) abd(mm,j) = t end if call daxpy ( lm, t, abd(m+1,k), 1, abd(mm+1,j), 1 ) end do end if end do ipvt(n) = n if ( abd(m,n) == 0.0D+00 ) then info = n end if return end subroutine dgbsl ( abd, lda, n, ml, mu, ipvt, b, job ) !*****************************************************************************80 ! !! DGBSL solves a real banded system factored by DGBCO or DGBFA. ! ! Discussion: ! ! DGBSL can solve either A * X = B or A' * X = B. ! ! A division by zero will occur if the input factor contains a ! zero on the diagonal. Technically this indicates singularity ! but it is often caused by improper arguments or improper ! setting of LDA. It will not occur if the subroutines are ! called correctly and if DGBCO has set 0.0 < RCOND ! or DGBFA has set INFO == 0. ! ! To compute inverse(A) * C where C is a matrix with P columns: ! ! call dgbco ( abd, lda, n, ml, mu, ipvt, rcond, z ) ! ! if ( rcond is too small ) then ! exit ! end if ! ! do j = 1, p ! call dgbsl ( abd, lda, n, ml, mu, ipvt, c(1,j), 0 ) ! end do ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! FORTRAN90 translation by John Burkardt. ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Parameters: ! ! Input, real ( kind = 8 ) ABD(LDA,N), the output from DGBCO or DGBFA. ! ! Input, integer LDA, the leading dimension of the array ABD. ! ! Input, integer N, the order of the matrix. ! ! Input, integer ML, MU, the number of diagonals below and ! above the main diagonal. 0 <= ML < N, 0 <= MU < N. ! ! Input, integer IPVT(N), the pivot vector from DGBCO or DGBFA. ! ! Input/output, real ( kind = 8 ) B(N). On input, the right hand side. ! On output, the solution. ! ! Input, integer JOB, job choice. ! 0, solve A*X=B. ! nonzero, solve A'*X=B. ! implicit none integer lda integer n real ( kind = 8 ) abd(lda,n) real ( kind = 8 ) b(n) integer ipvt(n) integer job integer k integer l integer la integer lb integer lm integer m integer ml integer mu real ( kind = 8 ) ddot real ( kind = 8 ) t m = mu + ml + 1 ! ! JOB = 0, Solve A * x = b. ! ! First solve L * y = b. ! if ( job == 0 ) then if ( 0 < ml ) then do k = 1, n-1 lm = min ( ml, n-k ) l = ipvt(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if call daxpy ( lm, t, abd(m+1,k), 1, b(k+1), 1 ) end do end if ! ! Now solve U * x = y. ! do k = n, 1, -1 b(k) = b(k) / abd(m,k) lm = min ( 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 ! ! JOB nonzero, solve A' * x = b. ! ! First solve U' * y = b. ! else do k = 1, n lm = min ( 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 L' * x = y. ! if ( 0 < ml ) then do k = n-1, 1, -1 lm = min ( ml, n-k ) b(k) = b(k) + ddot ( lm, abd(m+1,k), 1, b(k+1), 1 ) l = ipvt(k) if ( l /= k ) then t = b(l) b(l) = b(k) b(k) = t end if end do end if end if return end subroutine dscal ( n, sa, x, incx ) !*****************************************************************************80 ! !! DSCAL scales a vector by a constant. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 1999 ! ! Author: ! ! Jack Dongarra ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch and Pete Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real ( kind = 8 ) SA, the multiplier. ! ! Input/output, real ( kind = 8 ) X(*), the vector to be scaled. ! ! Input, integer INCX, the increment between successive ! entries of X. ! implicit none integer i integer incx integer ix integer m integer n real ( kind = 8 ) sa real ( kind = 8 ) x(*) if ( n <= 0 ) then else if ( incx == 1 ) then m = mod ( n, 5 ) x(1:m) = sa * x(1:m) do i = m+1, n, 5 x(i) = sa * x(i) x(i+1) = sa * x(i+1) x(i+2) = sa * x(i+2) x(i+3) = sa * x(i+3) x(i+4) = sa * x(i+4) end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if do i = 1, n x(ix) = sa * x(ix) ix = ix + incx end do end if return end subroutine file_name_inc ( file_name ) !*****************************************************************************80 ! !! FILE_NAME_INC increments a partially numeric filename. ! ! Discussion: ! ! It is assumed that the digits in the name, whether scattered or ! connected, represent a number that is to be increased by 1 on ! each call. If this number is all 9's on input, the output number ! is all 0's. Non-numeric letters of the name are unaffected. ! ! If the name is empty, then the routine stops. ! ! If the name contains no digits, the empty string is returned. ! ! Example: ! ! Input Output ! ----- ------ ! 'a7to11.txt' 'a7to12.txt' ! 'a7to99.txt' 'a8to00.txt' ! 'a9to99.txt' 'a0to00.txt' ! 'cat.txt' ' ' ! ' ' STOP! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) FILE_NAME. ! On input, a character string to be incremented. ! On output, the incremented string. ! implicit none character c integer change integer digit character ( len = * ) file_name integer i integer lens lens = len_trim ( file_name ) if ( lens <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_NAME_INC - Fatal error!' write ( *, '(a)' ) ' The input string is empty.' stop end if change = 0 do i = lens, 1, -1 c = file_name(i:i) if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then change = change + 1 digit = ichar ( c ) - 48 digit = digit + 1 if ( digit == 10 ) then digit = 0 end if c = char ( digit + 48 ) file_name(i:i) = c if ( c /= '0' ) then return end if end if end do if ( change == 0 ) then file_name = ' ' return end if return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT, the free unit number. ! implicit none integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine getg ( f, iline, my, neqn, u ) !*****************************************************************************80 ! !! GETG extracts the values of a quantity along the profile line. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) F(NEQN), the finite element coefficients. ! ! Input, integer ILINE(MY), the index of the finite element ! coefficient associated with each node on the profile. ! ! Input, integer MY, the number of nodes along the profile line. ! ! Input, integer NEQN, the number of finite element ! coefficients. ! ! Output, real ( kind = 8 ) U(MY), the finite element coefficients ! associated with each node of the profile. ! implicit none integer my integer neqn real ( kind = 8 ) f(neqn) integer i integer iline(my) integer j real ( kind = 8 ) u(my) do i = 1, my u(i) = 0.0D+00 j = iline(i) if ( j <= 0 ) then u(i) = 0.0D+00 else u(i) = f(j) end if end do return end subroutine gram ( gr, iline, indx, iwrite, my, nelemn, nnodes, node, & np, r, uprof, xc, xprof, yc ) !*****************************************************************************80 ! !! GRAM computes and stores the Gram matrix. ! ! Discussion: ! ! The routine computes the Gram matrix and the vector ! whose components are the line integral of UI*phi(j). ! ! A three-point Gauss quadrature rule is used to evaluate the line integral. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) GR(MY,MY), the Gram matrix. ! ! Input, integer ILINE(MY), the list of global unknown numbers ! along the profile line. ! ! Input, integer INDX(NP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer IWRITE, controls the amount of output printed. ! ! Input, integer MY, the number of nodes along the profile line. ! ! Input, integer NELEMN, the number of elements. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! ! Input, integer NODE(NELEMN,NNODES), contains, for each ! element, the global node indices of the element nodes. ! ! Input, integer NP, the number of nodes. ! ! Output, real ( kind = 8 ) R(MY), the line integral of UPROF * PHI. ! ! Input, real ( kind = 8 ) UPROF(MY), the horizontal velocity ! along the profile line. ! ! Input, real ( kind = 8 ) XC(NP), the X coordinates of the nodes. ! ! Input, real ( kind = 8 ) XPROF, the X coordinate at which the ! profile is measured. This value should be a grid value! ! ! Input, real ( kind = 8 ) YC(NP), the Y coordinates of the nodes. ! implicit none integer my integer nelemn integer nnodes integer np real ( kind = 8 ) ar real ( kind = 8 ) bb real ( kind = 8 ) bbb real ( kind = 8 ) bbx real ( kind = 8 ) bby real ( kind = 8 ) bma2 real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) gr(my,my) integer i integer igetl integer ii integer iline(my) integer indx(np,2) integer ip integer ipp integer iq integer iqq integer iquad integer it integer iun integer iwrite integer j integer jj integer k integer kk integer node(nelemn,nnodes) real ( kind = 8 ) r(my) real ( kind = 8 ) ubc real ( kind = 8 ) ubdry real ( kind = 8 ) uiqdpt real ( kind = 8 ) uprof(my) real ( kind = 8 ) wt(3) real ( kind = 8 ) x real ( kind = 8 ) xc(np) real ( kind = 8 ) xprof real ( kind = 8 ) y real ( kind = 8 ) yc(np) real ( kind = 8 ) yq(3) ! ! Values for 3 point Gauss quadrature. ! wt(1) = 5.0D+00 / 9.0D+00 wt(2) = 8.0D+00 / 9.0D+00 wt(3) = 5.0D+00 / 9.0D+00 yq(1) = -0.7745966692D+00 yq(2) = 0.0D+00 yq(3) = 0.7745966692D+00 r(1:my) = 0.0D+00 gr(1:my,1:my) = 0.0D+00 ! ! Compute the line integral by looping over intervals along line ! using three point Gauss quadrature ! do it = 1, nelemn ! ! Check to see if we are in a triangle with a side along line ! x = xprof ! k = node(it,1) kk = node(it,2) if ( 1.0D-04 < abs ( xc(k) - xprof ) .or. & 1.0D-04 < abs ( xc(kk) - xprof ) ) then go to 70 end if do iquad = 1, 3 bma2 = ( yc(kk) - yc(k) ) / 2.0D+00 ar = bma2 * wt(iquad) x = xprof y = yc(k) + bma2 * ( yq(iquad) + 1.0D+00 ) ! ! Compute U internal at quadrature points. ! uiqdpt = 0.0D+00 do iq = 1, nnodes if ( iq == 1 .or. iq == 2 .or. iq == 4 ) then call qbf ( x, y, it, iq, bb, bx, by, nelemn, nnodes, & node, np, xc, yc ) ip = node(it,iq) iun = indx(ip,1) if ( 0 < iun ) then ii = igetl(iun,iline,my) uiqdpt = uiqdpt + bb * uprof(ii) else if ( iun == -1) then ubc = ubdry(1,yc(ip)) uiqdpt = uiqdpt + bb * ubc end if end if end do ! ! Only loop over nodes lying on line x = xprof ! do iq = 1, nnodes if ( iq == 1 .or. iq == 2 .or. iq == 4 ) then ip = node(it,iq) call qbf(x,y,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) i = indx(ip,1) if ( 0 < i ) then ii = igetl(i,iline,my) r(ii) = r(ii) + bb * uiqdpt * ar do iqq = 1, nnodes if ( iqq == 1 .or. iqq == 2 .or. iqq == 4 ) then ipp = node(it,iqq) call qbf(x,y,it,iqq,bbb,bbx,bby,nelemn,nnodes,node,np,xc,yc) j = indx(ipp,1) if ( j /= 0 ) then jj = igetl(j,iline,my) gr(ii,jj) = gr(ii,jj)+bb*bbb*ar end if end if end do end if end if end do end do 70 continue end do if ( 3 <= iwrite ) then write(*,*)' ' write(*,*)'Gram matrix:' write(*,*)' ' do i = 1, my do j = 1, my write(*,*)i,j,gr(i,j) end do end do write(*,*)' ' write(*,*)'R vector:' write(*,*)' ' do i = 1, my write(*,*)r(i) end do end if return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of I4 division. ! ! Discussion: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! The MOD function computes a result with the same sign as the ! quantity being divided. Thus, suppose you had an angle A, ! and you wanted to ensure that it was between 0 and 360. ! Then mod(A,360) would do, if A was positive, but if A ! was negative, your result would be between -360 and 0. ! ! On the other hand, I4_MODP(A,360) is between 0 and 360, always. ! ! Example: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number to be divided. ! ! Input, integer J, the number that divides I. ! ! Output, integer I4_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none integer i integer i4_modp integer j integer value if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i8)' ) ' Illegal divisor J = ', j stop end if value = mod ( i, j ) if ( value < 0 ) then value = value + abs ( j ) end if i4_modp = value return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an I4 to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I Value ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IVAL, an integer value. ! ! Input, integer ILO, IHI, the desired bounds for the value. ! ! Output, integer I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer i4_modp integer i4_wrap integer ihi integer ilo integer ival integer jhi integer jlo integer value integer wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then value = jlo else value = jlo + i4_modp ( ival - jlo, wide ) end if i4_wrap = value return end function idamax ( n, dx, incx ) !*****************************************************************************80 ! !! IDAMAX finds the index of the vector element of maximum absolute value. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 April 1999 ! ! Reference: ! ! Jack Dongarra, Cleve Moler, Jim Bunch, Pete Stewart, ! LINPACK User's Guide, ! SIAM, (Society for Industrial and Applied Mathematics), ! 3600 University City Science Center, ! Philadelphia, PA, 19104-2688. ! ISBN 0-89871-172-X ! ! Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Basic Linear Algebra Subprograms for Fortran Usage, ! Algorithm 539, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! ! Parameters: ! ! Input, integer N, the number of entries in the vector. ! ! Input, real ( kind = 8 ) X(*), the vector to be examined. ! ! Input, integer INCX, the increment between successive ! entries of SX. ! ! Output, integer IDAMAX, the index of the element of SX of ! maximum absolute value. ! implicit none real ( kind = 8 ) dmax real ( kind = 8 ) dx(*) integer i integer idamax integer incx integer ix integer n idamax = 0 if ( n < 1 .or. incx <= 0 ) then return end if idamax = 1 if ( n == 1 ) then return end if if ( incx == 1 ) then dmax = abs ( dx(1) ) do i = 2, n if ( dmax < abs ( dx(i) ) ) then idamax = i dmax = abs ( dx(i) ) end if end do else ix = 1 dmax = abs ( dx(1) ) ix = ix + incx do i = 2, n if ( dmax < abs ( dx(ix) ) ) then idamax = i dmax = abs ( dx(ix) ) end if ix = ix + incx end do end if return end function igetl ( i, iline, my ) !*****************************************************************************80 ! !! IGETL gets the local unknown number along the profile line. ! ! Discussion: ! ! For our problem, the profile line is specified by X = XPROF. ! We are given a global unknown number, and need to determine the ! local unknown number. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the global unknown number. ! ! Input, integer ILINE(MY), the list of global unknown numbers ! along the profile line. ! ! Input, integer MY, the number of nodes on the profile line. ! ! Output, integer IGETL, the index of the node in the profile ! line at which global unknown I occurs, or -1 if there is no such entry. ! implicit none integer my integer i integer igetl integer iline(my) integer j igetl = -1 do j = 1, my if ( iline(j) == i ) then igetl = j return end if end do return end subroutine linsys ( a, area, f, g, indx, insc, isotri, itype, maxrow, & nband, nelemn, neqn, nlband, nnodes, node, np, nquad, nrow, & phi, psi, reynld, xc, xm, yc, ym ) !*****************************************************************************80 ! !! LINSYS solves the linearized Navier Stokes equation. ! ! Discussion: ! ! ITYPE = -1 for solutions of the Navier Stokes equation. ! ITYPE = -2 for solutions of the sensitivity equation. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) AREA(NELEMN), the area of the elements. ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! ! Local Parameters: ! ! Local, integer IPIVOT(NEQN), pivoting space. ! implicit none integer maxrow integer nelemn integer neqn integer nnodes integer np integer nquad real ( kind = 8 ) a(maxrow,neqn) real ( kind = 8 ) aij real ( kind = 8 ) ar real ( kind = 8 ) area(nelemn) real ( kind = 8 ) bb real ( kind = 8 ) bbb real ( kind = 8 ) bbbl real ( kind = 8 ) bbl real ( kind = 8 ) bbx real ( kind = 8 ) bby real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) det real ( kind = 8 ) etax real ( kind = 8 ) etay real ( kind = 8 ) f(neqn) real ( kind = 8 ) g(neqn) integer i integer ihor integer indx(np,2) integer info integer insc(np) integer ioff integer ip integer ipivot(neqn) integer ipp integer iprs integer iq integer iqq integer iquad integer isotri(nelemn) integer it integer itype integer iuse integer iver integer j integer job integer jp integer ju ! integer juse integer jv integer nband integer nlband integer node(nelemn,nnodes) integer nrow real ( kind = 8 ) phi(nelemn,nquad,nnodes,3) real ( kind = 8 ) psi(nelemn,nquad,nnodes) real ( kind = 8 ) reynld real ( kind = 8 ) ubc real ( kind = 8 ) ubdry real ( kind = 8 ) ubump real ( kind = 8 ) un(2) real ( kind = 8 ) unx(2) real ( kind = 8 ) uny(2) real ( kind = 8 ) visc real ( kind = 8 ) xc(np) real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) xm(nelemn,nquad) real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) ym(nelemn,nquad) real ( kind = 8 ) yq ioff = nlband + nlband + 1 visc = 1.0D+00 / reynld f(1:neqn) = 0.0D+00 a(1:nrow,1:neqn) = 0.0D+00 do it = 1, nelemn ar = area(it) / 3.0D+00 do iquad = 1, nquad yq = ym(it,iquad) xq = xm(it,iquad) if ( isotri(it) == 1 ) then call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) ar = det * area(it) / 3.0D+00 end if call uval ( etax, etay, g, indx, isotri, it, nelemn, neqn, & nnodes, node, np, un, uny, unx, xc, xix, xiy, xq, yc, yq ) ! ! For each basis function: ! do iq = 1, nnodes ip = node(it,iq) bb = phi(it,iquad,iq,1) bx = phi(it,iquad,iq,2) by = phi(it,iquad,iq,3) bbl = psi(it,iquad,iq) ihor = indx(ip,1) iver = indx(ip,2) iprs = insc(ip) if ( 0 < ihor ) then f(ihor) = f(ihor) + ar * bb*(un(1)*unx(1)+un(2)*uny(1)) end if if ( 0 < iver ) then f(iver) = f(iver) + ar * bb*(un(1)*unx(2)+un(2)*uny(2)) end if ! ! For another basis function, ! do iqq = 1, nnodes ipp = node(it,iqq) bbb = phi(it,iquad,iqq,1) bbx = phi(it,iquad,iqq,2) bby = phi(it,iquad,iqq,3) bbbl = psi(it,iquad,iqq) ju = indx(ipp,1) jv = indx(ipp,2) jp = insc(ipp) ! ! Horizontal velocity variable ! if ( 0 < ju ) then if ( 0 < ihor ) then iuse = ihor-ju+ioff a(iuse,ju) = a(iuse,ju)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))) end if if ( 0 < iver ) then iuse = iver-ju+ioff a(iuse,ju) = a(iuse,ju)+ar*bb*bbb*unx(2) end if if ( 0 < iprs ) then iuse = iprs-ju+ioff a(iuse,ju) = a(iuse,ju)+ar*bbx*bbl end if else if ( ju == itype ) then if ( ju == -1 ) then ubc = ubdry(1,yc(ipp)) else if ( ju == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,1,nelemn, & neqn,nnodes,node,np,xc,yc) end if if ( 0 < ihor ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))) f(ihor) = f(ihor)-ubc*aij end if if ( 0 < iver ) then aij = ar*bb*bbb*unx(2) f(iver) = f(iver)-ubc*aij end if if ( 0 < iprs ) then aij = ar*bbx*bbl f(iprs) = f(iprs)-ubc*aij end if end if ! ! Vertical velocity variable ! if ( 0 < jv ) then if ( 0 < ihor ) then iuse = ihor-jv+ioff a(iuse,jv) = a(iuse,jv)+ar*bb*bbb*uny(1) end if if ( 0 < iver ) then iuse = iver-jv+ioff a(iuse,jv) = a(iuse,jv)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))) end if if ( 0 < iprs ) then iuse = iprs-jv+ioff a(iuse,jv) = a(iuse,jv)+ar*bby*bbl end if else if ( jv == itype ) then if ( jv == -1 ) then ubc = ubdry(2,yc(ipp)) else if ( jv == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,2,nelemn, & neqn,nnodes,node,np,xc,yc) end if if ( 0 < ihor ) then aij = ar*bb*bbb*uny(1) f(ihor) = f(ihor)-ubc*aij end if if ( 0 < iver ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))) f(iver) = f(iver)-ubc*aij end if if ( 0 < iprs ) then aij = ar * bby * bbl f(iprs) = f(iprs) - ubc * aij end if end if ! ! Pressure variable ! if ( 0 < jp ) then if ( 0 < ihor ) then iuse = ihor-jp+ioff a(iuse,jp) = a(iuse,jp) - ar * bx * bbbl end if if ( 0 < iver ) then iuse = iver-jp+ioff a(iuse,jp) = a(iuse,jp) - ar * by * bbbl end if end if end do end do end do end do ! ! The last equation is "reset" to require that the last pressure ! be zero. ! f(neqn) = 0.0D+00 do j = neqn-nlband, neqn-1 i = neqn-j+ioff a(i,j) = 0.0D+00 end do a(nband,neqn) = 1.0D+00 ! ! Factor the matrix ! call dgbfa ( a, maxrow, neqn, nlband, nlband, ipivot, info ) if (info /= 0) then write (*,*) ' ' write (*,*) 'LINSYS - fatal error!' write (*,*) 'DGBFA returns INFO = ',info stop end if ! ! Solve the linear system ! job = 0 call dgbsl ( a, maxrow, neqn, nlband, nlband, ipivot, f, job ) return end subroutine nstoke ( a, area, f, g, indx, insc, isotri, maxnew, & maxrow, nband, nelemn, neqn, nlband, nnodes, node, np, nquad, & nrow, numnew, phi, psi, reynld, tolnew, xc, xm, yc, ym ) !*****************************************************************************80 ! !! NSTOKE solves the Navier Stokes equation using Taylor-Hood elements. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) AREA(NELEMN), the area of the elements. ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer maxrow integer nelemn integer neqn integer nnodes integer np integer nquad real ( kind = 8 ) a(maxrow,neqn) real ( kind = 8 ) area(nelemn) real ( kind = 8 ) diff real ( kind = 8 ) f(neqn) real ( kind = 8 ) g(neqn) integer idamax integer indx(np,2) integer insc(np) integer isotri(nelemn) integer iter integer itype integer maxnew integer nband integer nlband integer node(nelemn,nnodes) integer nrow integer numnew real ( kind = 8 ) phi(nelemn,nquad,nnodes,3) real ( kind = 8 ) psi(nelemn,nquad,nnodes) real ( kind = 8 ) reynld real ( kind = 8 ) tolnew real ( kind = 8 ) xc(np) real ( kind = 8 ) xm(nelemn,nquad) real ( kind = 8 ) yc(np) real ( kind = 8 ) ym(nelemn,nquad) ! ! G contains an initial estimate of the solution. ! do iter = 1, maxnew numnew = numnew + 1 itype = -1 call linsys ( a, area, f, g, indx, insc, isotri, itype, maxrow, & nband, nelemn, neqn, nlband, nnodes, node, np, nquad, nrow, & phi, psi, reynld, xc, xm, yc, ym ) ! ! Check for convergence. ! g(1:neqn) = g(1:neqn) - f(1:neqn) diff = abs ( g(idamax(neqn,g,1)) ) write(*,*) 'NSTOKE: Iteration ', iter, ' MaxNorm(diff) = ', diff g(1:neqn) = f(1:neqn) if ( diff <= tolnew ) then write(*,*) 'NSTOKE converged.' exit end if if ( iter == maxnew ) then write(*,*) 'NSTOKE failed!' stop end if end do return end subroutine pval ( g, insc, long, mx, my, nelemn, neqn, nnodes, node, & np, press ) !*****************************************************************************80 ! !! PVAL computes a table of pressures at all the nodes. ! ! Discussion: ! ! This table is needed in order to prepare output data for PLOT3D. ! ! This routine does not handle general isoparametric elements ! accurately. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MX, the number of nodes in the X direction. ! ! Input, integer MY, the number of nodes in the Y direction. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer mx integer my integer nelemn integer neqn integer nnodes integer np real ( kind = 8 ) g(neqn) integer i integer insc(np) integer ip integer iq integer it integer ivar integer j logical long integer node(nelemn,nnodes) real ( kind = 8 ) press(mx,my) press(1:mx,1:my) = 0.0D+00 ! ! Read the pressures where they are computed. ! These are "(odd, odd)" points. ! do it = 1, nelemn do iq = 1, 3 ip = node(it,iq) ivar = insc(ip) if ( long ) then i = ((ip-1)/my)+1 j = mod(ip-1,my)+1 else i = mod(ip-1,mx)+1 j = ((ip-1)/mx)+1 end if if ( 0 < ivar ) then press(i,j) = g(ivar) end if end do end do ! ! Interpolate the pressures at points (even, odd) and (odd, even). ! do i = 2, mx-1, 2 do j = 1, my, 2 press(i,j) = 0.5D+00 * (press(i-1,j)+press(i+1,j)) end do end do do j = 2, my-1, 2 do i = 1, mx, 2 press(i,j) = 0.5D+00 * (press(i,j-1)+press(i,j+1)) end do end do ! ! Interpolate the pressures at points (even,even). ! do j = 2, my-1, 2 do i = 2, mx-1, 2 press(i,j) = 0.5D+00 * (press(i-1,j-1)+press(i+1,j+1)) end do end do return end subroutine qbf ( xq, yq, it, in, bb, bx, by, nelemn, nnodes, node, np, xc, yc ) !*****************************************************************************80 ! !! QBF evaluates the quadratic basis functions. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer nelemn integer nnodes integer np real ( kind = 8 ) bb real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) c real ( kind = 8 ) d integer i1 integer i2 integer i3 integer in integer in1 integer in2 integer in3 integer inn integer it integer j1 integer j2 integer j3 integer node(nelemn,nnodes) real ( kind = 8 ) s real ( kind = 8 ) t real ( kind = 8 ) xc(np) real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) yq if ( in <= 3 ) then in1 = in in2 = mod(in,3)+1 in3 = mod(in+1,3)+1 i1 = node(it,in1) i2 = node(it,in2) i3 = node(it,in3) d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) t = 1.0D+00 +((yc(i2)-yc(i3))*(xq-xc(i1))+(xc(i3)-xc(i2))*(yq-yc(i1)))/d bb = t*(2.0D+00*t-1.0D+00) bx = (yc(i2)-yc(i3))*(4.0D+00*t-1.0D+00)/d by = (xc(i3)-xc(i2))*(4.0D+00*t-1.0D+00)/d else inn = in-3 in1 = inn in2 = mod(inn,3)+1 in3 = mod(inn+1,3)+1 i1 = node(it,in1) i2 = node(it,in2) i3 = node(it,in3) j1 = i2 j2 = i3 j3 = i1 d = (xc(i2)-xc(i1))*(yc(i3)-yc(i1))-(xc(i3)-xc(i1))*(yc(i2)-yc(i1)) c = (xc(j2)-xc(j1))*(yc(j3)-yc(j1))-(xc(j3)-xc(j1))*(yc(j2)-yc(j1)) t = 1.0D+00+((yc(i2)-yc(i3))*(xq-xc(i1))+(xc(i3)-xc(i2))*(yq-yc(i1)))/d s = 1.0D+00+((yc(j2)-yc(j3))*(xq-xc(j1))+(xc(j3)-xc(j2))*(yq-yc(j1)))/c bb = 4.0D+00 * s*t bx = 4.0D+00 * (t*(yc(j2)-yc(j3))/c+s*(yc(i2)-yc(i3))/d) by = 4.0D+00 * (t*(xc(j3)-xc(j2))/c+s*(xc(i3)-xc(i2))/d) end if return end function refbsp ( xq, yq, iq ) !*****************************************************************************80 ! !! REFBSP evaluates the linear basis functions in a reference triangle. ! ! Discussion: ! ! The reference triangle given here is not the best. The nodes ! are not given in the usual positions, and are listed in ! clockwise order instead of counterclockwise order! ! ! Diagram: ! ! 2 ! /| ! / | ! / | ! 1---3 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 June 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) XQ, YQ, the coordinates of a point in ! the reference triangle. ! ! Input, integer IQ, the index of a basis function ! in the reference triangle. ! ! Output, real ( kind = 8 ) REFBSP, the value of the IQ-th basis function ! at the point (XQ,YQ) in the reference triangle. ! implicit none integer iq real ( kind = 8 ) refbsp real ( kind = 8 ) xq real ( kind = 8 ) yq if ( iq == 1 ) then refbsp = 1.0D+00 - xq else if ( iq == 2 ) then refbsp = yq else if ( iq == 3 ) then refbsp = xq - yq end if return end subroutine refqbf ( x, y, in, bb, bx, by, etax, etay, xix, xiy ) !*****************************************************************************80 ! !! REFQBF evaluates the quadratic basis functions on reference triangle ! ! Diagram: ! ! 3 ! |\ ! 6 5 ! | \ ! 1-4-2 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 June 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, double precision X, Y, the coordinates of a point in ! the reference triangle. ! ! Input, integer IN, the index of a basis function ! in the reference triangle. ! ! Output, double precision BB, BX, BY, the value of the basis function ! and its X and Y derivatives at the point (X,Y) in the reference triangle. ! ! Input, real ( kind = 8 ) ETAX, ETAY, XIX, XIY, the values of ! dETA/dX, dETA/dY, dXI/dX and dXI/dY at (X,Y). ! implicit none real ( kind = 8 ) bb real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) etax real ( kind = 8 ) etay integer in real ( kind = 8 ) tbx real ( kind = 8 ) tby real ( kind = 8 ) x real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) y if ( in == 1) then bb = 1.0D+00-3.0D+00*x+2.0D+00*x*x tbx = -3.0D+00+4.0D+00*x tby = 0.0D+00 else if ( in == 2) then bb = -y+2.0D+00*y*y tbx = 0.0D+00 tby = -1.0D+00+4.0D+00*y else if (in == 3) then bb = -x+2.0D+00*x*x+y-4.0D+00*x*y+2*y*y tbx = -1.0D+00+4.0D+00*x-4.0D+00*y tby = 1.0D+00-4.0D+00*x+4.0D+00*y else if ( in == 4) then bb = 4.0D+00*y-4.0D+00*x*y tbx = -4.0D+00*y tby = 4.0D+00-4.0D+00*x else if ( in == 5) then bb = 4.0D+00*x*y-4.0D+00*y*y tbx = 4.0D+00*y tby = 4.0D+00*x-8.0D+00*y else if ( in == 6) then bb = 4.0D+00*x-4.0D+00*x*x-4.0D+00*y+4.0D+00*x*y tbx = 4.0D+00-8.0D+00*x+4.0D+00*y tby = -4.0D+00+4.0D+00*x else write(*,*)'REFQBF - Illegal value of IN = ',in stop end if bx = tbx * xix + tby * etax by = tbx * xiy + tby * etay return end subroutine resid ( area, g, indx, insc, isotri, iwrite, nelemn, neqn, & nnodes, node, np, nquad, phi, psi, res, reynld, xc, xm, yc, ym ) !*****************************************************************************80 ! !! RESID computes the residual. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) AREA(NELEMN), the area of the elements. ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer nelemn integer neqn integer nnodes integer np integer nquad real ( kind = 8 ) aij real ( kind = 8 ) ar real ( kind = 8 ) area(nelemn) real ( kind = 8 ) bb real ( kind = 8 ) bbb real ( kind = 8 ) bbbl real ( kind = 8 ) bbl real ( kind = 8 ) bbx real ( kind = 8 ) bby real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) det real ( kind = 8 ) etax real ( kind = 8 ) etay real ( kind = 8 ) g(neqn) integer i integer ibad integer ihor integer imax integer indx(np,2) integer insc(np) integer ip integer ipp integer iprs integer iq integer iqq integer iquad integer isotri(nelemn) integer it integer itype integer iver integer iwrite integer j integer jp integer ju integer jv integer node(nelemn,nnodes) real ( kind = 8 ) phi(nelemn,nquad,nnodes,3) real ( kind = 8 ) psi(nelemn,nquad,nnodes) real ( kind = 8 ) res(neqn) real ( kind = 8 ) reynld real ( kind = 8 ) rmax real ( kind = 8 ) test real ( kind = 8 ) ubc real ( kind = 8 ) ubdry real ( kind = 8 ) ubump real ( kind = 8 ) un(2) real ( kind = 8 ) unx(2) real ( kind = 8 ) uny(2) real ( kind = 8 ) visc real ( kind = 8 ) xc(np) real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) xm(nelemn,nquad) real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) ym(nelemn,nquad) real ( kind = 8 ) yq itype = -1 visc = 1.0D+00 / reynld res(1:neqn) = 0.0D+00 do it = 1, nelemn ar = area(it) / 3.0D+00 do iquad = 1, nquad yq = ym(it,iquad) xq = xm(it,iquad) if ( isotri(it) == 1 ) then call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) ar = det * area(it) / 3.0D+00 end if call uval(etax,etay,g,indx,isotri,it,nelemn,neqn, & nnodes,node,np,un,uny,unx,xc,xix,xiy,xq,yc,yq) ! ! For each basis function: ! do iq = 1, nnodes ip = node(it,iq) bb = phi(it,iquad,iq,1) bx = phi(it,iquad,iq,2) by = phi(it,iquad,iq,3) bbl = psi(it,iquad,iq) iprs = insc(ip) ihor = indx(ip,1) iver = indx(ip,2) if ( 0 < ihor ) then res(ihor) = res(ihor)-ar*bb*(un(1)*unx(1)+un(2)*uny(1)) end if if ( 0 < iver ) then res(iver) = res(iver)-ar*bb*(un(1)*unx(2)+un(2)*uny(2)) end if ! ! For another basis function, ! do iqq = 1, nnodes ipp = node(it,iqq) bbb = phi(it,iquad,iqq,1) bbx = phi(it,iquad,iqq,2) bby = phi(it,iquad,iqq,3) bbbl = psi(it,iquad,iqq) ju = indx(ipp,1) jv = indx(ipp,2) jp = insc(ipp) ! ! Horizontal velocity variable ! if ( 0 < ju ) then if ( 0 < ihor ) then res(ihor) = res(ihor)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2)))*g(ju) end if if ( 0 < iver ) then res(iver) = res(iver)+ar*bb*bbb*unx(2)*g(ju) end if if ( 0 < iprs ) then res(iprs) = res(iprs)+ar*bbx*bbl*g(ju) end if else if ( ju == itype ) then if ( ju == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,1,nelemn, & neqn,nnodes,node,np,xc,yc) else if ( ju == -1 ) then ubc = ubdry(1,yc(ipp)) end if if ( 0 < ihor ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*unx(1)+bbx*un(1)+bby*un(2))) res(ihor) = res(ihor)+ubc*aij end if if ( 0 < iver ) then aij = ar*bb*bbb*unx(2) res(iver) = res(iver)+ubc*aij end if if ( 0 < iprs ) then aij = ar*bbx*bbl res(iprs) = res(iprs)+ubc*aij end if end if ! ! Vertical velocity variable ! if ( 0 < jv ) then if ( 0 < ihor ) then res(ihor) = res(ihor)+ar*bb*bbb*uny(1)*g(jv) end if if ( 0 < iver ) then res(iver) = res(iver)+ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1)))*g(jv) end if if ( 0 < iprs ) then res(iprs) = res(iprs)+ar*bby*bbl*g(jv) end if else if ( jv == itype ) then if ( jv == -2 ) then ubc = ubump(g,indx,ipp,iqq,isotri,it,2,nelemn, & neqn,nnodes,node,np,xc,yc) else if ( jv == -1 ) then ubc = ubdry(2,yc(ipp)) end if if ( 0 < ihor ) then aij = ar*bb*bbb*uny(1) res(ihor) = res(ihor)+ubc*aij end if if ( 0 < iver ) then aij = ar*(visc*(by*bby+bx*bbx) & +bb*(bbb*uny(2)+bby*un(2)+bbx*un(1))) res(iver) = res(iver)+ubc*aij end if if ( 0 < iprs ) then aij = ar*bby*bbl res(iprs) = res(iprs)+ubc*aij end if end if ! ! Pressure variable ! if ( 0 < jp ) then if ( 0 < ihor ) then res(ihor) = res(ihor)-ar*bx*bbbl*g(jp) end if if ( 0 < iver ) then res(iver) = res(iver)-ar*by*bbbl*g(jp) end if end if end do end do end do end do ! ! The last equation is "reset" to require that the last pressure ! be zero. ! res(neqn) = g(neqn) rmax = 0.0D+00 imax = 0 ibad = 0 do i = 1, neqn test = abs(res(i)) if ( rmax < test ) then rmax = test imax = i end if if ( 1.0D-03 < test ) then ibad = ibad+1 end if end do if ( 1 <= iwrite ) then write(*,*)' ' write(*,*)'RESIDUAL INFORMATION:' write(*,*)' ' write(*,*)'Worst residual is number ',IMAX write(*,*)'of magnitude ',RMAX write(*,*)' ' write(*,*)'Number of "bad" residuals is ',IBAD,' out of ',NEQN write(*,*)' ' end if if ( 2 <= iwrite ) then write(*,*)'Raw residuals:' write(*,*)' ' i = 0 do j = 1, np if ( 0 < indx(j,1) ) then i = i+1 if ( abs(res(i)) <= 1.0D-03 ) then write(*,'(1x,a1,2i5,g14.6)')'U',i,j,res(i) else write(*,'(a1,a1,2i5,g14.6)')'*','U',i,j,res(i) end if end if if ( 0 < indx(j,2) ) then i = i+1 if ( abs(res(i)) <= 1.0D-03 ) then write(*,'(1x,a1,2i5,g14.6)')'V',i,j,res(i) else write(*,'(a1,a1,2i5,g14.6)')'*','V',i,j,res(i) end if end if if ( 0 < insc(j) ) then i = i+1 if ( abs(res(i)) <= 1.0D-03 ) then write(*,'(1x,a1,2i5,g14.6)')'P',i,j,res(i) else write(*,'(a1,a1,2i5,g14.6)')'*','P',i,j,res(i) end if end if end do end if return end subroutine setban ( indx, insc, maxrow, nband, nelemn, nlband, nnodes, & node, np, nrow ) !*****************************************************************************80 ! !! SETBAN computes the half band width. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer nelemn integer nnodes integer np integer i integer indx(np,2) integer insc(np) integer ip integer ipp integer iq integer iqq integer it integer iuk integer iukk integer j integer maxrow integer nband integer nlband integer node(nelemn,nnodes) integer nrow nlband = 0 do it = 1, nelemn do iq = 1, nnodes ip = node(it,iq) do iuk = 1, 3 if (iuk == 3) then i = insc(ip) else i = indx(ip,iuk) end if if ( 0 < i ) then do iqq = 1, nnodes ipp = node(it,iqq) do iukk = 1, 3 if (iukk == 3) then j = insc(ipp) else j = indx(ipp,iukk) end if if ( 0 < j ) then nlband = max(nlband,j-i) end if end do end do end if end do end do end do nband = nlband+nlband+1 nrow = nlband+nlband+nlband+1 write(*,*)' ' write(*,*)'SETBAN:' write(*,*)' ' write(*,*)' Lower bandwidth = ',nlband write(*,*)' Total bandwidth = ',nband write(*,*)' Required matrix rows = ',nrow if ( maxrow < nrow ) then write(*,*)'SETBAN - NROW is too large!' write(*,*)'The maximum allowed is ',maxrow stop end if return end subroutine setbas ( isotri, nelemn, nnodes, node, np, nquad, phi, psi, & xc, xm, yc, ym ) !*****************************************************************************80 ! !! SETBAS evaluates the basis functions at each integration point. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer nelemn integer nnodes integer np integer nquad real ( kind = 8 ) bb real ( kind = 8 ) bsp real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) det real ( kind = 8 ) etax real ( kind = 8 ) etay integer iq integer isotri(nelemn) integer it integer j integer node(nelemn,nnodes) real ( kind = 8 ) phi(nelemn,nquad,nnodes,3) real ( kind = 8 ) psi(nelemn,nquad,nnodes) real ( kind = 8 ) refbsp real ( kind = 8 ) xc(np) real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) xm(nelemn,nquad) real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) ym(nelemn,nquad) real ( kind = 8 ) yq do it = 1, nelemn do j = 1, nquad xq = xm(it,j) yq = ym(it,j) call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) do iq = 1, nnodes if ( isotri(it) == 0 ) then psi(it,j,iq) = bsp(it,iq,1,nelemn,nnodes,node,np,xc,xq,yc,yq) call qbf(xq,yq,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) else call refqbf(xq,yq,iq,bb,bx,by,etax,etay,xix,xiy) psi(it,j,iq) = refbsp(xq,yq,iq) end if phi(it,j,iq,1) = bb phi(it,j,iq,2) = bx phi(it,j,iq,3) = by end do end do end do return end subroutine setgrd ( ibump, indx, insc, isotri, iwrite, long, maxeqn, mx, my, & nelemn, neqn, nnodes, node, np, nx, ny, xbleft, xbrite, xlngth ) !*****************************************************************************80 ! !! SETGRD sets up the geometric grid. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer MX, the number of nodes in the X direction. ! ! Input, integer MY, the number of nodes in the Y direction. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer my integer nelemn integer nnodes integer np integer i integer ibump integer ic integer icnt integer ielemn integer indx(np,2) integer insc(np) integer ip integer ip1 integer ip2 integer isotri(nelemn) integer it integer iwrite integer jc integer jcnt logical long integer maxeqn integer mx integer nbleft integer nbrite integer neqn integer node(nelemn,nnodes) integer nx integer ny real ( kind = 8 ) xbleft real ( kind = 8 ) xbrite real ( kind = 8 ) xlngth write(*,*)' ' write(*,*)'SETGRD:' write(*,*)' ' ! ! Determine whether region is long or skinny. This will determine ! how we number the nodes and elements. ! if ( ny < nx ) then long = .true. write(*,*)'Using vertical ordering.' else long = .false. write(*,*)'Using horizontal ordering.' end if ! ! Report on isoparametric strategy: ! if ( ibump == 0 ) then write(*,*)'No isoparametric elements will be used.' else if ( ibump == 1 ) then write(*,*)'Isoparametric elements directly on bump.' else if ( ibump == 2 ) then write(*,*)'All elements above bump are isoparametric.' else if ( ibump == 3 ) then write(*,*)'All elements are isoparametric.' else write(*,*)'Unexpected value of IBUMP = ',ibump stop end if ! ! Compute node locations of bump corners based on X coordinates ! nbleft = nint(xbleft*(mx-1)/xlngth)+1 nbrite = nint(xbrite*(mx-1)/xlngth)+1 write(*,*)'Bump extends from ',xbleft,' at node ',nbleft write(*,*)' to ',xbrite,' at node ',nbrite ! ! Assign nodes to elements. ! neqn = 0 ielemn = 0 do ip = 1, np if ( long ) then ic = ((ip-1)/my)+1 jc = mod((ip-1),my)+1 else ic = mod((ip-1),mx)+1 jc = ((ip-1)/mx)+1 end if icnt = mod(ic,2) jcnt = mod(jc,2) ! ! If both the row count and the column count are odd, ! and we're not in the last row or top column, ! then we can define two new triangular elements based at the node. ! ! For horizontal ordering, ! given the following arrangement of nodes, for instance: ! ! 21 22 23 24 25 ! 16 17 18 19 20 ! 11 12 13 14 15 ! 06 07 08 09 10 ! 01 02 03 04 05 ! ! when we arrive at node 13, we will define ! ! element 7: (13, 23, 25, 18, 24, 19) ! element 8: (13, 25, 15, 19, 20, 14) ! ! ! For vertical ordering, ! given the following arrangement of nodes, for instance: ! ! 05 10 15 20 25 ! 04 09 14 19 24 ! 03 08 13 18 23 ! 02 07 12 17 22 ! 01 06 11 16 21 ! ! when we arrive at node 13, we will define ! ! element 7: (13, 25, 23, 19, 24, 18) ! element 8: (13, 15, 25, 14, 20, 19) ! if ( (icnt == 1.and.jcnt == 1).and.(ic /= mx).and.(jc /= my) ) then if ( long ) then ip1 = ip+my ip2 = ip+my+my ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip+2 node(ielemn,3) = ip2+2 node(ielemn,4) = ip+1 node(ielemn,5) = ip1+2 node(ielemn,6) = ip1+1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then isotri(ielemn) = 0 else if ( ibump == 2 ) then if ( nbleft <= ic .and. ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip2+2 node(ielemn,3) = ip2 node(ielemn,4) = ip1+1 node(ielemn,5) = ip2+1 node(ielemn,6) = ip1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then if ( jc == 1 .and. nbleft <= ic .and. ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else if ( ibump == 2 ) then if ( nbleft <= ic .and. ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if else ip1 = ip+mx ip2 = ip+mx+mx ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip2 node(ielemn,3) = ip2+2 node(ielemn,4) = ip1 node(ielemn,5) = ip2+1 node(ielemn,6) = ip1+1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then isotri(ielemn) = 0 else if ( ibump == 2 ) then if ( nbleft <= ic .and. ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if ielemn = ielemn+1 node(ielemn,1) = ip node(ielemn,2) = ip2+2 node(ielemn,3) = ip+2 node(ielemn,4) = ip1+1 node(ielemn,5) = ip1+2 node(ielemn,6) = ip+1 if ( ibump == 0 ) then isotri(ielemn) = 0 else if ( ibump == 1 ) then if ( jc == 1 .and. nbleft <= ic .and. ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else if ( ibump == 2 ) then if ( nbleft <= ic .and. ic < nbrite ) then isotri(ielemn) = 1 else isotri(ielemn) = 0 end if else isotri(ielemn) = 1 end if end if end if ! ! Left hand boundary, horizontal and vertical velocities specified. ! if ( ic == 1.and.1 < jc .and. jc < my ) then indx(ip,1) = -1 indx(ip,2) = -1 ! ! Right hand boundary, horizontal velocities unknown, vertical ! velocities specified ! else if ( ic == mx.and.1 < jc .and. jc < my) then neqn = neqn+1 indx(ip,1) = neqn indx(ip,2) = 0 ! ! Lower boundary, with isoperimetric triangle ! else if ( jc == 1 .and. isotri(ielemn) == 1 ) then indx(ip,1) = -2 indx(ip,2) = -2 ! ! Otherwise, just a wall ! else if ( ic == 1 .or. ic == mx .or. jc == 1 .or. jc == my ) then indx(ip,1) = 0 indx(ip,2) = 0 ! ! Otherwise, a normal interior node where both velocities are unknown ! else neqn = neqn+2 indx(ip,1) = neqn-1 indx(ip,2) = neqn end if if ( jcnt == 1 .and. icnt == 1 ) then neqn = neqn+1 insc(ip) = neqn else insc(ip) = 0 end if end do if ( 1 <= iwrite ) then write(*,*)' ' write(*,*)' I INDX 1, INDX 2, INSC' write(*,*)' ' do i = 1, np write(*,'(4i5)')i,indx(i,1),indx(i,2),insc(i) end do write(*,*)' ' write(*,*)'Isoparametric triangles:' write(*,*)' ' do i = 1, nelemn if ( isotri(i) == 1)write(*,*)i end do write(*,*)' ' write(*,*)' IT NODE(IT,*)' write(*,*)' ' do it = 1, nelemn write(*,'(7i6)') it,(node(it,i),i = 1,6) end do end if write(*,*)' ' write(*,*)'SETGRD: Number of unknowns = ',neqn if ( maxeqn < neqn ) then write(*,*)'SETGRD - Too many unknowns!' write(*,*)'The maximum allowed is MAXEQN = ',maxeqn write(*,*)'This problem requires NEQN = ',neqn stop end if return end subroutine setlin ( iline, indx, iwrite, long, mx, my, np, nx, ny, & xlngth, xprof ) !*****************************************************************************80 ! !! SETLIN determines unknown numbers along the profile line. ! ! Discussion: ! ! For our problem, the profile line has the equation X = XPROF. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer MX, the number of nodes in the X direction. ! ! Input, integer MY, the number of nodes in the Y direction. ! implicit none integer my integer np integer i integer iline(my) integer indx(np,2) integer ip integer itemp integer iwrite logical long integer mx integer nodex0 integer nx integer ny real ( kind = 8 ) xlngth real ( kind = 8 ) xprof ! ! Determine the number of a node on the profile line. ! itemp = nint ( ( 2.0D+00 * real ( nx - 1, kind = 8 ) * xprof ) / xlngth ) if ( long ) then nodex0 = itemp * ( 2 * ny - 1 ) + 1 else nodex0 = itemp + 1 end if write(*,*)' ' write(*,*)'SETLIN:' write(*,*)' ' write(*,*)' Profile generated at X = ',xprof write(*,*)' which is above node = ',nodex0 do i = 1, my if ( long ) then ip = nodex0+(i-1) else ip = nodex0+mx*(i-1) end if iline(i) = indx(ip,1) end do if ( 1 <= iwrite ) then write(*,*)' ' write(*,*)' Indices of unknowns along the profile line:' write(*,*)' ' write(*,'(5i5)') iline(1:my) end if return end subroutine setqud ( area, isotri, iwrite, nelemn, nnodes, node, np, nquad, & xc, xm, yc, ym ) !*****************************************************************************80 ! !! SETQUD sets midpoint quadrature rule information. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) AREA(NELEMN), the area of the elements. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer nelemn integer nnodes integer np integer nquad real ( kind = 8 ) area(nelemn) integer i integer ip1 integer ip2 integer ip3 integer isotri(nelemn) integer it integer iwrite integer j integer node(nelemn,nnodes) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) xc(np) real ( kind = 8 ) xm(nelemn,nquad) real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) yc(np) real ( kind = 8 ) ym(nelemn,nquad) do it = 1, nelemn ip1 = node(it,1) ip2 = node(it,2) ip3 = node(it,3) x1 = xc(ip1) x2 = xc(ip2) x3 = xc(ip3) y1 = yc(ip1) y2 = yc(ip2) y3 = yc(ip3) if ( isotri(it) == 0 ) then xm(it,1) = 0.5D+00*(x1+x2) xm(it,2) = 0.5D+00*(x2+x3) xm(it,3) = 0.5D+00*(x3+x1) ym(it,1) = 0.5D+00*(y1+y2) ym(it,2) = 0.5D+00*(y2+y3) ym(it,3) = 0.5D+00*(y3+y1) area(it) = 0.5D+00*abs((y1+y2)*(x2-x1)+(y2+y3)*(x3-x2)+(y3+y1)*(x1-x3)) else xm(it,1) = 0.5D+00 ym(it,1) = 0.5D+00 xm(it,2) = 1.0D+00 ym(it,2) = 0.5D+00 xm(it,3) = 0.5D+00 ym(it,3) = 0.0D+00 area(it) = 0.5D+00 end if end do if ( 3 <= iwrite ) then write(*,*)' ' write(*,*)'SETQUD: Element Areas and Quadrature points:' write(*,*)' ' do i = 1, nelemn write(*,*)i,area(i) do j = 1, nquad write(*,*)i,j,xm(i,j),ym(i,j) end do end do end if return end subroutine setxy ( iwrite, long, mx, my, np, nx, ny, xc, xlngth, yc, & ylngth, ypert ) !*****************************************************************************80 ! !! SETXY sets the grid coordinates based on the value of the parameter. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MX, the number of nodes in the X direction. ! ! Input, integer MY, the number of nodes in the Y direction. ! implicit none integer my integer np integer i integer ic integer ip integer iwrite integer jc logical long integer mx integer nx integer ny real ( kind = 8 ) xc(np) real ( kind = 8 ) xlngth real ( kind = 8 ) ybot real ( kind = 8 ) yc(np) real ( kind = 8 ) ylngth real ( kind = 8 ) ylo real ( kind = 8 ) ypert do ip = 1, np if ( long ) then ic = ( ( ip - 1 ) / my ) + 1 jc = mod ( ( ip - 1 ), my ) + 1 else ic = mod ( ( ip - 1 ), mx ) + 1 jc = ( ( ip - 1 ) / mx ) + 1 end if xc(ip) = (ic-1) * xlngth / (2*nx-2) ybot = -ypert * ( xc(ip) - 3.0D+00 ) * ( xc(ip) - 1.0D+00 ) ylo = max ( 0.0D+00, ybot ) yc(ip) = ( (my-jc)*ylo + (jc-1)*ylngth ) / (2*ny-2) end do if ( 2 <= iwrite ) then write(*,*)' ' write(*,*)'SETGRD:' write(*,*)' ' write(*,*)' I XC YC' write(*,*)' ' do i = 1, np write(*,'(i5,2f12.5)') i, xc(i), yc(i) end do end if return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) 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, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine trans ( det, etax, etay, it, nelemn, nnodes, node, np, xc, & xix, xiy, xq, yc, yq ) !*****************************************************************************80 ! !! TRANS calculates the element transformation mapping. ! ! Discussion: ! ! The element transformation mapping maps the reference element ! to a particular isoparametric element. The routine also outputs ! the determinant and partial derivatives of the transformation. ! ! Diagram of the mapping from reference to isoparametric element: ! ! 2 2 ! E /| / \ ! t 4 5 ===== > Y 4 5 ! a / | / \ ! 1-6-3 1--6----3 ! ! Xi X ! ! The form of the quadratic mapping is: ! ! x = a1 * xi^2 + b1 * xi * eta + c1 * eta^2 + d1 * xi + e1 * eta + f1 ! y = a2 * xi^2 + b2 * xi * eta + c2 * eta^2 + d2 * xi + e2 * eta + f2 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real ( kind = 8 ) DET, the determinant of the transformation, ! evaluated at (XQ,YQ). ! ! Output, real ( kind = 8 ) ETAX, ETAY, the Jacobian matrix entries ! dETA/dX and dETA/dY evaluated at (XQ,YQ). ! ! Input, integer IT, the index of the element ! containing (XQ,YQ). ! ! Input, integer NELEMN, the number of elements. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! ! Input, integer NODE(NELEMN,NNODES), contains the indices of ! the nodes that make up each element. ! ! Input, integer NP, the number of nodes. ! ! Input, real ( kind = 8 ) XC(NP), the X coordinates of nodes. ! ! Output, real ( kind = 8 ) XIX, XIY, the Jacobian matrix entries ! dXI/dX and dXI/dY evaluated at (XQ,YQ). ! ! Input, real ( kind = 8 ) XQ, the X coordinate of the point at ! which the mapping is to be evaluated. ! ! Input, real ( kind = 8 ) YC(NP), the Y coordinates of nodes. ! ! Input, real ( kind = 8 ) YQ, the Y coordinate of the point at ! which the mapping is to be evaluated. ! implicit none integer nelemn integer nnodes integer np real ( kind = 8 ) a1 real ( kind = 8 ) a2 real ( kind = 8 ) b1 real ( kind = 8 ) b2 real ( kind = 8 ) c1 real ( kind = 8 ) c2 real ( kind = 8 ) d1 real ( kind = 8 ) d2 real ( kind = 8 ) det real ( kind = 8 ) dxdeta real ( kind = 8 ) dxdxi real ( kind = 8 ) dydeta real ( kind = 8 ) dydxi real ( kind = 8 ) e1 real ( kind = 8 ) e2 real ( kind = 8 ) etax real ( kind = 8 ) etay integer i1 integer i2 integer i3 integer i4 integer i5 integer i6 integer it integer node(nelemn,nnodes) real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) x3 real ( kind = 8 ) x4 real ( kind = 8 ) x5 real ( kind = 8 ) x6 real ( kind = 8 ) xc(np) real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) xq real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) y3 real ( kind = 8 ) y4 real ( kind = 8 ) y5 real ( kind = 8 ) y6 real ( kind = 8 ) yc(np) real ( kind = 8 ) yq 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) ! ! Set the coefficients in the transformation: ! ! X = X(XI,ETA) ! Y = Y(XI,ETA) ! a1 = 2.0D+00*x3-4.0D+00*x6+2.0D+00*x1 b1 = -4.0D+00*x3-4.0D+00*x4+4.0D+00*x5+4.0D+00*x6 c1 = 2.0D+00*x2+2.0D+00*x3-4.0D+00*x5 d1 = -3.0D+00*x1-x3+4.0D+00*x6 e1 = -x2+x3+4.0D+00*x4-4.0D+00*x6 a2 = 2.0D+00*y3-4.0D+00*y6+2.0D+00*y1 b2 = -4.0D+00*y3-4.0D+00*y4+4.0D+00*y5+4.0D+00*y6 c2 = 2.0D+00*y2+2.0D+00*y3-4.0D+00*y5 d2 = -3.0D+00*y1-y3+4.0D+00*y6 e2 = -y2+y3+4.0D+00*y4-4.0D+00*y6 ! ! Compute partial derivatives d x/deta, d x/dxi, d y/deta/ dy/d xi, ! at point (xq,yq) in reference triangle. This is the Jacobian: ! ! ( dx/dxi dx/deta ) ! ( dy/dxi dy/deta ) ! dxdxi = 2.0D+00*a1*xq+b1*yq+d1 dxdeta = b1*xq+2.0D+00*c1*yq+e1 dydxi = 2.0D+00*a2*xq+b2*yq+d2 dydeta = b2*xq+2.0D+00*c2*yq+e2 ! ! Compute the determinant of the transformation. ! det = (2.0D+00*a1*b2-2.0D+00*a2*b1)*xq*xq & +(4.0D+00*a1*c2-4.0D+00*a2*c1)*xq*yq & +(2.0D+00*b1*c2-2.0D+00*b2*c1)*yq*yq & +(2.0D+00*a1*e2+b2*d1-b1*d2-2.0D+00*a2*e1)*xq & +(2.0D+00*c2*d1+b1*e2-b2*e1-2.0D+00*c1*d2)*yq+d1*e2-d2*e1 ! ! Compute the inverse jacobian ! ! ( dxi/dx dxi/dy ) ! ( deta/dx deta/dy ) ! xix = dydeta / det xiy = - dxdeta / det etax = - dydxi / det etay = dxdxi / det return end function ubdry ( iuk, yy ) !*****************************************************************************80 ! !! UBDRY sets the parabolic inflow in terms of the value of the parameter ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IUK, the index of the unknown. ! 1, the horizontal velocity. ! 2, the vertical velocity. ! ! Input, real ( kind = 8 ) YY, the Y coordinate of the boundary point. ! ! Output, real ( kind = 8 ) UBDRY, the value of the prescribed boundary ! flow component at this coordinate on the boundary. ! implicit none integer iuk real ( kind = 8 ) ubdry real ( kind = 8 ) yy if ( iuk == 1 ) then ubdry = ( -2.0D+00 * yy + 6.0D+00 ) * yy / 9.0D+00 else ubdry = 0.0D+00 end if return end function ubump ( g, indx, ip, iqq, isotri, it, iukk, nelemn, neqn, & nnodes, node, np, xc, yc ) !*****************************************************************************80 ! !! UBUMP calculates the sensitivity dU/dA on the bump. ! ! Discussion: ! ! This routine sets ! ! dU/dA = -uy * phi ! dV/dA = -vy * phi ! ! where PHI is the shape of the bump. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! implicit none integer nelemn integer neqn integer nnodes integer np real ( kind = 8 ) det real ( kind = 8 ) etax real ( kind = 8 ) etay real ( kind = 8 ) g(neqn) integer indx(np,2) integer ip integer iqq integer isotri(nelemn) integer it integer iukk integer node(nelemn,nnodes) real ( kind = 8 ) ubump real ( kind = 8 ) un(2) real ( kind = 8 ) unx(2) real ( kind = 8 ) uny(2) real ( kind = 8 ) xc(np) real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) yq if ( isotri(it) == 0 ) then xq = xc(ip) yq = yc(ip) else if ( iqq == 1 ) then xq = 0.0D+00 yq = 0.0D+00 else if ( iqq == 2 ) then xq = 1.0D+00 yq = 1.0D+00 else if ( iqq == 3 ) then xq = 1.0D+00 yq = 0.0D+00 else if ( iqq == 4 ) then xq = 0.5D+00 yq = 0.5D+00 else if ( iqq == 5 ) then xq = 1.0D+00 yq = 0.5D+00 else if ( iqq == 6 ) then xq = 0.5D+00 yq = 0.0D+00 end if call trans(det,etax,etay,it,nelemn,nnodes,node,np,xc,xix,xiy,xq,yc,yq) end if ! ! Calculate value of uy and vy (old solutions) at node ! call uval(etax,etay,g,indx,isotri,it,nelemn,neqn, & nnodes,node,np,un,uny,unx,xc,xix,xiy,xq,yc,yq) if ( iukk == 1 ) then ubump = uny(1) * (xc(ip)-1.0D+00) * (xc(ip)-3.0D+00) else if ( iukk == 2 ) then ubump = uny(2) * (xc(ip)-1.0D+00) * (xc(ip)-3.0D+00) else write(*,*)'UBUMP called for iukk = ',iukk stop end if return end subroutine uval ( etax, etay, g, indx, isotri, it, nelemn, neqn, nnodes, & node, np, un, uny, unx, xc, xix, xiy, xq, yc, yq ) !*****************************************************************************80 ! !! UVAL evaluates the velocities at a given quadrature point. ! ! Discussion: ! ! This routine also computes the spatial derivatives of the velocities. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer NNODES, the number of nodes per element, which ! is 6 for the quadratic triangular elements in use here. ! implicit none integer nelemn integer neqn integer nnodes integer np real ( kind = 8 ) bb real ( kind = 8 ) bx real ( kind = 8 ) by real ( kind = 8 ) etax real ( kind = 8 ) etay real ( kind = 8 ) g(neqn) ! integer i integer indx(np,2) integer ip integer iq integer isotri(nelemn) integer it integer iuk integer iun integer node(nelemn,nnodes) real ( kind = 8 ) ubc real ( kind = 8 ) ubdry real ( kind = 8 ) un(2) real ( kind = 8 ) unx(2) real ( kind = 8 ) uny(2) real ( kind = 8 ) xc(np) real ( kind = 8 ) xix real ( kind = 8 ) xiy real ( kind = 8 ) xq real ( kind = 8 ) yc(np) real ( kind = 8 ) yq un(1:2) = 0.0D+00 unx(1:2) = 0.0D+00 uny(1:2) = 0.0D+00 do iq = 1, nnodes if ( isotri(it) == 1 ) then call refqbf(xq,yq,iq,bb,bx,by,etax,etay,xix,xiy) else call qbf(xq,yq,it,iq,bb,bx,by,nelemn,nnodes,node,np,xc,yc) end if 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) else if ( iun == -1 ) then ubc = ubdry(iuk,yc(ip)) un(iuk) = un(iuk)+bb*ubc unx(iuk) = unx(iuk)+bx*ubc uny(iuk) = uny(iuk)+by*ubc end if end do end do return end subroutine uv_write ( f, indx, uv_unit, neqn, np, yc ) !*****************************************************************************80 ! !! UV_WRITE writes a velocity file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INDX(MAXNP,2), contains, for each node I, the ! index of the U and V velocities at that node, or 0. ! ! Input, integer UV_UNIT, the FORTRAN unit number associated ! with the output file. ! ! Input, integer NEQN, the number of equations. ! ! Input, integer NP, the number of nodes. ! ! Input, real ( kind = 8 ) YC(NP), the Y coordinates of nodes. ! implicit none integer neqn integer np real ( kind = 8 ) f(neqn) integer indx(np,2) integer ip integer k real ( kind = 8 ) u real ( kind = 8 ) ubdry integer uv_unit real ( kind = 8 ) v real ( kind = 8 ) yc(np) do ip = 1, np k = indx(ip,1) if ( k < 0 ) then u = ubdry ( 1, yc(ip) ) else if ( k == 0 ) then u = 0.0D+00 else u = f(k) end if k = indx(ip,2) if ( k < 0 ) then v = ubdry ( 2, yc(ip) ) else if ( k == 0 ) then v = 0.0D+00 else v = f(k) end if write ( uv_unit, '(2x,g14.6,2x,g14.6)' ) u, v end do return end subroutine xy_write ( xy_unit, np, xc, yc ) !*****************************************************************************80 ! !! XY_WRITE creates node coordinate data. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer XY_UNIT, the FORTRAN unit number associated ! with the output file. ! ! Input, integer NP, the number of nodes. ! ! Input, real ( kind = 8 ) XC(NP), YC(NP), the coordinates of nodes. ! implicit none integer np integer ip real ( kind = 8 ) xc(np) integer xy_unit real ( kind = 8 ) yc(np) do ip = 1, np write ( xy_unit, '(2x,g14.6,2x,g14.6)' ) xc(ip), yc(ip) end do return end