program femod1 ! !******************************************************************************* ! !! FEMOD1 is a one-dimensional interactive finite element program. ! ! ! Discussion: ! ! This version is fully functional but not very user friendly. ! In particular, you have to set up the problem by giving ! certain commands in the right order, but the code does not ! check that you do this, and doesn't tell you you need to do ! so. And the solution and output stages aren't easily ! controllable. ! ! This is an evolving version of the code. The original code ! simply asked you to define each item, according to an order ! defined internally. ! ! The new version, when complete, will allow you to change ! anything, everything, or nothing, in any order you like. ! It will, presumably, refuse to solve a system without enough ! information, and will demand that if some things are changed, ! others be changed as well. For instance, if you specified ! 10 nodes, and cubic elements, that's fine, because you need ! 3*n+1 nodes. But changing to quartic elements will require ! that you also specify a different number of nodes. ! ! I haven't gotten that all done yet, but I did manage to haul ! this code five or so years forward in time, plugging in the ! new COMRPN, and changing from the rigid procedure to a more ! friendly operation where the user is more in control... ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! implicit none ! integer, parameter :: maxelm = 100 integer, parameter :: maxfix = 80 integer, parameter :: maxgauss = 10 integer, parameter :: maxmat = 4 integer, parameter :: maxnpe = 6 integer, parameter :: maxpoi = 5 ! integer, parameter :: ldgk = 3 * ( maxnpe - 1 ) + 1 integer, parameter :: maxnod = maxelm * ( maxnpe - 1 ) + 1 integer, parameter :: maxrpn = 80 * ( 4 * maxmat + 2 ) ! character com character ( len = 40 ) command real dpsi(maxnpe) real ef(maxnpe) real ek(maxnpe,maxnpe) real gf(maxnod) real gk(ldgk,maxnod) real gklu(ldgk,maxnod) real gu(maxnod) integer ierror integer iexact integer ifrm character ( len = maxfix ) infix character ( len = maxrpn ) inform integer ipivot(maxnod) integer irhs integer irpn(maxrpn) integer kbc(2) integer kind(maxelm) integer kpoint(maxpoi) integer lenc logical s_eqi integer mat(maxelm) integer ml integer mu character ( len = 10 ) namvar integer nelem integer ngauss integer nmat integer nnode integer nodes(maxnpe,maxelm) integer npoint integer nrow real point(maxpoi) real psi(maxnpe) integer s_indexi character ( len = 80 ) title real value real vbc(2,2) real wgauss(maxgauss) real x(maxnod) real xgauss(maxgauss) ! call timestamp ( ) ! ! Initialize program data. ! call inifem ( iexact, inform, kbc, kind, kpoint, mat, maxelm, maxgauss, & maxnod, maxnpe, maxpoi, maxrpn, ml, mu, nelem, ngauss, nmat, nnode, & nodes, npoint, point, title, vbc, wgauss, x, xgauss ) ! ! Initialize the interactive compiler. ! com = 'I' ifrm = 1 namvar = ' ' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) ! ! Declare the name 'X' as a variable. ! com = 'A' namvar = 'X' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) ! ! Say hello. ! call hello ( maxelm, maxgauss, maxmat, maxnod, maxnpe, maxpoi ) ! ! Read the next user command. ! command = ' ' do if ( command(1:1) == '#' ) then command = ' ' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter command (HELP for menu)' end if read ( *, '(a)' ) command ! ! HELP = Print list of commands. ! if ( s_eqi ( command(1:1), 'H') ) then call help ! ! INIT = Initialize ! else if ( s_eqi ( command(1:1), 'I' ) ) then ! ! PREP = Preprocess (do the input in proper order). ! else if ( s_eqi ( command(1:4), 'PREP' ) ) then ierror = 0 call prep ( ierror, infix, inform, irpn, kbc, kind, kpoint, mat, maxelm, & maxfix, maxgauss, maxmat, maxnod, maxnpe, maxpoi, maxrpn, nelem, & ngauss, nmat, nnode, nodes, npoint, point, title, vbc, x ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' call bc_write ( kbc, vbc ) call element_write ( kind, mat, maxelm, maxnpe, nelem, nodes ) call gauss_write ( maxgauss, ngauss, wgauss, xgauss ) call node_write ( nnode, x ) call point_write ( kpoint, maxpoi, npoint, point ) call material_write ( inform, nmat ) ! ! QUIT = Quit. ! else if ( s_eqi ( command(1:1), 'Q' ) ) then write ( *, '(a)' ) 'Enter "Yes" to confirm you want to quit.' read ( *, '(a)' ) command if ( s_eqi ( command(1:1), 'Y' ) ) then exit end if ! ! SET = Set something. ! else if ( s_eqi ( command(1:3), 'SET' ) ) then if ( s_indexi ( command(5:), 'BC' ) > 0 ) then call bc_read ( kbc, vbc ) else if ( s_indexi ( command(5:), 'ELEM' ) > 0 ) then call element_read ( ierror, kind, mat, maxelm, maxnod, maxnpe, nelem, & nmat, nnode, nodes ) else if ( s_indexi ( command(5:), 'EXACT' ) /= 0 ) then call exact_read ( ierror, iexact, infix, inform, irpn, maxfix, maxmat, & maxrpn ) else if ( s_indexi ( command(5:), 'GAUSS' ) /= 0 ) then call gauss_read ( ierror, maxgauss, ngauss ) if ( ierror == 0 ) then call gauss_set ( ngauss, wgauss, xgauss ) end if else if ( s_indexi ( command(5:), 'MAT' ) /= 0 ) then call material_read ( ierror, infix, inform, irpn, maxfix, maxmat, & maxrpn, nmat ) else if ( s_indexi ( command(5:), 'NODE' ) /= 0 ) then call node_read ( kind, maxelm, maxnod, maxnpe, nelem, nnode, nodes, x ) else if ( s_indexi ( command(5:), 'POINT' )/=0 ) then call point_read ( kpoint, maxpoi, nnode, npoint, point ) else if ( s_indexi ( command(5:), 'TITLE' ) /= 0 ) then write ( *, '(a)' ) 'Enter a title for the problem' read ( *, '(a)' ) title end if ! ! PROC or SOLVE = Solve the system. ! else if ( s_eqi ( command(1:5), 'SOLVE' ) .or. & s_eqi ( command(1:4), 'PROC' ) ) then call proc ( dpsi, ef, ek, gf, gk, gklu, gu, ierror, infix, ipivot, irpn, & kbc, kind, kpoint, ldgk, mat, maxelm, maxfix, maxgauss, maxnod, & maxnpe, maxpoi, maxrpn, ml, mu, nelem, ngauss, nnode, nodes, & npoint, nrow, point, psi, vbc, wgauss, x, xgauss ) ! ! POST, print out stuff. ! else if ( s_eqi ( command(1:4), 'POST' ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' call solution_write ( dpsi, gu, iexact, infix, irpn, kind, maxelm, & maxfix, maxmat, maxnod, maxnpe, maxrpn, nelem, nodes, psi, title, x ) call exact_write ( iexact, inform, maxmat ) call error_write ( dpsi, gu, ierror, iexact, infix, irpn, kind, mat, & maxelm, maxfix, maxgauss, maxmat, maxnod, maxnpe, maxrpn, nelem, & ngauss, nodes, psi, wgauss, x, xgauss ) ! ! TYPE = Print out something. ! else if ( s_eqi ( command(1:4), 'TYPE' ) ) then if ( s_indexi ( command(5:), 'BC' ) > 0 ) then call bc_write ( kbc, vbc ) else if ( s_indexi ( command(5:), 'ELEM' ) > 0 ) then call element_write ( kind, mat, maxelm, maxnpe, nelem, nodes ) else if ( s_indexi ( command(5:), 'ERROR' ) /= 0 ) then call error_write ( dpsi, gu, ierror, iexact, infix, irpn, kind, mat, & maxelm, maxfix, maxgauss, maxmat, maxnod, maxnpe, maxrpn, nelem, & ngauss, nodes, psi, wgauss, x, xgauss ) else if ( s_indexi ( command(5:), 'EXACT' ) > 0 ) then call exact_write ( iexact, inform, maxmat ) else if ( s_indexi ( command(5:), 'GAUSS' ) > 0 ) then call gauss_write ( maxgauss, ngauss, wgauss, xgauss ) else if ( s_indexi ( command(5:), 'MAT' ) > 0 ) then call material_write ( inform, nmat ) else if ( s_indexi ( command(5:), 'NODE' ) > 0 ) then call node_write ( nnode, x ) else if ( s_indexi ( command(5:), 'POINT') > 0 ) then call point_write ( kpoint, maxpoi, npoint, point ) else if ( s_indexi ( command(5:), 'SOL' ) > 0 ) then call solution_write ( dpsi, gu, iexact, infix, irpn, kind, maxelm, & maxfix, maxmat, maxnod, maxnpe, maxrpn, nelem, nodes, psi, title, x ) else if ( s_indexi ( command(5:), 'SYS' ) > 0 ) then irhs = 1 call gkf_write ( gf, gk, ml, mu, nnode, nrow, irhs ) else if ( s_indexi ( command(5:), 'TITLE' ) > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' end if ! ! # Comment ! else if ( command(1:1) == '#' ) then else if ( len_trim ( command ) > 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD1 ignored the command:' write ( *, '(a)' ) ' "' // trim ( command ) //'"' end if end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD1:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine bc_apply ( gf, gk, ierror, kbc, maxnod, ml, mu, nnode, nrow, vbc ) ! !******************************************************************************* ! !! BC_APPLY modifies the global arrays to account for boundary conditions. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input/output, real GF(MAXNOD), the global force vector. ! ! Input/output, real GK(NROW,NNODE), the global stiffness matrix. ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, there was an illegal boundary condition code, or other problem. ! ! Input, integer KBC(2). KBC(1) and KBC(2) are the boundary condition ! codes at the left and right endpoints of the region. ! 1: U = value1; ! 2: K dU/dx = value1 * U - value2 ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NROW, the leading dimension of the global stiffness matrix. ! ! Input, real VBC(2,2), contains constants for boundary conditions. ! VBC(1,1) and VBC(2,1) are value1 and value2 for the left endpoint; ! VBC(1,2) and VBC(2,2) are value1 and value2 for the right endpoint. ! implicit none ! integer maxnod integer nnode integer nrow ! real gf(maxnod) real gk(nrow,nnode) integer i integer ierror integer ihi integer ilo integer irow integer kbc(2) integer ml integer mu real vbc(2,2) ! ! Apply the boundary condition at node 1. ! ! If the boundary condition is ! U(1) = VBC(1,1), ! subtract GK(I,1) * U(1) from every equation, ! and replace the first equation by the boundary condition. ! if ( kbc(1) == 1 ) then ihi = min ( 1+ml, nnode ) do i = 1, ihi irow = i-1+ml+mu+1 gf(i) = gf(i) - gk(irow,1) * vbc(1,1) gk(irow,1) = 0.0E+00 irow = 1-i+ml+mu+1 gk(irow,i) = 0.0E+00 end do irow = 1-1+ml+mu+1 gk(irow,1) = 1.0E+00 gf(1) = vbc(1,1) else if ( kbc(1) == 2 ) then irow = 1-1+ml+mu+1 gk(irow,1) = gk(irow,1) + vbc(1,1) gf(1) = gf(1) + vbc(2,1) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BC_APPLY - Fatal error!' write ( *, '(a,i6)' ) ' Illegal boundary condition type = ', kbc(1) return end if ! ! Apply the boundary condition at node NNODE. ! if ( kbc(2) == 1 ) then ilo = max(nnode-ml,1) do i = ilo, nnode irow = i-nnode+ml+mu+1 gf(i) = gf(i)-gk(irow,nnode)*vbc(1,2) gk(irow,nnode) = 0.0E+00 irow = nnode-i+ml+mu+1 gk(irow,i) = 0.0E+00 end do irow = nnode-nnode+ml+mu+1 gk(irow,nnode) = 1.0E+00 gf(nnode) = vbc(1,2) else if ( kbc(2)==2 ) then irow = nnode-nnode+ml+mu+1 gk(irow,nnode) = gk(irow,nnode) + vbc(1,2) gf(nnode) = gf(nnode) + vbc(2,2) else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BC_APPLY - Fatal error!' write ( *, '(a,i6)' ) ' Illegal boundary condition type = ', kbc(2) return end if return end subroutine bc_read ( kbc, vbc ) ! !******************************************************************************* ! !! BC_READ reads the boundary conditions from the user. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Output, integer KBC(2). KBC(1) and KBC(2) are the boundary condition ! codes at the left and right endpoints of the region. ! 1: U = value1; ! 2: K dU/dx = value1 * U - value2 ! ! Output, real VBC(2,2), contains constants for boundary conditions. ! VBC(1,1) and VBC(2,1) are value1 and value2 for the left endpoint; ! VBC(1,2) and VBC(2,2) are value1 and value2 for the right endpoint. ! implicit none ! integer i integer kbc(2) integer kval real val1 real val2 real vbc(2,2) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BC_READ:' write ( *, '(a)' ) ' Set up the boundary conditions.' do i = 1, 2 if ( i == 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Boundary conditions on the LEFT:' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Boundary conditions on the RIGHT:' end if do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Choose your boundary condition type:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' 1 for the condition " U = V1",' write ( *, '(a)' ) ' 2 for the condition "K dU/dx = V1 U + V2".' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter the boundary condition type:' read ( *, * ) kval if ( kval == 1 ) then write ( *, '(a)' ) 'Enter the value V1 in condition "U = V1"' read ( *, * ) val1 kbc(i) = kval vbc(1,i) = val1 vbc(2,i) = 0.0E+00 exit else if ( kval == 2 ) then write ( *, '(a)' ) 'Enter the values V1, V2 in BC "K dU/dx = V1 U + V2"' read ( *, * ) val1, val2 kbc(i) = kval vbc(1,i) = val1 vbc(2,i) = val2 exit else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BC_READ - Warning!' write ( *, '(a)' ) ' Your input was not acceptable:' write ( *, '(a,i6)' ) ' KVAL = ', kval write ( *, '(a)' ) ' but values of 1 or 2 were expected.' end if end do end do return end subroutine bc_write ( kbc, vbc ) ! !******************************************************************************* ! !! BC_WRITE prints the boundary conditions. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer KBC(2). KBC(1) and KBC(2) are the boundary condition ! codes at the left and right endpoints of the region. ! 1: U = value1; ! 2: K dU/dx = value1 * U - value2 ! ! Input, real VBC(2,2), contains constants for boundary conditions. ! VBC(1,1) and VBC(2,1) are value1 and value2 for the left endpoint; ! VBC(1,2) and VBC(2,2) are value1 and value2 for the right endpoint. ! implicit none ! integer i integer kbc(2) real vbc(2,2) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BC_WRITE:' write ( *, '(a)' ) ' Boundary conditions:' do i = 1, 2 write ( *, '(a)' ) ' ' if ( i == 1 ) then write ( *, '(a)' ) 'At the left boundary point,' else write ( *, '(a)' ) 'At the right boundary point,' end if write ( *, '(a,i6)' ) 'the boundary condition is of type ', kbc(i) if ( kbc(i) == 1 ) then write ( *, '(a,g14.6)' ) 'U = ', vbc(1,i) else write ( *, '(a,g14.6,a,g14.6)' ) 'K * dU/dx = ', vbc(1,i), & ' * U - ', vbc(2,i) end if end do return end subroutine ch_cap ( c ) ! !******************************************************************************* ! !! CH_CAP capitalizes a single character. ! ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none ! character c integer itemp ! itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) ! !******************************************************************************* ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! ! Examples: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! Modified: ! ! 14 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C1, C2, the characters to compare. ! ! Output, logical CH_EQI, the result of the comparison. ! implicit none ! logical ch_eqi character c1 character c2 character cc1 character cc2 ! cc1 = c1 cc2 = c2 call ch_cap ( cc1 ) call ch_cap ( cc2 ) if ( cc1 == cc2 ) then ch_eqi = .true. else ch_eqi = .false. end if return end function ch_is_alpha ( c ) ! !******************************************************************************* ! !! CH_IS_ALPHA returns TRUE if C is an alphabetic character. ! ! ! Modified: ! ! 05 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, a character to check. ! ! Output, logical CH_IS_ALPHA is TRUE if C is an alphabetic character. ! implicit none ! character c logical ch_is_alpha ! if ( ( lle ( 'a', c ) .and. lle ( c, 'z' ) ) .or. & ( lle ( 'A', c ) .and. lle ( c, 'Z' ) ) ) then ch_is_alpha = .true. else ch_is_alpha = .false. end if return end subroutine ch_to_digit ( c, digit ) ! !******************************************************************************* ! !! CH_TO_DIGIT returns the integer value of a base 10 digit. ! ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer DIGIT, the corresponding integer value. If C was ! 'illegal', then DIGIT is -1. ! implicit none ! character c integer digit ! if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c == ' ' ) then digit = 0 else digit = -1 end if return end subroutine comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) ! !******************************************************************************* ! !! COMRPN can translate formulas you type in and evaluate them. ! ! ! This means that you can design interactive FORTRAN programs ! which can input their formulas at run time, rather than ! recompiling every time you want a new formula. ! ! If you declare vectors or matrices, you may reference entries ! such as X(3) or A(3,2) or even X(2+3). ! ! ! Formulas are made from operators, constants, punctuation, ! variables, and functions. ! ! ! The list of legal operators includes: ! ! + - * / ** ^ = \ ! ! * means multiplication, and is standard matrix multiplication if ! both arguments are (conformable) matrices. If one argument is a ! scalar, it multiplies all the entries of the other argument. ! If both arguments are row or column vectors, the dot product ! is taken. ! ! Thus, ! ! x is a scalar equal to 2, ! y is a scalar equal to 3, ! u is a vector equal to (1, 2, 3), ! v is a vector equal to (1, 0, 2), ! ! then ! ! 2*3 results in 6. ! ! x*y results in 6. ! ! u*v results in 7. ! ! x*v would result in the vector (2, 0, 4). ! ! ! / means division, as in A/B, but for matrices the only allowed ! form requires that B be a scalar, in which case each element of ! A is divided by B. ! ! ! + and - are allowed for pairs of scalars, vectors or matrices of ! the same order, and also for a square matrix and a scalar, ! in which case A+2 is interpreted as A+2*I. ! ! ! ** or ^ means exponentiation. In the scalar case, any ! nonnegative scalar can be taken to any power, positive, zero, or ! negative. A negative scalar may only be taken to an integer ! power. ! ! ** is also supported for square matrices, which can be taken ! to any nonnegative integer power. Nonsingular matrices may ! also be taken to negative integer powers. ! ! ! = is used to assign values. For vectors or matrices, this may ! also be used to assign a single entry, as in A(3,2) = 7. ! ! ! / is standard scalar division. It is used in a formula like ! ! X/Y ! ! or ! ! (A*B)/(X+1) ! ! Normally, X would be a scalar quantity. However, as NONSTANDARD ! shorthand, you may 'divide' a matrix by a scalar, in which case ! each entry of the matrix is divided by the scalar, e.g. ! ! A/2 = (1/2)*A. ! ! MATALG will NOT allow you to use the "/" operator to represent ! multiplication by an inverse matrix. Thus, if A is a matrix, ! the statement ! ! B/A ! ! is illegal. However, the statement ! ! INV(A)*B ! ! will compute the inverse of A and multiply by B, and the ! statement ! ! A \ B ! ! will compute the LU factorization of A, and use that to ! compute the inverse of A times B. ! ! ! Constants: ! ! ! Real and integer constants may appear in your formulas, as well ! as the symbol 'PI' and the machine unit roundoff number 'EPS' ! which is the power of two such that 1+EPS>1 but 1+(EPS/2) = 1. ! ! ! Punctuation: ! ! ! Blanks may be used anywhere. ! ! Commas are used to separate arguments in a function, such as ! MAX(X,Y). ! ! Parentheses are used: ! ! to group quantities: (X+Y)*Z ! ! to reference an entry of a vector or matrix: A(2,2) ! ! to mark the argument of a function: SIN(X), MAX(X,Y) ! ! ! Functions and operators: ! ! ! S, S1, S2: Arguments which may only be scalar. ! V : Arguments may only be a scalar or vector. ! M : Arguments may only be scalar or square matrix. ! MV : Arguments may only be scalar, vector, or square matrix. ! * : Arguments which may be scalar, vector, or matrix. ! I : Arguments which may only be positive integers. ! ! ABS(*) Absolute value of *. ! ! ACOS(S) The arc cosine of S. ! -1 < = S <=1. ! ! ALOG(S) Natural logarithm of S. ! S must be positive. ! ! ALOG10(S) Logarithm base 10 of S. ! S must be positive. ! ! ASIN(S) The arc sine of S. ! -1 < = S <=1. ! ! ATAN(S) The arc tangent of S. ! ! ATAN2(S1,S2) Arc tangent of (S1/S2) ! Correctly computes ATAN2(1.0,0.0) = PI/2. ! ! COS(S) The cosine of S, with S measured in radians. ! ! COSD(S) The cosine of S, with S measured in degrees. ! ! COSH(S) Hyperbolic cosine of S. ! ! DET(M) Determinant of square matrix M. ! ! DIAG(*) The diagonal entries of *, stored in a vector. ! ! E The base of the natural logarithm system. ! E = 2.71828... ! You may not change the value of E. ! ! EPS The machine epsilon, or unit roundoff number. ! EPS is the power of 2 such that ! ! 1 < 1+EPS and 1 = 1+(EPS/2). ! ! You may not change the value of EPS. ! ! EVAL(M) Real and imaginary parts of eigenvalues of matrix ! M, stored in a matrix of N rows and 2 columns, ! with the real parts in column 1, and the ! imaginary parts in column 2. ! ! EVEC(M) Eigenvectors of square matrix M, stored in a ! square matrix of same size as M. ! ! (Eigenvectors corresponding to a complex pair: ! If eigenvalues j and j+1 are a complex pair, ! then the eigenvector for eigenvalue j is ! column j + i*column j+1, and the eigenvector ! for eigenvalue j+1 is column j - i*column j+1.) ! ! EXP(S) Exponential of S. ! ! GCF(I,J) Greatest common factor of two integers. ! ! HILBERT(I) The Hilbert matrix of order I. ! ! HILBINV(I) The inverse of the Hilbert matrix of order I. ! ! HOUSE(V) The Householder elementary reflector matrix H, ! defined as ! ! H = I - 2*(V * TRANSPOSE(V)) / (TRANS(V)*V) ! ! ID(I) The matrix identity of order I. ! ID(3) is the 3 by 3 identity, for instance. ! ! INT(*) Truncates real values to their integer part. ! INT(1.0) = INT(1.1) = INT(1.9) = 1. ! ! INV(M) The inverse matrix of the square matrix M. ! If M is singular, there will be no inverse. ! ! IVAL(M) Imaginary parts of eigenvalues of square matrix M ! stored in a column vector. ! ! KRYLOV(M,V,I) Returns a rectangular matrix whose columns are ! a sequence of I Krylov vectors: ! ! V, M*V, M*M*V, ..., M**(I-1)*V. ! ! LCM(I,J) Least common multiple of two integers. ! ! LN(S) The natural logarithm of S. ! S must be greater than 0. ! ! LOG(S) The natural logarithm of S. ! S must be greater than 0. ! ! LOG10(S) The logarithm base 10 of S. ! S must be greater than 0. ! ! LOG2(S) The logarithm base 2 of S. ! S must be greater than 0. ! ! MAX(S1,S2) The maximum of S1 and S2. ! ! MIN(S1,S2) Minimum of S1 and S2. ! ! NEG(*) Changes sign of *. ! ! NINT(*) Nearest integer value to entries of *. ! ! NORM0(*) Maximum or infinity norm of a vector or matrix. ! NORM0(S) = ABS(S) ! NORM0(V) returns the maximum of the absolute ! values of the entries of V. ! NORM0(M) returns the maximum of the sum of the ! absolute values of the entries in each row. ! ! NORM1(*) The L1-norm of a vector or matrix. ! NORM1(S) = ABS(S). ! NORM1(V) returns the sum of the absolute values ! of the entries of V. ! NORM0(M) returns the maximum of the sum of the ! absolute values of the entries in each column. ! ! NORM2(MV) L2-norm, Euclidean norm or root-mean-square norm ! of a vector or a square matrix M. NORM2(V) returns ! the square root of the sum of the squares of the ! entries of V. NORM2(M) returns the maximum magnitude ! of the eigenvalues of Transpose(M)*M. ! ! NORMF(*) The Frobenius norm. NORMF(M) returns the square ! root of the sum of the squares of all the entries ! of a matrix M. NORMF may also be applied to a ! vector, giving the same results as NORM2(V). ! ! NORMS(MV) The spectral "norm". NORMS(M) returns the ! value of the maximum of the absolute values of ! all the eigenvalues of the square matrix M. ! ! Note that the spectral norm is NOT a true ! norm, although it is frequently used as though ! it were. ! ! PI 3.14159265358979... ! You may not change the value of PI. ! ! POLY(V,M) Polynomial evaluation. V contains the ! coefficients of the polynomial, with V(1) the ! coefficient of M**(N-1) and V(N) the constant ! term. ! ! M is the scalar or square matrix argument of the ! polynomial. ! ! RAN(*) Fills * with random numbers between 0 and 1. ! ! RCOND(M) RCOND(M) is the estimate of the reciprocal ! condition number of the matrix M, as computed ! by the LINPACK routine SGECO. ! ! The condition number of a matrix M is defined as ! ! COND(M) = NORM(M) * NORM(INVERSE(M)) ! ! RCOND assumes that the L-1 norm is used to define ! COND(M), and then RCOND(M) is defined by ! ! RCOND(M) = 1 / COND(M). ! ! If the matrix M is singular, then the condition ! number is infinite, and RCOND will be zero. ! ! ROUND(*) Replaces all "small" entries of * by 0. ! ! RVAL(M) Real parts of eigenvalues of square matrix M, ! stored in a column vector. ! ! SIN(S) The sine of S, with S measured in radians. ! ! SIND(S) The sine of S, with S measured in degrees. ! ! SINH(S) Hyperbolic sine of S. ! ! SQRT(S) Square root of S. S must be non-negative. ! ! STEP(S) Heavyside step function. ! ! STEP = 0 if S < 0, ! STEP = 1 if 0 <= S. ! ! TAN(S) Tangent of S. ! ! TAND(S) Tangent of S measured in degrees. ! ! TANH(S) Hyperbolic tangent of S. ! ! TRACE(*) Sum of diagonal elements of *. ! ! TRANS(*) Transpose of matrix or vector. ! ! A = TRANS(A) is a legal formula if A is square. ! ! ZERO(I) Matrix zero of order I. ! ZERO(4) is the 4 by 4 zero matrix. ! ! I! Factorial of I. I should be an integer from 0 to ! 25. 0! = 1, 1!=1, 2!=1*2=2, 3!=1*2*3=6, and so on. ! For large I, the exact value is not returned. ! !******************************************************************************* ! ! To add a new function to COMRPN: ! ! Update the information in the IOPSYM, IPRSYM, and SYMBOL ! arrays. ! ! Insert text into FUNVAL or FUNSCL which evaluates the ! function, and determines its dimensions. ! !******************************************************************************* ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, CHARACTER COM, the command to be executed, which determines ! what other input must be supplied. ! ! ! COM = 'A' Add (=declare) a variable: ! ! Causes the program to add the variable named in NAMVAR ! to its list of symbols. From now on, it may be used in ! formulas, and values may be assigned to it. If you ! declare a name which already exists, then if that name ! was already declared by the program, your command is ! ignored. Otherwise, your new definition overrides your old one. ! ! Input for the "A" command includes NAMVAR, NROW, and NCOL. ! ! ! COM = 'D' delete a variable: ! ! This command can be used to free up memory. You ! must also give the name of the variable to be deleted. ! The variable must be one you have already defined. ! ! Input for the "D" command includes NAMVAR. ! ! ! COM = 'E' evaluate a formula: ! ! Using current values of variables, evaluate the formula. ! The formula number is given by IFRM. The value of the ! formula is returned in VALUE, and its dimensions in NROW and NCOL. ! ! Input to the "E" command includes IFRM. ! Output from the "E" command includes VALUE. ! ! ! COM = 'F', Formula is to be compiled: ! ! The user assigns an index, IFRM, to the formula, so that ! it can be referred to later. The formula is stored in ! INFIX on input. The formula may refer to variables which ! have not yet been declared. ! ! Output from the "F" command includes IFRM. ! ! ! COM = 'G', Formula is to be compiled, all variables ! in the formula have already been declared: ! ! The same as 'F' except that the formula may not refer to ! any variables which have not been declared. ! ! Output from the "G" command includes IFRM. ! ! ! COM = 'I' Initialize: ! ! This must be the first call, to initialize data. Also, ! if old data is to be flushed out, this command can be used. ! ! ! COM = 'L' List symbols and values: ! ! This command prints out the internal symbols, variables ! and values, and the formulas you have supplied. ! ! The value of NAMVAR determines what is printed. ! ! If NAMVAR = ' ', all user variables are printed. ! If NAMVAR = 'ALL', internal variables are also printed. ! If NAMVAR = 'DEBUG', all that and more is printed. ! ! Otherwise, NAMVAR is assumed to be the name of a ! particular variable and only that is printed. ! ! If NAMVAR = '#', the contents of VALUE are listed. ! ! ! COM = 'R' Read value of symbol: ! ! The program returns the current value of the variable ! NAMVAR. VALUE will contain the value, and NROW and NCOL ! the dimensions. ! ! ! COM = 'V' assign value to symbol: ! ! This command assigns the value of the numbers in VALUE ! to the variable NAMVAR, with NROW and NCOL specifying ! the dimension of the variable. ! ! Output, INTEGER IERROR, an error flag. ! ! 0 Means no error occurred. ! ! 1 means some kind of error. These errors are usually ! signaled by a printed message. They can include the ! following problems: ! ! On adding a variable, if the name was already in use, ! either by you or the program, the program returns ! with this error warning. ! ! On setting a variable, if the name supplied does ! not represent any recognized variable, the program ! ignores the command and returns. ! ! On compiling a formula, the program will return with ! an error warning if any of the following occur: ! ! The formula does not compile. An unknown variable, ! missing parentheses, garbled information or mistyping ! may be responsible. For that matter, the compiler ! may fail to compile perfectly good, but complicated ! formulas. ! ! There is not enough space in IRPN to store the RPN code. ! In this case, you must either re-initialize, or increase ! the storage available in IRPN and signal this with ! an increased value of MAXRPN. ! ! On evaluating a formula, an error can occur if there ! is no formula corresponding to the given index IFRM. ! ! Certain arithmetic errors, such as division ! by zero, may be caught and flagged by the compiler. ! ! Input/output, INTEGER IFRM. ! ! COMRPN can store more than one formula internally, at ! one time. COMRPN refers to the formulas by an index ! number. When a formula is entered, with the "F" or "G" ! command, it is assigned an index number, whose value is ! returned to the user. Then, whenever the formula is to ! be evaluated, the user must input the value of IFRM, so ! that COMRPN knows which formula to evaluate. In many ! cases, only one formula is being used, so that IFRM is ! usually just 1. ! ! For the "E" command, IFRM is the index of the formula to evaluate. ! ! For the "F" or "G" commands, the index assigned to the formula. ! ! Input, CHARACTER ( len = * ) INFIX, for the F or G commands, ! contains the formula to be evaluated. ! ! Input/output, INTEGER IRPN(MAXRPN), for the F, G or E commands, ! a vector used to store the compiled versions of the input formulas. ! ! Input, INTEGER MAXROW, the maximum number of rows and columns for matrices. ! If you aren't using matrices, set MAXROW to 1. ! ! Input, INTEGER MAXRPN, the maximum length of IRPN. 80 should be enough. ! ! Input, CHARACTER*(*) NAMVAR, used for the A, D, L, R and ! V commands, containing the name of the variable to be ! added, deleted, listed, read or valued. ! ! NAMVAR should not be a blank, nor should it be the name 'VOID'. ! ! For the L command only, using NAMVAR = 'DEBUG' produces ! extra output. ! ! Input/output, INTEGER NCOL, the number of columns of the variable ! specified by NAMVAR. ! ! Input/output, INTEGER NROW, the number of rows of the variable ! specified by NAMVAR. ! ! Workspace, CHARACTER*100 OUTPUT. ! ! Input/output, REAL VALUE(MAXROW,MAXROW), used to pass values in and out ! of the program. The actual size of the object in VALUE ! is (NROW,NCOL). Thus to pass a scalar, set VALUE(1,1) ! to the value and set NROW = NCOL=1. To pass a 2 by 3 ! matrix, set the values of VALUE and set NROW = 2, NCOL=3. ! ! When the program is asked to return a value in VALUE, it ! sets NROW and NCOL as well. ! ! For the V command, VALUE is the input value to be ! assigned to a variable. ! ! For the E command, VALUE is the output value of the formula. ! ! For the R command, VALUE is the output value of the variable. ! !******************************************************************************* ! ! Internal variables: ! ! Internal, INTEGER IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Internal, INTEGER IPVAL(MAXSYM), contains, for each ! symbol, the address in VALSYM where the value of the ! symbol is stored, if it is a scalar. If the symbol ! represents a vector or matrix, IPVAL(IARG) points to ! the location of the first entry. ! ! Internal, INTEGER MAXSYM, controls the maximum number ! of symbols allowed. This includes permanent symbols ! defined by the program, as well as symbols the user ! declares while running the program. ! implicit none ! integer, parameter :: maxsym = 150 integer, parameter :: maxval = 3000 integer, parameter :: maxlen = 20 ! integer maxrpn ! character com integer ierror integer, save :: ifinis integer, save :: ifree integer ifreesv integer ifrm integer ihi integer ilo integer imin integer implic integer, save :: indx1 integer, save :: indx2 integer, save :: ineg character ( len = * ) infix integer, dimension ( maxsym ), save :: intsym integer, dimension ( maxsym ), save :: iopsym integer, dimension ( maxsym ), save :: iprsym integer, dimension ( maxsym ), save :: ipval integer, save :: irad integer irpn(maxrpn) integer, dimension ( maxsym ), save :: istack integer maxfix integer maxfrm character ( len = * ) namvar character ( len = maxlen ) namvr integer ncol integer, save :: nints integer nrow integer nrpn integer, save :: nsym integer, save :: nsymp integer, save :: nsyms integer nuvar character ( len = 100 ) output logical s_eqi logical s_paren_check character ( len = maxlen ), dimension ( maxsym ), save :: symbol real value real, dimension ( maxval ), save :: valsym ! ierror = 0 implic = 0 namvr = namvar call s_blank_delete ( namvr ) ncol = 1 nrow = 1 maxfix = len ( infix ) if ( maxfix > 0 ) then maxfrm = ( maxrpn / maxfix ) else maxfrm = 0 end if if ( s_eqi ( com,'E') .or. s_eqi ( com, 'G' ) ) then if ( ifrm <= 0 .or. ifrm > maxfrm ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a,i6)' ) ' Illegal formula index = ', ifrm ierror = 1 return end if end if ! ! COM = I initialize. ! if ( s_eqi ( com, 'I' ) ) then call inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, istack, maxrpn, maxsym, & maxval, nints, nsym, nsymp, nsyms, symbol, valsym, value ) return ! ! COM = L, List symbols. ! else if ( s_eqi ( com,'l') ) then if ( namvr == '#' ) then namvr = ' ' call matwrt ( ierror, namvr, ncol, nrow, value ) else call symlst ( ifree, intsym, iopsym, iprsym, ipval, irpn, & maxfrm, maxrpn, maxsym, maxval, namvr, nints, nsym, nsymp, nsyms, & symbol, valsym ) end if return ! ! COM = A, Add variable. ! else if ( s_eqi ( com, 'A' ) ) then call symadd ( ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvr, nsym, nsymp, symbol, valsym ) return ! ! COM = D, Delete variable. ! COM = V, Set variable value. ! COM = R, Get variable value. ! else if ( s_eqi ( com, 'D') .or. s_eqi ( com, 'R' ) .or. & s_eqi ( com, 'V' ) ) then call symbol_value ( com, ierror, ifree, iopsym, iprsym, ipval, maxsym, & maxval, namvr, ncol, nrow, nsym, nsymp, symbol, valsym, value ) return ! ! COM = F or G, compile formula. ! else if ( s_eqi ( com, 'F' ) .or. s_eqi ( com, 'G' ) ) then if ( s_eqi ( com, 'F' ) ) then nuvar = 1 else if ( s_eqi ( com, 'G' ) ) then nuvar = 0 end if ! ! Remove all blanks from the formula. ! call s_blank_delete ( infix ) if ( len_trim ( infix ) <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' The formula has zero length.' return end if ! ! Do a simple check on parentheses. ! if ( .not. s_paren_check ( infix ) ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' Formula does not pass parentheses check.' return end if ! ! Convert string of characters to array of integers which ! are indices of tokens. ! call tokens ( ierror, ifinis, ifree, implic, indx1, indx2, ineg, infix, & intsym, iopsym, iprsym, ipval, maxsym, maxval, & nints, nsym, nsymp, nuvar, symbol, valsym ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' Could not convert formula into tokens.' if ( implic /= 0 ) then write ( *, '(a)' ) & ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' if ( implic == nsym ) then nsym = nsym - 1 end if end if return end if ! ! Convert to RPN. ! imin = ( ifrm - 1 ) * 80 call rpnset ( ierror, ifinis, imin, intsym, iopsym, iprsym, irpn, & istack, maxrpn, maxsym, nints, nrpn, symbol ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' Could not compile formula.' if ( implic /= 0 ) then write ( *, '(a)' ) & ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if infix = ' ' ! ! Check that operators and arguments are balanced. ! ihi = imin + nrpn call rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' RPNCHK fails.' ierror = 1 if ( implic /= 0 ) then write ( *, '(a)' ) & ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if if ( ilo /= imin + 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' RPNCHK fails.' write ( *, '(a,i6)' ) ' IFRM = ', ifrm write ( *, '(a,i6)' ) ' NRPN = ', nrpn write ( *, '(a,i6)' ) ' ILO = ', ilo write ( *, '(a,i6)' ) ' IHI = ', ihi ierror = 1 if ( implic /= 0 ) then write ( *, '(a)' ) & ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' end if return end if ! ! For implicit definition via equality, evaluate formula to ! get dimensions. ! if ( implic == 0 ) then return end if else if ( s_eqi ( com, 'E' ) ) then else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' Unknown command "' // trim ( com ) // '"' ierror = 1 return end if ! ! COM = E, Evaluate formula. ! 150 continue imin = (ifrm-1) * 80 + 1 nrpn = 80 ifreesv = ifree call rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, & irad, irpn, istack, maxrpn, maxsym, maxval, & nrpn, nsym, nsyms, symbol, valsym, value ) ifree = ifreesv if ( ierror /= 0 ) then value = 0.0E+00 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' The formula could not be evaluated!' if ( implic /= 0 ) then write ( *, '(a)' ) ' Cancelling the variable ' // trim ( symbol(implic) ) symbol(implic) = 'VOID' if ( implic == nsym ) then nsym = nsym - 1 end if end if end if return end subroutine digit_to_ch ( digit, c ) ! !******************************************************************************* ! !! DIGIT_TO_CH returns the character representation of a decimal digit. ! ! ! Example: ! ! DIGIT C ! ----- --- ! 0 '0' ! 1 '1' ! ... ... ! 9 '9' ! 17 '*' ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer DIGIT, the digit value between 0 and 9. ! ! Output, character C, the corresponding character, or '*' if DIGIT ! was illegal. ! implicit none ! character c integer digit ! if ( 0 <= digit .and. digit <= 9 ) then c = char ( digit + 48 ) else c = '*' end if return end subroutine ekf_add ( ef, ek, gf, gk, ielem, kind, maxelm, maxnod, maxnpe, ml, & mu, nnode, nodes, nrow ) ! !******************************************************************************* ! !! EKF_ADD adds the element arrays to the global arrays. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, real EF(MAXNPE), the element force vector. ! ! Input, real EK(MAXNPE,MAXNPE), the element stiffness matrix. ! ! Input/output, real GF(MAXNOD), the global force vector. ! ! Input/output, real GK(NROW,NNODE), the global stiffness matrix. ! ! Input, integer IELEM, the index of the element. ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Input, integer NROW, the leading dimension of the global stiffness matrix. ! implicit none ! integer maxelm integer maxnod integer maxnpe integer nnode integer nrow ! real ef(maxnpe) real ek(maxnpe,maxnpe) real gf(maxnod) real gk(nrow,nnode) integer ie integer ielem integer ig integer irow integer je integer jg integer kind(maxelm) integer ml integer mu integer nodes(maxnpe,maxelm) integer norder ! ! The order of the local stiffness matrix and the local force vector ! is one more than the degree of the element. ! norder = kind(ielem) + 1 do ie = 1, norder ! ! IG is the global index of the unknown whose local index is IE. ! ig = nodes(ie,ielem) ! ! Add the local right hand side to the global right hand side. ! gf(ig) = gf(ig) + ef(ie) ! ! Add the local stiffness to the global stiffness. ! do je = 1, norder jg = nodes(je,ielem) irow = ig - jg + ml + mu + 1 gk(irow,jg) = gk(irow,jg) + ek(ie,je) end do end do return end subroutine ekf_set ( dpsi, ef, ek, ielem, infix, irpn, kind, mat, maxelm, & maxfix, maxgauss, maxnod, maxnpe, maxrpn, ngauss, nodes, psi, wgauss, x, & xgauss ) ! !******************************************************************************* ! !! EKF_SET computes the local stiffness matrix EK, and force vector EF. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Workspace, real DPSI(MAXNPE), shape function derivatives. ! ! Output, real EF(MAXNPE), the element force vector. ! ! Output, real EK(MAXNPE,MAXNPE), the element stiffness matrix. ! ! Input, integer IELEM, the index of the element being processed. ! ! Input, character ( len = * ) INFIX, ? ! ! Input, integer IRPN(MAXRPN), ? ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Input, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Workspace, real PSI(MAXNPE), shape function values. ! ! Workspace, real WGAUSS(MAXGAUSS), quadrature weights. ! ! Input, real X(MAXNOD), the X coordinates of the nodes. ! ! Workspace, real XGAUSS(MAXGAUSS), quadrature abscissas. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxnpe integer maxrpn ! real dpsi(maxnpe) real dxdxi real ef(maxnpe) real ek(maxnpe,maxnpe) integer ie integer ielem integer ig1 integer ig2 character ( len = * ) infix integer irpn(maxrpn) integer je integer kind(maxelm) integer l integer mat(maxelm) integer mater integer ngauss integer nodes(maxnpe,maxelm) integer norder real psi(maxnpe) real wgauss(maxgauss) real x(maxnod) real xb real xc real xf real xgauss(maxgauss) real xi real xk real xx ! norder = kind(ielem) + 1 ! ! Get the endpoints of the element. ! ig1 = nodes(1,ielem) ig2 = nodes(norder,ielem) ! ! Get the "stretch" factor for going from the reference interval ! [ -1.0, 1.0E+00 ] to the physical interval [ X(IG1), X(IG2) ]. ! dxdxi = 0.5E+00 * ( x(ig2) - x(ig1) ) ! ! Initialize the local element arrays. ! ef(1:norder) = 0.0E+00 ek(1:norder,1:norder) = 0.0E+00 ! ! Begin the integration point loop. ! do l = 1, ngauss xi = xgauss(l) xx = x(ig1) + ( 1.0E+00 + xi ) * dxdxi mater = mat(ielem) ! ! Evaluate the material functions. ! call material_evaluate ( infix, irpn, mater, maxfix, maxrpn, xb, xc, xf, & xk, xx ) ! ! Evaluate the shape functions. ! call shape ( dpsi, norder, psi, xi ) do ie = 1, norder ef(ie) = ef(ie) + psi(ie) * xf * wgauss(l) * dxdxi do je = 1, norder ek(ie,je) = ek(ie,je) + wgauss(l) * dxdxi * & ( xk * dpsi(ie) / dxdxi * dpsi(je) / dxdxi & + xc * psi(ie) * dpsi(je) / dxdxi & + xb * psi(ie) * psi(je) ) end do end do end do return end subroutine ekf_write ( ef, ek, maxelm, maxnpe, nel, nodes, norder ) ! !******************************************************************************* ! !! EKF_WRITE writes out an element stiffness matrix and right hand side. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real EF(MAXNPE), the element force vector. ! ! Input, real EK(MAXNPE,MAXNPE), the element stiffness matrix. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NEL, the index of the element. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Input, integer NORDER, the polynomial order of the element. ! implicit none ! integer maxelm integer maxnpe integer norder ! real ef(maxnpe) real ek(maxnpe,maxnpe) integer i integer j integer nel integer nodes(maxnpe,maxelm) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EKF_WRITE:' write ( *, '(a,i6)' ) ' Element matrix and force vector for element ', nel write ( *, '(a)' ) ' ' write ( *, '(9(3x,i6,3x))' ) nodes(1:norder,nel) write ( *, '(a)' ) ' ' do i = 1, norder write ( *, '(i3,10g12.4)' ) nodes(i,nel), ( ek(i,j), j = 1, norder), ef(i) end do return end subroutine element_read ( ierror, kind, mat, maxelm, maxnod, maxnpe, nelem, & nmat, nnode, nodes ) ! !******************************************************************************* ! !! ELEMENT_READ reads the element data. ! ! ! Modified: ! ! 26 May 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred and the data was not properly set. ! ! Output, integer KIND(MAXELM), the degree of each element. ! ! Output, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Output, integer NELEM, the number of elements. ! ! Output, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Output, integer NNODE, the number of nodes. ! ! Output, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! implicit none ! integer maxelm integer maxnpe ! integer i integer ierror integer ihi integer ilo character isay integer j integer jhi integer jlo integer kind(maxelm) integer kindi integer lastn integer mat(maxelm) integer mati integer maxnod integer nelem integer nelem1 integer nlast integer nmat integer nnode integer nodes(maxnpe,maxelm) integer numnod logical s_eqi ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENT_READ:' write ( *, '(a)' ) ' Set up the elements.' write ( *, '(a)' ) ' ' nelem = 0 lastn = 1 do nelem = nelem + 1 write ( *, '(a,i6)' ) 'Set up element ', nelem ! ! Get the polynomial degree. ! do write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Enter polynomial degree between 1 and ', maxnpe-1 read ( *, * ) kindi if ( 1 <= kindi .and. kindi <= maxnpe-1 ) then exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENT_READ - Warning!' write ( *, '(a)' ) ' Your input was not acceptable!' end do ! ! Get the material number. ! if ( nmat == 1 ) then mati = 1 else do write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Enter material number (between 1 and ', nmat read ( *, * ) mati if ( 1 <= mati .and. mati <= nmat ) then exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENT_READ - Warning!' write ( *, '(a)' ) ' Your input was not acceptable!' write ( *, '(a)' ) ' Please try again!' end do end if ! ! Store the information. ! numnod = kindi + 1 if ( ( lastn + kindi ) > maxnod ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENT_READ - Fatal error!' write ( *, '(a)' ) ' Maximum node number exceeded.' ierror = 1 return end if do i = 1, numnod nodes(i,nelem) = lastn - 1 + i end do ihi = maxnpe if ( numnod < ihi ) then ilo = numnod + 1 do i = ilo, ihi nodes(i,nelem) = 0 end do end if lastn = lastn - 1 + numnod kind(nelem) = kindi mat(nelem) = mati ! ! See if there are more elements just like this one. ! if ( nelem >= maxelm ) then exit end if nelem1 = nelem + 1 70 continue write ( *, '(a)' ) 'Enter Y if all the elements have been defined.' read ( *, '(a)' ) isay if ( s_eqi ( isay, 'Y' ) ) then exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'If there are more elements in a row, ' write ( *, '(a)' ) 'like this one,' write ( *, '(a)' ) 'type the number of the LAST such element.' write ( *, '(a,i6,a,i6)' ) 'between ', nelem1,' and ', maxelm write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'If there are no more elements, type ', nelem write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter element number.' read ( *, * ) nlast if ( nlast <= nelem ) then cycle end if if ( nlast > maxelm ) then write ( *, '(a)' ) 'Too many elements!' go to 70 end if do i = nelem1, nlast nelem = i kind(nelem) = kindi mat(nelem) = mati if ( (lastn+kindi) > maxnod ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Maximum node number exceeded.' ierror = 1 return end if do j = 1, numnod nodes(j,nelem) = nodes(j,nelem-1)+kindi end do jhi = kindi + 1 if ( numnod < jhi ) then jlo = numnod + 1 do j = jlo, jhi nodes(j,nelem) = 0 end do end if lastn = lastn + kindi end do if ( nelem >= maxelm ) then exit end if write ( *, '(a)' ) 'Enter Y if all the elements have been defined.' read ( *, '(a)' ) isay if ( s_eqi ( isay, 'Y' ) ) then exit end if end do ! ! Compute the number of nodes. ! nnode = 1 + sum ( kind(1:nelem) ) return end subroutine element_write ( kind, mat, maxelm, maxnpe, nelem, nodes ) ! !******************************************************************************* ! !! ELEMENT_WRITE prints out information about the elements. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! implicit none ! integer maxelm integer maxnpe ! integer i integer j integer kind(maxelm) integer mat(maxelm) integer nelem integer nodes(maxnpe,maxelm) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) & 'Element Degree Material Node1 Node2 Node3 Node4 Node5 Node6' write ( *, '(a)' ) ' ' do i = 1, nelem write ( *, '(i5,i7,i9,6i6)' ) i, kind(i), mat(i), & ( nodes(j,i), j = 1, kind(i) + 1 ) end do return end subroutine error_write ( dpsi, gu, ierror, iexact, infix, irpn, kind, mat, & maxelm, maxfix, maxgauss, maxmat, maxnod, maxnpe, maxrpn, nelem, ngauss, & nodes, psi, wgauss, x, xgauss ) ! !******************************************************************************* ! !! ERROR_WRITE prints the RMS and energy norm errors. ! ! ! Discussion: ! ! It is assumed that the user has supplied formulas for the ! exact solution U(X) and its derivative dU/dX. These will be ! compares witht he computed solution Uh(X) and its derivative ! dUh/dX. ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Workspace, real DPSI(MAXNPE), the derivatives of shape functions. ! ! Input, integer IEXACT, records what is known: ! 0, no formulas for the exact solution are known; ! 1, a formula for the exact value of U(X) is known; ! 2, formulas for the exact values of U(X) and dU/dX are known. ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXMAT, the maximum number of materials. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Input, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Workspace, real WGAUSS(MAXGAUSS), quadrature weights. ! ! Workspace, real XGAUSS(MAXGAUSS), quadrature abscissas. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxnpe integer maxrpn ! character com real dpsi(maxnpe) real dudx real duhdx real eint real enorm real gu(maxnod) integer i integer ierror integer iexact integer ifrm integer indx character ( len = * ) infix integer irpn(maxrpn) integer k integer kind(maxelm) integer mat(maxelm) integer mater integer maxmat integer n character ( len = 10 ) namvar integer nel integer nelem integer ngauss integer nodes(maxnpe,maxelm) real psi(maxnpe) real sqint real sqnorm real uh real ux real value real xgauss(maxgauss) real x(maxnod) real x1 real x2 real xb real xc real xf real xi real xk real xx real wgauss(maxgauss) ! if ( iexact < 2 ) then ierror = 1 write ( *, '(a)' ) 'ERROR_WRITE: Error!' write ( *, '(a)' ) ' You must define formulas for ' write ( *, '(a)' ) ' U, the exact solution, and' write ( *, '(a)' ) ' dUdX, the derivative of U,' write ( *, '(a)' ) ' before you can compute the errror.' return end if ifrm = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Error computations' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Element L2-norm Energy-norm' write ( *, '(a)' ) ' ' sqnorm = 0.0E+00 enorm = 0.0E+00 do i = 1, nelem sqint = 0.0E+00 eint = 0.0E+00 xi = -1.0E+00 indx = nodes(1,i) x1 = x(indx) n = kind(i) + 1 indx = nodes(n,i) x2 = x(indx) nel = i do k = 1, ngauss xi = xgauss(k) xx = x1 + 0.5E+00 * ( 1.0E+00 + xi ) * ( x2 - x1 ) ! ! Set the value of X. ! com = 'V' namvar = 'X' value = xx call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ! ! Evaluate the exact solution UX. ! com = 'E' ifrm = maxmat * 4 + 1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) ux = value ! ! Evaluate the exact solution derivative DUDX. ! com = 'E' ifrm = maxmat * 4 + 2 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) dudx = value ! ! Evaluate the computed solution UH and DUHDX. ! call solution_evaluate ( dpsi, duhdx, gu, kind, maxelm, maxnod, maxnpe, & nel, nodes, psi, uh, x, xi ) ! ! Evaluate the material functions B, C, F, K. ! mater = mat(i) call material_evaluate ( infix, irpn, mater, maxfix, maxrpn, xb, xc, & xf, xk, xx ) eint = eint + 0.5E+00 * wgauss(k) * ( x2 - x1 ) * abs ( -xk * & ( dudx - duhdx )**2 + xb * ( ux - uh )**2 ) sqint = sqint + 0.5E+00 * wgauss(k) * ( x2 - x1 ) * ( ux - uh )**2 end do enorm = enorm + eint eint = sqrt ( eint ) sqnorm = sqnorm + sqint sqint = sqrt ( sqint ) write ( *, '(i4,1x,2g14.6)' ) i, sqint, eint end do sqnorm = sqrt ( sqnorm ) enorm = sqrt ( enorm ) write ( *, '(a)' ) ' ' write ( *, '(a,2g14.6)' ) 'Total ', sqnorm, enorm return end subroutine exact_read ( ierror, iexact, infix, inform, irpn, maxfix, maxmat, & maxrpn ) ! !******************************************************************************* ! !! EXACT_READ reads formulas for the exact solution U(x) and dU/dx. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error, the formulas were read properly. ! nonzero, an error occurred. ! ! Output, integer IEXACT, records what is known: ! 0, no formulas for the exact solution are known; ! 1, a formula for the exact value of U(X) is known; ! 2, formulas for the exact values of U(X) and dU/dX are known. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXMAT, the maximum number of materials. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! implicit none ! integer maxfix integer maxrpn ! character com integer ierror integer iexact integer ifrm integer ihi integer ilo character ( len = * ) infix character ( len = * ) inform integer irpn(maxrpn) integer maxmat character ( len = 10 ) namvar real value ! iexact = 0 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Set up the exact solution U and derivative.' write ( *, '(a)' ) 'Enter a formula for the exact solution U.' read ( *, '(a)' ) infix ilo = 4 * maxmat*80+1 ihi = 4 * maxmat*80+80 inform(ilo:ihi) = infix if ( infix == ' ' ) return com = 'F' ifrm = maxmat * 4 + 1 namvar = ' ' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) if ( ierror /= 0 ) then return end if iexact = 1 ! ! Get the formula for dUdX. ! write ( *, '(a)' ) 'Enter a formula for dUdX' write ( *, '(a)' ) 'or RETURN to leave it blank.' read ( *, '(a)' ) infix ilo = (4*maxmat+1)*80+1 ihi = (4*maxmat+1)*80+80 inform(ilo:ihi) = infix if ( infix == ' ' ) then return end if com = 'F' ifrm = maxmat*4+2 namvar = ' ' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, namvar, value ) if ( ierror /= 0 ) then return end if iexact = 2 return end subroutine exact_write ( iexact, inform, maxmat ) ! !******************************************************************************* ! !! EXACT_WRITE prints the formulas for the exact solution U and its derivative. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IEXACT, records what is known: ! 0, no formulas for the exact solution are known; ! 1, a formula for the exact value of U(X) is known; ! 2, formulas for the exact values of U(X) and dU/dX are known. ! ! Input, integer MAXMAT, the maximum number of materials. ! implicit none ! integer iexact integer ihi integer ilo character ( len = * ) inform integer maxmat ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Exact solution information.' write ( *, '(a)' ) ' ' if ( iexact == 0 ) then write ( *, '(a)' ) 'No information has been specified.' return end if ilo = (4*maxmat*80)+ 1 ihi = (4*maxmat*80)+ 80 write ( *, '(a,a)' ) 'U(X): ', inform(ilo:ihi) if ( iexact >= 2 ) then ilo = (4*maxmat+1)*80+ 1 ihi = (4*maxmat+1)*80+ 80 write ( *, '(a,a)' ) 'U''(X): ', inform(ilo:ihi) end if return end subroutine funscl ( angle_to_radian, arg1, arg2, ierror, result, sym ) ! !******************************************************************************* ! !! FUNSCL evaluates a scalar function of one or more scalar arguments. ! ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ARG2, ARG2, the values of the arguments of the function. ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Output, real RESULT, the value of the scalar function. ! ! Input, character ( len = * ) SYM, the symbolic name of the ! function to be evaluated. ! implicit none ! real, parameter :: degrad = 3.14159265358979E+00 / 180.0E+00 ! real angle_to_radian real arg1 real arg2 integer i_lcm integer iarg integer iarg1 integer iarg2 integer ierror integer i_gcd real result logical s_eqi real sarg character ( len = * ) sym real temp ! ierror = 0 result = 0.0E+00 ! ! Addition. ! if ( sym == '+' ) then result = arg1 + arg2 ! ! Subtraction. ! else if ( sym == '-' ) then result = arg1 - arg2 ! ! Division. ! else if ( sym == '/' ) then if ( arg2 == 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Attempt to divide by 0!' return end if result = arg1 / arg2 ! ! Multiplication. ! else if ( sym == '*' ) then result = arg1 * arg2 ! ! Exponentiation. ! else if ( sym == '**' .or. sym == '^' ) then if ( arg1 == 0.0E+00 .and. arg2 == 0.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Attempt to compute 0**0 !' ierror = 1 else if ( arg1 < 0.0E+00 .and. real(int(arg2)) /= arg2 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Illegal exponentiation:' write ( *, '(a,g14.6,a,g14.6)' ) ' ', arg1, ' ** ', arg2 ierror = 1 else sarg = 1.0E+00 if ( arg1 < 0.0E+00 ) then iarg = int ( arg2 ) if ( 2 * ( iarg / 2 ) /= iarg ) then sarg = - 1.0E+00 end if end if result = sarg * abs ( arg1 )**arg2 end if ! ! Assignment. ! else if ( sym == '=' ) then result = arg2 ! ! Absolute value. ! else if ( s_eqi ( sym, 'ABS' ) ) then result = abs ( arg1 ) ! ! Arc Cosine. ! else if ( s_eqi ( sym, 'ACOS' ) ) then if ( arg1 < -1.0E+00 .or. arg1 > 1.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal inverse cosine of ', arg1 ierror = 1 return end if result = acos ( arg1 ) / angle_to_radian ! ! Arc Sine. ! else if ( s_eqi ( sym, 'ASIN' ) ) then if ( arg1 < - 1.0E+00 .or. arg1 > 1.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal inverse sine of ', arg1 ierror = 1 end if result = asin ( arg1 ) / angle_to_radian ! ! Arc Tangent. ! else if ( s_eqi ( sym, 'ATAN' ) ) then result = atan ( arg1 ) / angle_to_radian ! ! Arc Tangent of ratio. ! else if ( s_eqi ( sym, 'ATAN2' ) ) then if ( arg1 == 0.0E+00 .and. arg2 == 0.0E+00 ) then result = 0.0E+00 else result = atan2 ( arg1, arg2 ) / angle_to_radian end if ! ! Cosine. ! else if ( s_eqi ( sym, 'COS' ) ) then result = cos ( angle_to_radian * arg1 ) ! ! Cotangent of angle. ! else if ( s_eqi ( sym, 'COT' ) ) then result = cos ( angle_to_radian * arg1 ) / sin ( angle_to_radian * arg1 ) ! ! Hyperbolic cosine. ! else if ( s_eqi ( sym, 'COSH' ) ) then result = cosh ( angle_to_radian * arg1 ) ! ! Cosecant of angle. ! else if ( s_eqi ( sym, 'CSC' ) ) then result = 1.0E+00 / sin ( angle_to_radian * arg1 ) ! ! Determinant (of a scalar). ! else if ( s_eqi ( sym, 'DET' ) ) then result = arg1 ! ! Diagonal (of a scalar). ! else if ( s_eqi ( sym, 'DIAG' ) ) then result = arg1 ! ! Exponential. ! else if ( s_eqi ( sym, 'EXP' ) ) then result = exp ( arg1 ) ! ! Greatest common factor. ! else if ( s_eqi ( sym, 'GCD' ) ) then iarg1 = nint ( arg1 ) iarg2 = nint ( arg2 ) result = real ( i_gcd ( iarg1, iarg2 ) ) ! ! Hilbert. ! else if ( s_eqi ( sym, 'HILBERT' ) ) then result = 1.0E+00 ! ! Hilbert inverse. ! else if ( s_eqi ( sym, 'HILBINV' ) ) then result = 1.0E+00 ! ! Identity. ! else if ( s_eqi ( sym, 'ID' ) ) then result = 1.0E+00 ! ! Integer value. ! else if ( s_eqi ( sym, 'INT' ) ) then result = real ( int ( arg1 ) ) ! ! Inverse (of a scalar). ! else if ( s_eqi ( sym, 'INV' ) ) then if ( arg1 == 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Attempt to compute 1/0!' return end if result = 1.0E+00 / arg1 ! ! Least common multiple. ! else if ( s_eqi ( sym, 'LCM' ) ) then iarg1 = nint ( arg1 ) iarg2 = nint ( arg2 ) result = real ( i_lcm ( iarg1, iarg2 ) ) ! ! Natural logarithm. ! else if ( s_eqi ( sym, 'ALOG' ) .or. s_eqi ( sym, 'LN' ) .or. & s_eqi ( sym, 'LOG' ) ) then if ( arg1 <= 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal LOG of ', arg1 return end if result = log ( arg1 ) ! ! Logarithm base 10. ! else if ( s_eqi ( sym, 'ALOG10' ) .or. s_eqi ( sym, 'LOG10' ) ) then if ( arg1 <= 0.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal LOG10 of ', arg1 ierror = 1 return end if result = log10 ( arg1 ) ! ! Logarithm base 2. ! else if ( s_eqi ( sym, 'LOG2' ) ) then if ( arg1 <= 0.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal LOG2 of ', arg1 ierror = 1 return end if result = log ( arg1 ) / log ( 2.0E+00 ) ! ! Maximum of two values. ! else if ( s_eqi ( sym, 'MAX' ) ) then if ( arg1 > arg2 ) then result = arg1 else result = arg2 end if ! ! Minimum of two values. ! else if ( s_eqi ( sym, 'MIN' ) ) then if ( arg1 < arg2 ) then result = arg1 else result = arg2 end if ! ! Negation ! else if ( s_eqi ( sym, 'NEG' ) ) then result = - arg1 ! ! Nearest integer value ! else if ( s_eqi ( sym, 'NINT' ) ) then result = anint ( arg1 ) ! ! NORM0, NORM1, NORM2, NORME or NORMF of a scalar. ! else if ( s_eqi ( sym, 'NORM0' ) .or. s_eqi ( sym, 'NORM1' ) .or. & s_eqi ( sym, 'NORM2' ) .or. s_eqi ( sym, 'NORME' ) .or. & s_eqi ( sym, 'NORMF' ) ) then result = abs ( arg1 ) ! ! POLY (scalar coefficient array, means constant polynomial). ! else if ( s_eqi ( sym, 'POLY' ) ) then result = arg1 ! ! Random value. ! else if ( s_eqi ( sym, 'RAN' ) ) then call random_number ( result ) ! ! Rounding. ! else if ( s_eqi ( sym, 'ROUND' ) ) then if ( abs ( arg1 ) < epsilon ( arg1 ) ) then result = 0.0E+00 else result = arg1 end if ! ! Secant. ! else if ( s_eqi ( sym, 'SEC' ) ) then result = 1.0E+00 / cos ( angle_to_radian * arg1 ) ! ! Sine. ! else if ( s_eqi ( sym, 'SIN' ) ) then result = sin ( angle_to_radian * arg1 ) ! ! Hyperbolic sine. ! else if ( s_eqi ( sym, 'SINH' ) ) then result = sinh ( angle_to_radian * arg1 ) ! ! Square root. ! else if ( s_eqi ( sym, 'SQRT' ) ) then if ( arg1 < 0.0E+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a,g14.6)' ) ' Illegal SQRT of ', arg1 return end if result = sqrt ( arg1 ) ! ! Step function. ! else if ( s_eqi ( sym, 'STEP' ) ) then if ( arg1 < 0.0E+00 ) then result = 0.0E+00 else result = 1.0E+00 end if ! ! Tangent. ! else if ( s_eqi ( sym, 'TAN' ) ) then result = tan ( angle_to_radian * arg1 ) ! ! Hyperbolic tangent. ! else if ( s_eqi ( sym, 'TANH' ) ) then result = tanh ( angle_to_radian * arg1 ) ! ! Trace (of a scalar) ! else if ( s_eqi ( sym, 'TRACE' ) ) then result = arg1 ! ! Transpose (of a scalar) ! else if ( s_eqi ( sym, 'TRANS' ) ) then result = arg1 ! ! Zero ! else if ( s_eqi ( sym, 'ZEROS' ) ) then result = 0.0E+00 ! ! Factorial function. ! else if ( sym == '!' ) then if ( arg1 < 0.0E+00 .or. arg1 > 25.0E+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Argument out of range for factorial!' ierror = 1 return end if temp = arg1 result = 1.0E+00 do while ( temp > 1.0E+00 ) result = result * temp temp = temp - 1.0E+00 end do ! ! Unknown function. ! else result = 0.0E+00 ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FUNSCL - Error!' write ( *, '(a)' ) ' Unrecognized function name: "' // trim ( sym ) // '"' end if return end subroutine funval ( iarg1, iarg2, ierror, ifree, iopsym, iprsym, ipset, & ipval, irad, itemp, maxsym, maxval, nsyms, sym, symbol, valsym ) ! !******************************************************************************* ! !! FUNVAL evaluates a function given an argument and an RPN formula. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, INTEGER IARG1, IARG2, IARG3, are the indices ! of the arguments to the current function. Actually, ! the function might have just one or two arguments, ! in which case IARG2 and IARG3 may be zero. ! ! Output, INTEGER IERROR. ! ! IERROR is an error flag. If it is zero on return, ! then no error was detected, and the formula was ! evaluated successfully. If it is nonzero on return, ! then an error was found, and the evaluation of the ! formula could not be carried out. ! ! Input, INTEGER IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! implicit none ! integer, parameter :: maxlen = 10 ! integer maxsym integer maxval ! real angle_to_radian real arg1 real arg2 character ( len = 3 ) ctemp integer i integer iarg1 integer iarg2 integer ierror integer ifree integer index1 integer index4 integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer itemp integer jseed logical s_eqi integer ncol1 integer ndim1 integer ndim2 integer nrow1 integer nsyms real result character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! ipset = 0 ierror = 0 ! ! If the operator is assignment, then the user may be trying to ! reset the values of certain reserved variables, including ! PI and EPS. ! if ( sym == ' = ' ) then if ( s_eqi ( symbol(iarg1),'E').or. & s_eqi ( symbol(iarg1),'EPS').or. & s_eqi ( symbol(iarg1),'PI') ) then ierror = 1 write ( *, '(a)' ) 'Error! You may not change the value of ' & // trim ( symbol(iarg1) ) return end if end if ! ! Set dimensions for HILBERT, HILBINV, ID and ZERO functions, ! which otherwise would look like scalar functions and get ! trapped in FUNSCL. ! if ( & s_eqi ( sym,'HILBERT').or. & s_eqi ( sym,'HILBINV').or. & s_eqi ( sym,'ID') .or. & s_eqi ( sym,'ZERO') ) then nrow1 = int(valsym(ipval(iarg1))) ncol1 = nrow1 else nrow1 = 1 ncol1 = 1 end if ndim1 = nrow1 * ncol1 index1 = ipval(iarg1) ndim2 = 0 if ( nsyms >= maxsym ) then ierror = 1 write ( *, '(a)' ) 'Error! Not enough free memory left!' write ( *, '(a)' ) 'Perhaps the KILL or INIT commands are needed!' return end if nsyms = nsyms + 1 iopsym(nsyms) = 0 iprsym(nsyms) = 10 itemp = itemp + 1 ctemp = ' ' call i_to_s_zero ( itemp, ctemp ) symbol(nsyms) = 'STK000' symbol(nsyms)(4:6) = ctemp ipval(nsyms) = ifree index4 = ipval(nsyms) ! ! Determine the index of the SEED function, in case the ! RANDOM function is called. ! do i = 1, maxsym if ( s_eqi ( symbol(i),'seed' ) ) then jseed = i end if end do ! ! Check for simple scalar arithmetic. ! if ( (ndim1<=1 .and. iarg2 == 0).or.(ndim1<=1 .and. ndim2<=1) ) then arg1 = valsym(index1) arg2 = 0.0E+00 if ( iarg2 /= 0 ) then arg2 = valsym(ipval(iarg2)) end if angle_to_radian = valsym(irad) call funscl ( angle_to_radian, arg1, arg2, ierror, result, sym ) ! ! Store the current value of ISEED, which will have changed ! if the RANDOM function was invoked. ! ! valsym(ipval(jseed)) = real(iseed) valsym(index4) = result if ( sym == ' = ' ) then arg1 = arg2 valsym(index1) = arg2 end if return else write ( *, '(a)' ) 'Illegal matrix operation (I think)!' write ( *, '(a,i6)' ) 'NDIM1 = ', ndim1 write ( *, '(a,i6)' ) 'NDIM2 = ', ndim2 write ( *, '(a,i6)' ) 'IARG1 = ', iarg1 write ( *, '(a,i6)' ) 'IARG2 = ', iarg2 stop end if end subroutine gauss_read ( ierror, maxgauss, ngauss ) ! !******************************************************************************* ! !! GAUSS_READ reads the desired quadrature rule order. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order allowed. ! ! Output, integer NGAUSS, the desired quadrature rule order. ! implicit none ! integer ierror integer maxgauss integer ngauss integer ntemp ! write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Choose a quadrature rule order between 1 and ', & maxgauss read ( *, * ) ntemp if ( ntemp < 1 .or. ntemp > maxgauss ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GAUSS_READ - Fatal error!' write ( *, '(a)' ) ' Your choice was not acceptable.' else ierror = 0 ngauss = ntemp end if return end subroutine gauss_set ( ngauss, wgauss, xgauss ) ! !******************************************************************************* ! !! GAUSS_SET sets the weights and abscissas for Gaussian quadrature. ! ! ! Modified: ! ! 17 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NGAUSS, the number of points to use in the quadrature. ! 1 <= NGAUSS <= 10. ! ! Output, real WGAUSS(NGAUSS), XGAUSS(NGAUSS), the weights and abscissas ! for the Gauss quadrature rule. ! implicit none ! integer ngauss ! real wgauss(ngauss) real xgauss(ngauss) ! if ( ngauss == 1 ) then xgauss(1) = 0.0E+00 wgauss(1) = 2.0E+00 else if ( ngauss == 2 ) then xgauss(1) = - 0.577350269189625764509148780502E+00 xgauss(2) = 0.577350269189625764509148780502E+00 wgauss(1) = 1.0E+00 wgauss(2) = 1.0E+00 else if ( ngauss == 3 ) then xgauss(1) = - 0.774596669241483377035853079956E+00 xgauss(2) = 0.0E+00 xgauss(3) = 0.774596669241483377035853079956E+00 wgauss(1) = 5.0E+00 / 9.0E+00 wgauss(2) = 8.0E+00 / 9.0E+00 wgauss(3) = 5.0E+00 / 9.0E+00 else if ( ngauss == 4 ) then xgauss(1) = - 0.861136311594052575223946488893E+00 xgauss(2) = - 0.339981043584856264802665759103E+00 xgauss(3) = 0.339981043584856264802665759103E+00 xgauss(4) = 0.861136311594052575223946488893E+00 wgauss(1) = 0.347854845137453857373063949222E+00 wgauss(2) = 0.652145154862546142626936050778E+00 wgauss(3) = 0.652145154862546142626936050778E+00 wgauss(4) = 0.347854845137453857373063949222E+00 else if ( ngauss == 5 ) then xgauss(1) = - 0.906179845938663992797626878299E+00 xgauss(2) = - 0.538469310105683091036314420700E+00 xgauss(3) = 0.0E+00 xgauss(4) = 0.538469310105683091036314420700E+00 xgauss(5) = 0.906179845938663992797626878299E+00 wgauss(1) = 0.236926885056189087514264040720E+00 wgauss(2) = 0.478628670499366468041291514836E+00 wgauss(3) = 0.568888888888888888888888888889E+00 wgauss(4) = 0.478628670499366468041291514836E+00 wgauss(5) = 0.236926885056189087514264040720E+00 else if ( ngauss == 6 ) then xgauss(1) = - 0.932469514203152027812301554494E+00 xgauss(2) = - 0.661209386466264513661399595020E+00 xgauss(3) = - 0.238619186083196908630501721681E+00 xgauss(4) = 0.238619186083196908630501721681E+00 xgauss(5) = 0.661209386466264513661399595020E+00 xgauss(6) = 0.932469514203152027812301554494E+00 wgauss(1) = 0.171324492379170345040296142173E+00 wgauss(2) = 0.360761573048138607569833513838E+00 wgauss(3) = 0.467913934572691047389870343990E+00 wgauss(4) = 0.467913934572691047389870343990E+00 wgauss(5) = 0.360761573048138607569833513838E+00 wgauss(6) = 0.171324492379170345040296142173E+00 else if ( ngauss == 7 ) then xgauss(1) = - 0.949107912342758524526189684048E+00 xgauss(2) = - 0.741531185599394439863864773281E+00 xgauss(3) = - 0.405845151377397166906606412077E+00 xgauss(4) = 0.0E+00 xgauss(5) = 0.405845151377397166906606412077E+00 xgauss(6) = 0.741531185599394439863864773281E+00 xgauss(7) = 0.949107912342758524526189684048E+00 wgauss(1) = 0.129484966168869693270611432679E+00 wgauss(2) = 0.279705391489276667901467771424E+00 wgauss(3) = 0.381830050505118944950369775489E+00 wgauss(4) = 0.417959183673469387755102040816E+00 wgauss(5) = 0.381830050505118944950369775489E+00 wgauss(6) = 0.279705391489276667901467771424E+00 wgauss(7) = 0.129484966168869693270611432679E+00 else if ( ngauss == 8 ) then xgauss(1) = - 0.960289856497536231683560868569E+00 xgauss(2) = - 0.796666477413626739591553936476E+00 xgauss(3) = - 0.525532409916328985817739049189E+00 xgauss(4) = - 0.183434642495649804939476142360E+00 xgauss(5) = 0.183434642495649804939476142360E+00 xgauss(6) = 0.525532409916328985817739049189E+00 xgauss(7) = 0.796666477413626739591553936476E+00 xgauss(8) = 0.960289856497536231683560868569E+00 wgauss(1) = 0.101228536290376259152531354310E+00 wgauss(2) = 0.222381034453374470544355994426E+00 wgauss(3) = 0.313706645877887287337962201987E+00 wgauss(4) = 0.362683783378361982965150449277E+00 wgauss(5) = 0.362683783378361982965150449277E+00 wgauss(6) = 0.313706645877887287337962201987E+00 wgauss(7) = 0.222381034453374470544355994426E+00 wgauss(8) = 0.101228536290376259152531354310E+00 else if ( ngauss == 9 ) then xgauss(1) = - 0.968160239507626089835576202904E+00 xgauss(2) = - 0.836031107326635794299429788070E+00 xgauss(3) = - 0.613371432700590397308702039341E+00 xgauss(4) = - 0.324253423403808929038538014643E+00 xgauss(5) = 0.0E+00 xgauss(6) = 0.324253423403808929038538014643E+00 xgauss(7) = 0.613371432700590397308702039341E+00 xgauss(8) = 0.836031107326635794299429788070E+00 xgauss(9) = 0.968160239507626089835576202904E+00 wgauss(1) = 0.0812743883615744119718921581105E+00 wgauss(2) = 0.180648160694857404058472031243E+00 wgauss(3) = 0.260610696402935462318742869419E+00 wgauss(4) = 0.312347077040002840068630406584E+00 wgauss(5) = 0.330239355001259763164525069287E+00 wgauss(6) = 0.312347077040002840068630406584E+00 wgauss(7) = 0.260610696402935462318742869419E+00 wgauss(8) = 0.180648160694857404058472031243E+00 wgauss(9) = 0.0812743883615744119718921581105E+00 else if ( ngauss == 10 ) then xgauss(1) = - 0.973906528517171720077964012084E+00 xgauss(2) = - 0.865063366688984510732096688423E+00 xgauss(3) = - 0.679409568299024406234327365115E+00 xgauss(4) = - 0.433395394129247290799265943166E+00 xgauss(5) = - 0.148874338981631210884826001130E+00 xgauss(6) = 0.148874338981631210884826001130E+00 xgauss(7) = 0.433395394129247290799265943166E+00 xgauss(8) = 0.679409568299024406234327365115E+00 xgauss(9) = 0.865063366688984510732096688423E+00 xgauss(10) = 0.973906528517171720077964012084E+00 wgauss(1) = 0.0666713443086881375935688098933E+00 wgauss(2) = 0.149451349150580593145776339658E+00 wgauss(3) = 0.219086362515982043995534934228E+00 wgauss(4) = 0.269266719309996355091226921569E+00 wgauss(5) = 0.295524224714752870173892994651E+00 wgauss(6) = 0.295524224714752870173892994651E+00 wgauss(7) = 0.269266719309996355091226921569E+00 wgauss(8) = 0.219086362515982043995534934228E+00 wgauss(9) = 0.149451349150580593145776339658E+00 wgauss(10) = 0.0666713443086881375935688098933E+00 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GAUSS_SET - Fatal error!' write ( *, '(a,i6)' ) ' Illegal Gauss quadrature order = ', ngauss stop end if return end subroutine gauss_write ( maxgauss, ngauss, wgauss, xgauss ) ! !******************************************************************************* ! !! GAUSS_WRITE prints out information about the quadrature rule. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer MAXGAUSS, the maximum quadrature rule order allowed. ! ! Input, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Input, real WGAUSS(NGAUSS), weights of the quadrature rule. ! ! Input, real XGAUSS(NGAUSS), abscissas of the quadrature rule. ! implicit none ! integer maxgauss ! integer i integer ngauss real wgauss(maxgauss) real xgauss(maxgauss) ! write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'The quadrature rule is of of order ', ngauss write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I Weight(I), X(I)' write ( *, '(a)' ) ' ' do i = 1, ngauss write ( *, '(i6,2g14.6)' ) i, wgauss(i), xgauss(i) end do return end subroutine gkf_assemble ( dpsi, ef, ek, gf, gk, infix, irpn, kind, & mat, maxelm, maxfix, maxgauss, maxnod, maxnpe, maxrpn, ml, mu, & nelem, ngauss, nnode, nodes, nrow, psi, wgauss, x, xgauss ) ! !******************************************************************************* ! !! GKF_ASSEMBLE forms the global stiffness matrix GK and force vector GF. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Workspace, real DPSI(MAXNPE), the derivatives of shape functions. ! ! Workspace, real EF(MAXNPE), the element force vector. ! ! Workspace, real EK(MAXNPE,MAXNPE), the element stiffness matrix. ! ! Output, real GF(MAXNOD), the global force vector. ! ! Output, real GK(NROW,NNODE), the global stiffness matrix. ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Input, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Input, integer NROW, the leading dimension of the global stiffness matrix. ! ! Workspace, real WGAUSS(MAXGAUSS), quadrature weights. ! ! Input, real X(MAXNOD), the X coordinates of the nodes. ! ! Workspace, real XGAUSS(MAXGAUSS), quadrature abscissas. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxnpe integer maxrpn integer nnode integer nrow ! real dpsi(maxnpe) real ef(maxnpe) real ek(maxnpe,maxnpe) real gf(maxnod) real gk(nrow,nnode) integer i integer ielem integer imax character ( len = * ) infix integer irpn(maxrpn) integer ishow integer j integer kind(maxelm) integer mat(maxelm) integer ml integer mu integer nel integer nelem integer ngauss integer nodes(maxnpe,maxelm) integer norder real psi(maxnpe) real wgauss(maxgauss) real x(maxnod) real xgauss(maxgauss) ! ! Zero out the global stiffness matrix and force vector. ! gf(1:maxnod) = 0.0E+00 gk(1:nrow,1:nnode) = 0.0E+00 ! ! Compute the global stiffness matrix EK and right hand side EF, ! by constructing the local quantities and adding them up. ! ishow = 0 imax = nelem+1 do nel = 1, nelem ielem = nel if ( ( ishow == 0 .and. nel /= 1 ) .or. ishow == imax .or. & ishow >= nel ) then go to 50 end if write ( *, 1020 ) write ( *, 1021 ) ishow, nelem write ( *, 1022 ) imax, ishow, nelem read ( *, *,end = 40, err = 40 ) ishow go to 50 40 continue ishow = 0 50 continue ! ! Form the element matrix EK and right hand side vector EF. ! call ekf_set ( dpsi, ef, ek, ielem, infix, irpn, kind, mat, & maxelm, maxfix, maxgauss, maxnod, maxnpe, maxrpn, ngauss, & nodes, psi, wgauss, x, xgauss ) ! ! Offer to print EK and EF. ! if ( ishow == nel .or. ishow == imax ) then norder = kind(ielem) + 1 call ekf_write ( ef, ek, maxelm, maxnpe, ielem, nodes, norder ) end if ! ! Add EK and EF to the global items GK and GF. ! call ekf_add ( ef, ek, gf, gk, ielem, kind, maxelm, & maxnod, maxnpe, ml, mu, nnode, nodes, nrow ) end do return 1020 format(' Enter 0 to display no more element matrices') 1021 format(' or number between ',i6,' and ',i6) 1022 format(' or enter ',i6,' for matrices ',i6,' to ',i6) end subroutine gkf_width ( kind, maxelm, maxnpe, ml, mu, nelem, nodes ) ! !******************************************************************************* ! !! GKF_WIDTH computes the bandwidth of the stiffness matrix. ! ! ! Modified: ! ! 16 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Output, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! implicit none ! integer maxelm integer maxnpe ! integer i integer irow integer j integer jcol integer kind(maxelm) integer ml integer mu integer nel integer nelem integer norder integer nodes(maxnpe,maxelm) ! ml = 0 mu = 0 do nel = 1, nelem norder = kind(nel) + 1 do i = 1, norder irow = nodes(i,nel) do j = 1, norder jcol = nodes(j,nel) ml = max ( ml, irow - jcol ) mu = max ( mu, jcol - irow ) end do end do end do return end subroutine gkf_write ( gf, gk, ml, mu, nnode, nrow, irhs ) ! !******************************************************************************* ! !! GKF_WRITE prints out the global stiffness matrix and force vector. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real GF(NNODE), the global force vector. ! ! Input, real GK(NROW,NNODE), the global stiffness matrix. ! ! Input, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NROW, the leading dimension of the global stiffness matrix. ! ! Input, integer IRHS, right hand side indicator. ! 0, print GK only. ! nonzero, print GK and GF. ! implicit none ! integer nnode integer nrow ! character ( len = 14 ) ctemp(10) real entry real gf(nnode) real gk(nrow,nnode) integer i integer icopy integer ihi integer iihi integer iilo integer ilo integer inc integer incp1 integer incx integer irhs integer j integer jcopy integer jhi integer jlo integer ml integer mu integer nonzer character ( len = 100 ) output ! if ( irhs == 0 ) then incx = 5 incp1 = incx write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'System matrix' else incx = 4 incp1 = incx + 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'System matrix and right hand side' end if write ( *, '(a)' ) ' ' do jlo = 1, nnode, incx jhi = min ( jlo+incx-1, nnode ) inc = min ( incx, nnode+1-jlo ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Column' ! ! Set up the column headers. ! iilo = 10 iihi = 15 output = ' ' do j = jlo, jhi write ( output(iilo:iihi), '(i6)' ) j iilo = iilo+14 iihi = iihi+14 end do output(68:70) = 'RHS' write ( *, '(a)' ) output ! ! Print the row index column header. ! write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' ilo = max ( jlo-mu, 1 ) ihi = min ( jhi+ml, nnode ) do i = ilo, ihi icopy = i nonzer = 0 do j = 1, inc jcopy = jlo-1+j if ( jcopy - icopy > mu ) then entry = 0 else if ( icopy - jcopy > ml ) then entry = 0 else entry = gk(icopy-jcopy+ml+mu+1,jcopy) end if if ( entry /= 0.0E+00 ) then nonzer = 1 end if if ( -mu > i-jcopy .or. i-jcopy > ml ) then ctemp(j) = ' ' else write ( ctemp(j), '(g14.6)' ) entry end if end do if ( inc < incx ) then do j = inc+1, incp1 ctemp(j) = ' ' end do end if if ( irhs /= 0 ) then write ( ctemp(incp1), '(g14.6)' ) gf(i) end if if ( nonzer == 1 ) then write ( *,'(i4,1x,8a14)') i, (ctemp(j),j = 1,incp1) end if end do end do write ( *, '(a)' ) ' ' return end subroutine hello ( maxelm, maxgauss, maxmat, maxnod, maxnpe, maxpoi ) ! !******************************************************************************* ! !! HELLO prints out an introductory message. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXMAT, the maximum number of materials. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! implicit none ! integer maxelm integer maxgauss integer maxmat integer maxnod integer maxnpe integer maxpoi ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD1' write ( *, '(a)' ) ' A program for solving one dimensional boundary' write ( *, '(a)' ) ' value problems using the finite element method.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Program maxima:' write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Elements = ', maxelm write ( *, '(a,i6)' ) ' Nodes = ', maxnod write ( *, '(a,i6)' ) ' Nodes per element = ', maxnpe write ( *, '(a,i6)' ) ' Materials = ', maxmat write ( *, '(a,i6)' ) ' Quadrature points = ', maxgauss write ( *, '(a,i6)' ) ' Point loads = ', maxpoi write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Last modified on 17 May 1999' write ( *, '(a)' ) ' ' return end subroutine help ! !******************************************************************************* ! !! HELP prints out a list of legal commands. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD1 commands:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Help Print list of commands.' write ( *, '(a)' ) 'Quit Stop the program.' write ( *, '(a)' ) 'Set BC, ELEM, EXACT, GAUSS, MAT, NODE, POINT, TITLE.' write ( *, '(a)' ) 'Solve Solve the system.' write ( *, '(a)' ) 'Type BC, ELEM, ERROR, EXACT, GAUSS, MAT, NODE, ' & // 'POINT, SOL, SYS, TITLE.' write ( *, '(a)' ) '# To make a comment.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PREP Input the data in an orderly fashion.' write ( *, '(a)' ) 'PROC Solve the system.' write ( *, '(a)' ) 'POST Print out the solution data.' write ( *, '(a)' ) ' ' return end function i_gcd ( i, j ) ! !******************************************************************************* ! !! I_GCD finds the greatest common divisor of I and J. ! ! ! Modified: ! ! 03 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, two numbers whose greatest common divisor ! is desired. ! ! Output, integer I_GCD, the greatest common divisor of I and J. ! ! Note that only the absolute values of I and J are ! considered, so that the result is always nonnegative. ! ! If I or J is 0, I_GCD is returned as max ( 1, abs ( I ), abs ( J ) ). ! ! If I and J have no common factor, I_GCD is returned as 1. ! ! Otherwise, using the Euclidean algorithm, I_GCD is the ! largest common factor of I and J. ! implicit none ! integer i integer i_gcd integer ip integer iq integer ir integer j ! i_gcd = 1 ! ! Return immediately if either I or J is zero. ! if ( i == 0 ) then i_gcd = max ( 1, abs ( j ) ) return else if ( j == 0 ) then i_gcd = max ( 1, abs ( i ) ) return end if ! ! Set IP to the larger of I and J, IQ to the smaller. ! This way, we can alter IP and IQ as we go. ! ip = max ( abs ( i ), abs ( j ) ) iq = min ( abs ( i ), abs ( j ) ) ! ! Carry out the Euclidean algorithm. ! do ir = mod ( ip, iq ) if ( ir == 0 ) then exit end if ip = iq iq = ir end do i_gcd = iq return end function i_lcm ( i, j ) ! !******************************************************************************* ! !! I_LCM computes the least common multiple of two integers. ! ! ! Definition: ! ! The least common multiple may be defined as ! ! LCM(I,J) = ABS( I * J ) / GCD(I,J) ! ! where GCD(I,J) is the greatest common divisor of I and J. ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, the integers whose I_LCM is desired. ! ! Output, integer I_LCM, the least common multiple of I and J. ! I_LCM is never negative. I_LCM is 0 if either I or J is zero. ! implicit none ! integer i integer i_gcd integer j integer i_lcm ! i_lcm = abs ( i * ( j / i_gcd ( i, j ) ) ) return end subroutine i_to_s_zero ( intval, s ) ! !******************************************************************************* ! !! I_TO_S_ZERO converts an integer to a string, with zero padding. ! ! ! Examples: ! ! Assume that S is 6 characters long: ! ! INTVAL S ! ! 1 000001 ! -1 -00001 ! 0 000000 ! 1952 001952 ! 123456 123456 ! 1234567 ****** <-- Not enough room! ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer INTVAL, an integer to be converted. ! ! Output, character ( len = * ) S, the representation of the integer. ! The integer will be right justified, and zero padded. ! If there is not enough space, the string will be filled with stars. ! implicit none ! character c integer i integer idig integer ihi integer ilo integer intval integer ipos integer ival character ( len = * ) s ! s = ' ' ilo = 1 ihi = len ( s ) if ( ihi <= 0 ) then return end if ! ! Make a copy of the integer. ! ival = intval ! ! Handle the negative sign. ! if ( ival < 0 ) then if ( ihi <= 1 ) then s(1:1) = '*' return end if ival = - ival s(1:1) = '-' ilo = 2 end if ! ! Working from right to left, strip off the digits of the integer ! and place them into S(ILO:IHI). ! ipos = ihi do while ( ival /= 0 .or. ipos == ihi ) idig = mod ( ival, 10 ) ival = ival / 10 if ( ipos < ilo ) then do i = 1, ihi s(i:i) = '*' end do return end if call digit_to_ch ( idig, c ) s(ipos:ipos) = c ipos = ipos - 1 end do ! ! Fill the empties with zeroes. ! do i = ilo, ipos s(i:i) = '0' end do return end subroutine inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, istack, maxrpn, maxsym, & maxval, nints, nsym, nsymp, nsyms, symbol, valsym, value ) ! !******************************************************************************* ! !! INICOM initializes data for COMRPN. ! ! ! Modified: ! ! 02 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Output, integer IFREE, the index of the first free address in VALSYM. ! ! Output, integer INDX1, the index of the symbol 'INDEX1'. ! ! Output, integer INDX2, the index of the symbol 'INDEX2'. ! ! Output, integer INEG, the index of the symbol 'NEG'. ! ! Output, character ( len = * ) INFIX, space for an infix formula. ! ! Output, integer INTSYM(80), a set of integers which ! are the indices of symbols, representing an infix formula. ! ! Output, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Output, integer IRAD, the index of the symbol 'ANGLE_TO_RADIAN'. ! ! Output, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Workspace, integer ISTACK(MAXSYM), workspace for interpreting the formula. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed in VALSYM. ! ! Output, integer NINTS, the number of useful entries in INTSYM. ! ! Output, integer NSYM, the total number of symbols, including temporaries. ! ! Output, integer NSYMP, the number of permanent symbols. ! ! Output, integer NYMS, the number of declared symbols. ! ! Output, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Output, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Output, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of the variable to be saved. ! implicit none ! integer, parameter :: maxlen = 20 ! integer maxrpn integer maxsym integer maxval ! integer i integer ifinis integer ifree integer indx1 integer indx2 integer ineg character ( len = * ) infix integer intsym(80) integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer irad integer irpn(maxrpn) integer istack(maxsym) integer j integer nints integer nsym integer nsymp integer nsyms logical s_eqi character ( len = maxlen ) symbol(maxsym) real temp real valsym(maxval) real value ! nsym = 0 nsym = nsym + 1 symbol(nsym) = 'ABS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ACOS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ALOG' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ALOG10' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ANGLE_TO_RADIAN' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = 1.0E+00 nsym = nsym + 1 symbol(nsym) = 'ASIN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ATAN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ATAN2' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COSH' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'COT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'CSC' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'DET' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'DIAG' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'E' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = exp ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'EPS' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = epsilon ( 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'EVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'EVEC' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'EXP' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'GCD' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HILBERT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HILBINV' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'HOUSE' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ID' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INV' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'IVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'KRYLOV' iopsym(nsym) = 3 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LCM' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG10' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'LOG2' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'MAX' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'MIN' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NEG' iopsym(nsym) = 1 iprsym(nsym) = 5 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NINT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM0' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM1' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORM2' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORME' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORMF' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'NORMS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'PI' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = 4.0E+00 * atan2 ( 1.0E+00, 1.0E+00 ) nsym = nsym + 1 symbol(nsym) = 'POLY' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RAN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RCOND' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ROUND' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'RVAL' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SEC' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SEED' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SIN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SINH' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'SQRT' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'STEP' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TAN' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TANH' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TRACE' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'TRANS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'ZEROS' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '**' iopsym(nsym) = 2 iprsym(nsym) = 8 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '^' iopsym(nsym) = 2 iprsym(nsym) = 8 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '+' iopsym(nsym) = 2 iprsym(nsym) = 5 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '-' iopsym(nsym) = 2 iprsym(nsym) = 5 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '/' iopsym(nsym) = 2 iprsym(nsym) = 7 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '*' iopsym(nsym) = 2 iprsym(nsym) = 6 valsym(nsym) = 0.0E+00 ! ! Take care of fact that ninny UNIX FORTRAN compilers won't let us ! explicitly use a backslash character! ! nsym = nsym + 1 symbol(nsym) = char ( 92 ) iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '(' iopsym(nsym) = -1 iprsym(nsym) = 1 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = ')' iopsym(nsym) = -1 iprsym(nsym) = 2 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '=' iopsym(nsym) = 2 iprsym(nsym) = 3 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = ',' iopsym(nsym) = -1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '!' iopsym(nsym) = 1 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INDEX1' iopsym(nsym) = 2 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = 'INDEX2' iopsym(nsym) = 3 iprsym(nsym) = 9 valsym(nsym) = 0.0E+00 nsym = nsym + 1 symbol(nsym) = '$' iopsym(nsym) = -1 iprsym(nsym) = 0 valsym(nsym) = 0.0E+00 nsymp = nsym nsyms = nsymp ifree = nsymp + 1 infix = ' ' intsym(1:80) = 0 irpn(1:maxrpn) = 0 istack(1:maxsym) = 0 nints = 0 do i = 1, nsymp ipval(i) = i end do do i = 1, nsymp if ( s_eqi ( symbol(i), 'ANGLE_TO_RADIAN' ) ) then irad = i else if ( s_eqi ( symbol(i), 'NEG' ) ) then ineg = i else if ( symbol(i) == '$' ) then ifinis = i else if ( s_eqi ( symbol(i), 'INDEX1' ) ) then indx1 = i else if ( s_eqi ( symbol(i), 'INDEX2' ) ) then indx2 = i end if end do call random_seed ( ) value = 0.0E+00 return end subroutine inifem ( iexact, inform, kbc, kind, kpoint, mat, maxelm, maxgauss, & maxnod, maxnpe, maxpoi, maxrpn, ml, mu, nelem, ngauss, nmat, nnode, nodes, & npoint, point, title, vbc, wgauss, x, xgauss ) ! !******************************************************************************* ! !! INIFEM initializes the finite element data. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IEXACT, records what is known: ! 0, no formulas for the exact solution are known; ! 1, a formula for the exact value of U(X) is known; ! 2, formulas for the exact values of U(X) and dU/dX are known. ! ! Output, integer KBC(2). KBC(1) and KBC(2) are the boundary condition ! codes at the left and right endpoints of the region. ! 1: U = value1; ! 2: K dU/dx = value1 * U - value2 ! ! Output, integer KIND(MAXELM), the degree of each element. ! ! Output, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXPOI, the maximum number of point loads. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Output, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Output, integer NELEM, the number of elements. ! ! Output, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Output, integer NNODE, the number of nodes. ! ! Output, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Output, integer NPOINT, the number of point loads. ! ! Output, real VBC(2,2), contains constants for boundary conditions. ! VBC(1,1) and VBC(2,1) are value1 and value2 for the left endpoint; ! VBC(1,2) and VBC(2,2) are value1 and value2 for the right endpoint. ! ! Output, real X(MAXNOD), the X coordinates of the nodes. ! implicit none ! integer maxelm integer maxgauss integer maxnod integer maxnpe integer maxpoi ! integer i integer iexact character ( len = * ) inform integer j integer kbc(2) integer kind(maxelm) integer kpoint(maxpoi) integer mat(maxelm) integer maxrpn integer ml integer mu integer nelem integer ngauss integer nmat integer nnode integer nodes(maxnpe,maxelm) integer npoint real point(maxpoi) character ( len = 80 ) title real vbc(2,2) real wgauss(maxgauss) real x(maxnod) real xgauss(maxgauss) ! iexact = 0 do i = 1, maxrpn inform(i:i) = ' ' end do kbc(1:2) = 0 kind(1:maxelm) = 0 kpoint(1:maxpoi) = 0 mat(1:maxelm) = 0 ml = 0 mu = 0 nelem = 0 ngauss = 0 nmat = 0 nnode = 0 nodes(1:maxnpe,1:maxelm) = 0 npoint = 0 point(1:maxpoi) = 0.0E+00 title = '-d/dx K(x) dU/dx + C(x) dU/dx + B(x)*U = F(x) + P(x)' vbc(1:2,1:2) = 0.0E+00 wgauss(1:maxgauss) = 0.0E+00 x(1:maxnod) = 0.0E+00 xgauss(1:maxgauss) = 0.0E+00 return end function lcm ( i, j ) ! !******************************************************************************* ! !! LCM computes the least common multiple of two integers. ! ! ! Definition: ! ! The LCM may be defined as ! ! LCM(I,J) = ABS( I * J ) / GCF(I,J) ! ! where GCF(I,J) is the greatest common factor of I and J. ! ! Modified: ! ! 01 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, J, the integers whose LCM is desired. ! ! Output, integer LCM, the least common multiple of I and J. ! LCM is never negative. LCM is 0 if either I or J is zero. ! implicit none ! integer i integer i_gcd integer j integer lcm ! lcm = abs ( i * ( j / i_gcd ( i, j ) ) ) return end subroutine material_evaluate ( infix, irpn, mater, maxfix, maxrpn, xb, xc, & xf, xk, xx ) ! !******************************************************************************* ! !! MATERIAL_EVALUATE evaluates the material function formulas. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MATER, the index of the material. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Output, real XB, XC, XF, XK, the value of the material functions ! B(X), C(X), F(X) and K(X). ! ! Input, real XX, the (physical) point at which the material functions ! are to be evaluated. ! implicit none ! integer maxfix integer maxrpn ! character com integer ierror integer ifrm character ( len = * ) infix integer irpn(maxrpn) integer mater character ( len = 10 ) namvar real xb real xc real xf real xk real xx real value ! ! Set the value of X. ! com = 'V' ifrm = 1 namvar = 'X' value = xx call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) ! ! Evaluate K(X). ! ifrm = (mater-1)*4+1 com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) xk = value ! ! Evaluate C(X). ! ifrm = (mater-1)*4+2 com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) xc = value ! ! Evaluate B(X). ! ifrm = (mater-1)*4+3 com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) xb = value ! ! Evaluate F(X). ! ifrm = (mater-1)*4+4 com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) xf = value return end subroutine material_read ( ierror, infix, inform, irpn, maxfix, maxmat, & maxrpn, nmat ) ! !******************************************************************************* ! !! MATERIAL_READ reads the material function formulas. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXMAT, the maximum number of materials. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Input, integer NMAT, the number of materials. ! implicit none ! integer maxfix integer maxrpn ! character com character ctemp integer i integer ierror integer ifrm integer ihi integer ilo character ( len = * ) infix character ( len = * ) inform integer irpn(maxrpn) integer j integer maxmat character ( len = 10 ) namvar integer nmat real value ! namvar = ' ' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Set up the material functions.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Set the functions K(x), C(x), B(x) and F(x) in' write ( *, '(a)' ) ' -d/dx K dU/dx + C dU/dx + B U = F + P' ! ! Get the number of materials. ! 10 continue write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Enter the number of materials, up to ', maxmat write ( *, '(a)' ) '(Simple problems use 1 material.)' read ( *, * ) nmat if ( nmat < 1 .or. nmat > maxmat ) then write ( *, '(a)' ) 'Your input was not acceptable!' go to 10 end if ! ! Get formulas for each material. ! ifrm = 0 do i = 1, nmat write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'For material ', i write ( *, '(a)' ) 'enter formulas for K(x), C(x), B(x) and F(x):' write ( *, '(a)' ) ' ' do j = 1, 4 if ( j == 1 ) then ctemp = 'K' else if ( j == 2 ) then ctemp = 'C' else if ( j == 3 ) then ctemp = 'B' else if ( j == 4 ) then ctemp = 'F' end if ifrm = ifrm+1 20 continue write ( *, '(a)' ) 'Enter the formula for the '// ctemp // ' function.' read ( *, '(a)' ) infix ilo = (ifrm-1)*80+1 ihi = (ifrm-1)*80+80 inform(ilo:ihi) = infix com = 'F' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) if ( ierror /= 0 ) then write ( *, '(a)' ) 'Your formula did not compile!' ierror = 0 go to 20 end if end do end do return end subroutine material_write ( inform, nmat ) ! !******************************************************************************* ! !! MATERIAL_WRITE prints the material function formulas. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NMAT, the number of materials. ! implicit none ! integer i integer ifrm integer ihi integer ilo character ( len = * ) inform integer j character label(4) integer nmat ! label(1) = 'K' label(2) = 'C' label(3) = 'B' label(4) = 'F' ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Material functions.' write ( *, '(a)' ) ' ' if ( nmat==1 ) then write ( *, '(a)' ) 'There is 1 material.' else write ( *, '(a)' ) 'There are ', nmat, ' separate materials.' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Each material is defined by functions' write ( *, '(a)' ) ' K(x), C(x), B(x) and F(x), which appear in' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' - d/dx K(x) dU/dx + C(x) dU/dx + B(x) U = F(x) + P' ifrm = 0 do i = 1, nmat write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Functions for material ',i write ( *, '(a)' ) ' ' do j = 1, 4 ifrm = ifrm+1 ilo = (ifrm-1)*80+1 ihi = (ifrm-1)*80+80 write ( *, '(a,a,a)' ) label(j), '(x) = ', inform(ilo:ihi) end do end do return end subroutine matwrt ( ierror, namvr, ncol, nrow, value ) ! !******************************************************************************* ! !! MATWRT prints out a scalar, vector or array. ! ! ! Discussion: ! ! The items to be printed are stored in VALUE, an array of logical ! first dimension MAXROW and used dimensions NROW by NCOL. ! ! It tries to print out integer valued reals in compact format. ! It transposes column vectors. It can print up to 5 columns ! of an array or vector per line. ! ! For an array, it will label the columns. ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred. ! ! Input, integer MAXROW, the leading dimension of the array holding the data. ! ! Input, CHARACTER*10 NAMVR, the name of the item to be printed. ! ! Input, integer NCOL, NROW, the number of columns and rows of data. ! ! Input, real VALUE(MAXROW,NCOL), the NROW by NCOL array of data. ! implicit none ! integer, parameter :: maxlen = 10 ! ! ! You don't want to know why I had to use two billion! ! ! real, parameter :: big = 2147483647.0E+00 real, parameter :: big = 2000000000.0E+00 ! integer ncol ! character ( len = 14 ) chrtmp integer i integer i1 integer i2 integer ierror integer ihi integer ii integer ilo integer iv integer j integer jhi integer jlo character ( len = maxlen ) namvr integer ndim integer nrow character ( len = 100 ) output real v real value ! ! Print a scalar value. ! v = value write ( chrtmp, '(g14.6)' ) v output = namvr//' '//chrtmp if ( abs(v)<=big ) then iv = int(v) if ( real(iv) == v ) then write ( output,'(a10,1x,i14)')namvr,iv end if end if call s_blanks_delete ( output ) write ( *, '(a)' ) trim ( output ) return end subroutine node_read ( kind, maxelm, maxnod, maxnpe, nelem, nnode, nodes, x ) ! !******************************************************************************* ! !! NODE_READ sets the X coordinates of the nodes. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Output, real X(NNODE), the (physical) X coordinates of the nodes. ! implicit none ! integer maxelm integer maxnod integer maxnpe ! integer i integer ilelem integer irelem integer kind(maxelm) integer kindp1 integer nelem integer nnode integer nodel integer noder integer nodes(maxnpe,maxelm) real x(maxnod) real x1 real x2 ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Set up the nodes.' write ( *, '(a)' ) ' ' ! ! Get the X coordinate of the left endpoint of the first element. ! ilelem = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter the location of the LEFT boundary node.' read ( *, * ) x1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Define equally spaced subintervals by giving' write ( *, '(a)' ) 'an element, and its rightmost X coordinate.' x(1) = x1 10 continue write ( *, '(a,i6)' ) 'The left element is ', ilelem write ( *, '(a,g14.6)' ) 'The left X value is ', x1 ! ! Request the right element, and its rightmost X coordinate. ! 20 continue write ( *, '(a,i6,a,i6)' ) 'Enter right element between ', ilelem, & ' and ', nelem read ( *, * ) irelem if ( irelem < ilelem ) then write ( *, '(a)' ) 'Warning, your input was not acceptable!' write ( *, '(a)' ) 'The right element number must be GREATER than' write ( *, '(a)' ) 'the left one.' go to 20 end if if ( irelem > nelem ) then write ( *, '(a)' ) 'Warning, your input was not acceptable!' write ( *, '(a)' ) 'Your element number was larger than NELEM.' go to 20 end if ! ! Request rightmost x coordinate of right element ! do write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) & 'Enter X coordinate of rightmost node in element ', irelem read ( *, * ) x2 if ( x2 > x1 ) then exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'No, the X position of the right element must be' write ( *, '(a)' ) 'GREATER than that of the left element.' end do ! ! Assign X coordinates for all nodes between left and right endpoints. ! nodel = nodes(1,ilelem) kindp1 = kind(irelem)+1 noder = nodes(kindp1,irelem) do i = nodel, noder x(i) = ( real ( noder - i ) * x1 + real ( i - nodel ) * x2 ) & / real ( noder - nodel ) end do ! ! Prepare for the next segment. ! if ( noder < nnode ) then nodel = noder ilelem = irelem + 1 x1 = x2 go to 10 end if return end subroutine node_write ( nnode, x ) ! !******************************************************************************* ! !! NODE_WRITE prints out the location of the nodes. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer NNODE, the number of nodes. ! ! Input, real X(NNODE), the locations of the nodes. ! implicit none ! integer nnode ! integer i real x(nnode) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Node X coordinate' write ( *, '(a)' ) ' ' do i = 1, nnode write ( *, '(i4,g14.6)' ) i, x(i) end do return end subroutine point_read ( kpoint, maxpoi, nnode, npoint, point ) ! !******************************************************************************* ! !! POINT_READ gets the point load data. ! ! ! Modified: ! ! 16 May 2001 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer MAXPOI, the maximum number of point loads. ! ! Input, integer NNODE, the number of nodes. ! ! Output, integer NPOINT, the number of point loads. ! implicit none ! integer maxpoi ! integer i integer kpoint(maxpoi) integer nnode integer npoint integer nval real point(maxpoi) real rval ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Set up the point loads.' write ( *, '(a)' ) '(Simple problems have 0 point loads).' npoint = 0 do i = 1, maxpoi do write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Enter a node number, between 1 and ', nnode write ( *, '(a)' ) 'where a point load applies,' write ( *, '(a)' ) 'and the value of the point load there.' write ( *, '(a)' ) '(0, 0.0E+00 to quit)' read ( *, * ) nval, rval if ( nval == 0 ) then return end if if ( 1 <= nval .and. nval <= nnode ) then exit end if write ( *, '(a)' ) 'That node number was out of bounds!' end do npoint = npoint + 1 kpoint(i) = nval point(i) = rval end do return end subroutine point_write ( kpoint, maxpoi, npoint, point ) ! !******************************************************************************* ! !! POINT_WRITE prints out the point load data. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Input, integer MAXPOI, the maximum number of point loads. ! ! Input, integer NPOINT, the number of point loads. ! implicit none ! integer maxpoi ! integer i integer kpoint(maxpoi) integer npoint real point(maxpoi) ! write ( *, '(a)' ) ' ' if ( npoint <= 0 ) then write ( *, '(a)' ) 'There are no point loads!' return end if write ( *, '(a)' ) 'Point loads:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Load at Node Value' write ( *, '(a)' ) ' ' do i = 1, npoint write ( *, '(i4,i8,g14.6)' ) i, kpoint(i), point(i) end do return end subroutine prep ( ierror, infix, inform, irpn, kbc, kind, kpoint, mat, & maxelm, maxfix, maxgauss, maxmat, maxnod, maxnpe, maxpoi, maxrpn, nelem, & ngauss, nmat, nnode, nodes, npoint, point, title, vbc, x ) ! !******************************************************************************* ! !! PREP does the preprocessing of the problem. ! ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error. ! nonzero, an error occurred, and the problem was not set up. ! ! Input, integer KBC(2). KBC(1) and KBC(2) are the boundary condition ! codes at the left and right endpoints of the region. ! 1: U = value1; ! 2: K dU/dx = value1 * U - value2 ! ! Output, integer KIND(MAXELM), the degree of each element. ! ! Output, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXPOI, the maximum number of point loads. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Output, integer NELEM, the number of elements. ! ! Output, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Output, integer NMAT, the number of materials. ! ! Output, integer NNODE, the number of nodes. ! ! Output, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Output, integer NPOINT, the number of point loads. ! ! Output, character ( len = * ) TITLE, the title of the problem. ! ! Output, real VBC(2,2), contains constants for boundary conditions. ! VBC(1,1) and VBC(2,1) are value1 and value2 for the left endpoint; ! VBC(1,2) and VBC(2,2) are value1 and value2 for the right endpoint. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxnpe integer maxpoi integer maxrpn ! integer ierror character ( len = * ) infix character ( len = * ) inform integer irpn(maxrpn) integer kbc(2) integer kind(maxelm) integer kpoint(maxpoi) integer mat(maxelm) integer maxmat integer nelem integer ngauss integer nmat integer nnode integer nodes(maxnpe,maxelm) integer npoint real point(maxpoi) character ( len = * ) title real vbc(2,2) real x(maxnod) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter a title for the problem' read ( *, '(a)' ) title ! ! Get the material properties K(x), C(x), B(x) and F(x). ! call material_read ( ierror, infix, inform, irpn, maxfix, maxmat, maxrpn, & nmat ) if ( ierror /= 0 ) then return end if ! ! Set up the elements. ! call element_read ( ierror, kind, mat, maxelm, maxnod, maxnpe, nelem, nmat, & nnode, nodes ) if ( ierror /= 0 ) then return end if ! ! Set the quadrature order. ! call gauss_read ( ierror, maxgauss, ngauss ) if ( ierror /= 0 ) then return end if ! ! Get the X coordinates of the nodes. ! call node_read ( kind, maxelm, maxnod, maxnpe, nelem, nnode, nodes, x ) ! ! Get the boundary conditions. ! call bc_read ( kbc, vbc ) ! ! Get the point loads. ! call point_read ( kpoint, maxpoi, nnode, npoint, point ) return end subroutine proc ( dpsi, ef, ek, gf, gk, gklu, gu, ierror, infix, ipivot, irpn, & kbc, kind, kpoint, ldgk, mat, maxelm, maxfix, maxgauss, maxnod, maxnpe, & maxpoi, maxrpn, ml, mu, nelem, ngauss, nnode, nodes, npoint, nrow, point, & psi, vbc, wgauss, x, xgauss ) ! !******************************************************************************* ! !! PROC assembles the linear system and solves it. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Workspace, real DPSI(MAXNPE), shape functio derivatives. ! ! Workspace, real EF(MAXNPE), the element force vector. ! ! Workspace, real EK(MAXNPE,MAXNPE), the element stiffness matrix. ! ! Workspace, real GF(MAXNOD), the global force vector. ! ! Workspace, real GK(NROW,NNODE), the global stiffness matrix. ! ! Input, integer KBC(2). KBC(1) and KBC(2) are the boundary condition ! codes at the left and right endpoints of the region. ! 1: U = value1; ! 2: K dU/dx = value1 * U - value2 ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAT(MAXELM), the material index of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXGAUSS, the maximum quadrature rule order. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXPOI, the maximum number of point loads. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Output, integer ML, MU, the lower and upper bandwidths of the ! global stiffness matrix. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NGAUSS, the order of the Gauss quadrature rule. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Input, integer NPOINT, the number of point loads. ! ! Output, integer NROW, the leading dimension of the global stiffness matrix. ! ! Workspace, real PSI(MAXNPE), shape function values. ! ! Input, real VBC(2,2), contains constants for boundary conditions. ! VBC(1,1) and VBC(2,1) are value1 and value2 for the left endpoint; ! VBC(1,2) and VBC(2,2) are value1 and value2 for the right endpoint. ! ! Workspace, real WGAUSS(MAXGAUSS), quadrature weights. ! ! Input, real X(MAXNOD), the X coordinates of the nodes. ! ! Workspace, real XGAUSS(MAXGAUSS), quadrature abscissas. ! implicit none ! integer ldgk integer maxelm integer maxfix integer maxgauss integer maxnod integer maxnpe integer maxpoi integer maxrpn ! real dpsi(maxnpe) real ef(maxnpe) real ek(maxnpe,maxnpe) real gf(maxnod) real gk(ldgk,maxnod) real gklu(ldgk,maxnod) real gu(maxnod) integer i integer ierror character ( len = * ) infix integer info integer inode integer ipivot(maxnod) integer irpn(maxrpn) integer j integer job integer kbc(2) integer kind(maxelm) integer kpoint(maxpoi) integer mat(maxelm) integer ml integer mu integer nelem integer ngauss integer nnode integer nodes(maxnpe,maxelm) integer npoint integer nrow real point(maxpoi) real psi(maxnpe) real vbc(2,2) real wgauss(maxgauss) real x(maxnod) real xgauss(maxgauss) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Set up and solve the linear system:' ! ! Set the Gauss weights and abscissas. ! call gauss_set ( ngauss, wgauss, xgauss ) ! ! Compute the bandwidth of the matrix. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) '1: Compute the bandwidth of the stiffness matrix GK.' call gkf_width ( kind, maxelm, maxnpe, ml, mu, nelem, nodes ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Lower bandwidth = ', ml write ( *, '(a,i6)' ) ' Upper bandwidth = ', mu ! ! Set up storage information to be used later. ! nrow = ml+ml+mu+1 ! ! Set the global stiffness matrix GK and force vector GF. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) '2: Assemble the stiffness matrix GK and right hand side GF.' call gkf_assemble ( dpsi, ef, ek, gf, gk, infix, irpn, kind, mat, & maxelm, maxfix, maxgauss, maxnod, maxnpe, maxrpn, ml, mu, & nelem, ngauss, nnode, nodes, nrow, psi, wgauss, x, xgauss ) ! ! Modify the global stiffness matrix GK and the right hand side GF ! to account for boundary conditions. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) '3: Modify the system for boundary conditions.' call bc_apply ( gf, gk, ierror, kbc, maxnod, ml, mu, nnode, nrow, vbc ) if ( ierror /= 0 ) then return end if ! ! Modify the force vector GF for point loads. ! write ( *, '(a)' ) '4: Modify the right hand side GF for point loads.' do i = 1, npoint inode = kpoint(i) gf(inode) = gf(inode) + point(i) end do ! ! Factor the matrix K in the linear system GK * U = GF. ! write ( *, '(a)' ) '5: Factor the stiffness matrix.' gklu(1:ldgk,1:nnode) = gk(1:ldgk,1:nnode) call sgb_fa ( gklu, nrow, nnode, ml, mu, ipivot, info ) if ( info /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PROC - Fatal error!' write ( *, '(a)' ) ' The stiffness matrix is singular.' write ( *, '(a,i6)' ) ' INFO = ', info ierror = 1 return end if ! ! Solve the linear system GK * U = GF. ! write ( *, '(a)' ) '6: Solve the system GK * U = GF.' gu(1:nnode) = gf(1:nnode) job = 0 call sgb_sl ( gklu, nrow, nnode, ml, mu, ipivot, gu, job ) return end subroutine rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) ! !******************************************************************************* ! !! RPNCHK examines an RPN formula, looking for a complete RPN expression. ! ! ! Discussion: ! ! The routine starts at location IHI, and finds the position ILO ! such that IRPN(ILO)...IRPN(IHI) represents a single argument. ! ! Modified: ! ! 24 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, 0, no error; 1 an error. ! ! Input, integer IHI, the location in IRPN where the search begins. ! ! Output, integer ILO, the location in IRPN such that IRPN(ILO) through ! IRPN(IHI) represents a single argument, or IHI+1 if no such location ! was found. ! ! Input, integer IOPSYM(MAXSYM), contains, for each symbol, the number ! of operands. In particular, if a symbol represents a constant, ! IOPSYM(I) is 0. If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! implicit none ! integer maxrpn integer maxsym ! integer ierror integer ihi integer ilo integer iopsym(maxsym) integer irpn(maxrpn) integer isum ! isum = 0 ierror = 0 ilo = ihi + 1 do ilo = ilo - 1 if ( ilo <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNCHK - Fatal error!' write ( *, '(a)' ) ' Reached beginning of RPN formula,' write ( *, '(a)' ) ' but still had operators missing arguments.' ierror = 1 return end if if ( iopsym(irpn(ilo)) < 0 .and. isum == 0 ) then cycle end if isum = isum + 1 - iopsym(irpn(ilo)) if ( isum == 1 ) then exit end if end do return end subroutine rpnset ( ierror, ifinis, imin, intsym, iopsym, iprsym, irpn, & istack, maxrpn, maxsym, nints, nrpn, symbol ) ! !******************************************************************************* ! !! RPNSET converts the infix formula into an RPN formula. ! ! ! Modified: ! ! 25 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Input, integer IMIN, an offset for accessing IRPN. ! ! Input, integer INTSYM(80), a set of integers which ! are the indices of symbols, representing an infix formula. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the priority of each symbol. ! ! Output, integer IRPN(MAXRPN), the RPN version of the infix formula. ! ! Workspace, integer ISTACK(MAXSYM). ! ! Input, integer MAXRPN, the maximum number of RPN symbols allowed. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer NINTS, the number of useful entries in INTSYM. ! ! Output, integer NRPN, the number of useful entries in IRPN. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! implicit none ! integer, parameter :: maxlen = 20 ! integer maxrpn integer maxsym ! integer ierror integer ifinis integer imin integer intsym(80) integer iopsym(maxsym) integer iprsym(maxsym) integer iread integer irpn(maxrpn) integer istack(maxsym) integer istak integer isym integer jsym integer nints integer nrpn character ( len = maxlen ) sym character ( len = maxlen ) sym2 character ( len = maxlen ) symbol(maxsym) ! ierror = 0 nrpn = 0 ! ! An error if the formula seems to have nothing in it. ! if ( nints <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNSET - Error!' write ( *, '(a)' ) ' This formula seems to be blank!' ierror = 1 return end if iread = 0 istak = 0 ! ! Read symbol number IREAD from INTSYM. ! do iread = iread + 1 if ( iread > nints ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNSET - Error!' write ( *, '(a)' ) ' This formula does not make sense!' return end if isym = intsym(iread) if ( isym == 0 ) then cycle end if sym = symbol(isym) if ( sym == ',' ) then cycle end if ! ! If the symbol is "$", then it's time to pop the stack. ! if ( sym == '$' ) then do if ( istak <= 0 ) then nrpn = nrpn + 1 irpn(nrpn+imin) = ifinis if ( iread < nints ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNSET - Error!' write ( *, '(a)' ) ' Some of the formula is left over!' ierror = 1 end if return end if nrpn = nrpn + 1 irpn(nrpn+imin) = istack(istak) istak = istak - 1 end do end if ! ! A left parenthesis goes onto the stack. ! if ( sym == '(' ) then istak = istak + 1 istack(istak) = isym cycle end if ! ! A variable or constant goes immediately into IRPN. ! if ( iopsym(isym) == 0 ) then nrpn = nrpn + 1 irpn(nrpn+imin) = isym cycle end if ! ! Put operators onto the stack. ! do if ( istak <= 0 ) then istak = istak + 1 istack(istak) = isym exit end if jsym = istack(istak) if ( iprsym(jsym) <= iprsym(isym) ) then jsym = istack(istak) sym2 = symbol(jsym) if ( sym == ')' .and. sym2 == '(' ) then istak = istak - 1 else istak = istak + 1 istack(istak) = isym end if exit end if nrpn = nrpn + 1 irpn(nrpn+imin) = istack(istak) istak = istak - 1 end do end do return end subroutine rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, & irad, irpn, istack, maxrpn, maxsym, maxval, & nrpn, nsym, nsyms, symbol, valsym, value ) ! !******************************************************************************* ! !! RPNVAL evaluates the symbolic functions in an RPN formula. ! ! ! Discussion: ! ! RPNVAL determines the number of arguments, gathering their values, ! and calling FUNVAL to evaluate the given function. ! ! Modified: ! ! 01 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real EPS, the value of the machine precision. ! ! Output, integer IERROR. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFREE, the index of the next free location in VALSYM. ! ! Input, integer IMIN, an offset for indexing IRPN. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Workspace, integer IRPN(MAXRPN), used to store the compiled ! versions of user formulas. ! ! Workspace, integer ISTACK(MAXSYM), workspace for interpreting ! the formula. ! ! Input, integer MAXRPN, specifies the length of IRPN. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed in VALSYM. ! ! Output, integer NCOL, the number of columns in the result. ! ! Output, integer NROW, the number of rows in the result. ! ! Input, integer NRPN, the number of useful entries in IRPN. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMS, the number of declared symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = * ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! ! Workspace, real VALUE, space used to hold ! the value of the variable to be saved. ! implicit none ! integer, parameter :: maxlen = 20 ! integer maxrpn integer maxsym ! character ( len = 3 ) ctemp integer iarg1 integer iarg2 integer ierror integer ifree integer imin integer index1 integer index4 integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer iread integer irpn(maxrpn) integer istack(maxsym) integer istak integer isym integer isymo integer itemp integer maxval integer nrpn integer nsym integer nsyms character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value ! value = 0.0E+00 if ( nrpn <= 1 ) then return end if istak = 0 isymo = 0 nsyms = nsym itemp = 0 iread=0 do iread = iread + 1 if ( iread > nrpn ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNVAL - Error!' write ( *, '(a)' ) ' The formula ended unexpectedly!' exit end if isym = irpn(iread+imin) if ( 0 < isym .and. isym <= maxsym ) then sym = symbol(isym) else sym = '$' end if if ( sym == '$' .or. iread > nrpn ) then index4 = ipval(isymo) value = valsym(index4) exit end if ! ! Constants and variables go into stack. ! if ( iopsym(isym) == 0 ) then if ( isym <= nsym ) then istak = istak+1 istack(istak) = isym isymo = isym else if ( nsyms >= maxsym ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNVAL - Error!' write ( *, '(a)' ) ' There is no memory left for more symbols!' write ( *, '(a)' ) ' Perhaps the KILL or INIT command would help.' ierror = 1 exit end if nsyms = nsyms+1 isymo = nsyms istak = istak+1 istack(istak) = nsyms iprsym(nsyms) = 10 iopsym(nsyms) = 0 itemp = itemp+1 ctemp = ' ' call i_to_s_zero ( itemp, ctemp ) symbol(nsyms) = 'STK000' symbol(nsyms)(4:6) = ctemp ipval(nsyms) = ifree ifree = ifree + 1 index1 = ipval(isym) index4 = ipval(nsyms) valsym(index4) = valsym(index1) end if cycle end if ! ! Pull off arguments. ! iarg1 = istack(istak) iarg2 = 0 if ( iopsym(isym) == 2 ) then iarg2 = istack(istak) istak = istak-1 iarg1 = istack(istak) else if (iopsym(isym) == 3 ) then istak = istak-1 iarg2 = istack(istak) istak = istak-1 iarg1 = istack(istak) end if sym = symbol(isym) call funval ( iarg1, iarg2, ierror, ifree, iopsym, iprsym, ipset, & ipval, irad, itemp, maxsym, maxval, nsyms, sym, symbol, valsym ) isymo = nsyms if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'RPNVAL - Warning!' write ( *, '(a)' ) ' Evaluation of this formula is abandoned.' value = 0.0E+00 exit end if if ( ipset == 0 ) then ipval(nsyms) = ifree ifree = ifree + 1 else ipval(nsyms) = ipset ipset = 0 end if istack(istak) = nsyms end do return end subroutine s_blank_delete ( s ) ! !******************************************************************************* ! !! S_BLANK_DELETE removes blanks from a string, left justifying the remainder. ! ! ! Comment: ! ! All TAB characters are also removed. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none ! character c integer iget integer iput character ( len = * ) s character, parameter :: TAB = char ( 9 ) ! iput = 0 do iget = 1, len ( s ) c = s(iget:iget) if ( c /= ' ' .and. c /= TAB ) then iput = iput + 1 s(iput:iput) = c end if end do s(iput+1:) = ' ' return end subroutine s_blanks_delete ( s ) ! !******************************************************************************* ! !! S_BLANKS_DELETE replaces consecutive blanks by one blank. ! ! ! Discussion: ! ! The remaining characters are left justified and right padded with blanks. ! TAB characters are converted to spaces. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none ! integer i integer j character newchr character oldchr character ( len = * ) s character, parameter :: TAB = char ( 9 ) ! j = 0 newchr = ' ' do i = 1, len ( s ) oldchr = newchr newchr = s(i:i) if ( newchr == TAB ) then newchr = ' ' end if s(i:i) = ' ' if ( oldchr /= ' ' .or. newchr /= ' ' ) then j = j + 1 s(j:j) = newchr end if end do return end subroutine s_cap ( string ) ! !******************************************************************************* ! !! S_CAP replaces any lowercase letters by uppercase ones in a string. ! ! ! Modified: ! ! 16 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) STRING, the string to be transformed. ! implicit none ! character c integer i integer nchar character ( len = * ) string ! nchar = len ( string ) do i = 1, nchar c = string(i:i) call ch_cap ( c ) string(i:i) = c end do return end function s_eqi ( strng1, strng2 ) ! !******************************************************************************* ! !! S_EQI is a case insensitive comparison of two strings for equality. ! ! ! Examples: ! ! S_EQI ( 'Anjana', 'ANJANA' ) is .TRUE. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) STRNG1, STRNG2, the strings to compare. ! ! Output, logical S_EQI, the result of the comparison. ! implicit none ! integer i integer len1 integer len2 integer lenc logical s_eqi character s1 character s2 character ( len = * ) strng1 character ( len = * ) strng2 ! len1 = len ( strng1 ) len2 = len ( strng2 ) lenc = min ( len1, len2 ) s_eqi = .false. do i = 1, lenc s1 = strng1(i:i) s2 = strng2(i:i) call ch_cap ( s1 ) call ch_cap ( s2 ) if ( s1 /= s2 ) then return end if end do do i = lenc + 1, len1 if ( strng1(i:i) /= ' ' ) then return end if end do do i = lenc + 1, len2 if ( strng2(i:i) /= ' ' ) then return end if end do s_eqi = .true. return end function s_indexi ( s, sub ) ! !******************************************************************************* ! !! S_INDEXI is a case-insensitive INDEX function. ! ! ! Discussion: ! ! The function returns the location in the string at which the ! substring SUB is first found, or 0 if the substring does not ! occur at all. ! ! The routine is also trailing blank insensitive. This is very ! important for those cases where you have stored information in ! larger variables. If S is of length 80, and SUB is of ! length 80, then if S = 'FRED' and SUB = 'RED', a match would ! not be reported by the standard FORTRAN INDEX, because it treats ! both variables as being 80 characters long! This routine assumes that ! trailing blanks represent garbage! ! ! Because of the suppression of trailing blanks, this routine cannot be ! used to find, say, the first occurrence of the two-character ! string 'A '. However, this routine treats as a special case the ! occurrence where S or SUB is entirely blank. Thus you can ! use this routine to search for occurrences of double or triple blanks ! in a string, for example, although INDEX itself would be just as ! suitable for that problem. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be searched. ! ! Input, character ( len = * ) SUB, the substring to search for. ! ! Output, integer S_INDEXI. 0 if SUB does not occur in ! the string. Otherwise S(S_INDEXI:S_INDEXI+LENS-1) = SUB, ! where LENS is the length of SUB, and is the first place ! this happens. However, note that this routine ignores case, ! unlike the standard FORTRAN INDEX function. ! implicit none ! integer i integer llen1 integer llen2 character ( len = * ) s logical s_eqi integer s_indexi character ( len = * ) sub ! s_indexi = 0 llen1 = len_trim ( s ) llen2 = len_trim ( sub ) ! ! In case S or SUB is blanks, use LEN. ! if ( llen1 == 0 ) then llen1 = len ( s ) end if if ( llen2 == 0 ) then llen2 = len ( sub ) end if if ( llen2 > llen1 ) then return end if do i = 1, llen1 + 1 - llen2 if ( s_eqi ( s(i:i+llen2-1), sub ) ) then s_indexi = i return end if end do return end function s_is_alpha ( s ) ! !******************************************************************************* ! !! S_IS_ALPHA returns .TRUE. if the string contains only alphabetic characters. ! ! ! Discussion: ! ! Here, alphabetic characters are 'A' through 'Z' and 'a' through 'z'. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be checked. ! ! Output, logical S_IS_ALPHA, .TRUE. if the string contains only ! alphabetic characters, .FALSE. otherwise. ! implicit none ! logical ch_is_alpha integer i character ( len = * ) s logical s_is_alpha ! s_is_alpha = .false. do i = 1, len ( s ) if ( .not. ch_is_alpha ( s(i:i) ) ) then return end if end do s_is_alpha = .true. return end function s_paren_check ( s ) ! !******************************************************************************* ! !! S_PAREN_CHECK checks the parentheses in a string. ! ! ! Discussion: ! ! Blanks are removed from the string, and then the following checks ! are made: ! ! 1) as we read the string left to right, there must never be more ! right parentheses than left ones; ! 2) there must be an equal number of left and right parentheses; ! 3) there must be no occurrences of the abutting packages '...)(...'. ! 4) there must be no occurrences of the empty package '()'. ! ! Modified: ! ! 20 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to check. ! ! Output, logical S_PAREN_CHECK is TRUE if the string passed the checks. ! implicit none ! integer i integer isum character ( len = * ) s character ( len = 256 ) s_copy integer s_len logical s_paren_check ! s_copy = s call s_blank_delete ( s_copy) s_len = len_trim ( s_copy ) ! ! 1) Letting '(' = +1 and ')' = -1, check that the running parentheses sum ! is always nonnegative. ! isum = 0 do i = 1, s_len if ( s_copy(i:i) == '(' ) then isum = isum + 1 end if if ( s_copy(i:i) == ')' ) then isum = isum - 1 if ( isum < 0 ) then s_paren_check = .false. return end if end if end do ! ! 2) Check that the final parentheses sum is zero. ! if ( isum /= 0 ) then s_paren_check = .false. return end if ! ! 3) Check that there are no ")(" pairs. ! do i = 2, s_len if ( s_copy(i-1:i) == ')(' ) then s_paren_check = .false. return end if end do ! ! 4) Check that there are no "()" pairs. ! do i = 2, s_len if ( s_copy(i-1:i) == '()' ) then s_paren_check = .false. return end if end do ! ! The checks were passed. ! s_paren_check = .true. return end subroutine s_to_r ( s, r, ierror, lchar ) ! !******************************************************************************* ! !! S_TO_R reads a real number from a string. ! ! ! Discussion: ! ! This routine will read as many characters as possible until it reaches ! the end of the string, or encounters a character which cannot be ! part of the real number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 spaces ! 3 integer part, ! 4 decimal point, ! 5 fraction part, ! 6 'E' or 'e' or 'D' or 'd', exponent marker, ! 7 exponent sign, ! 8 exponent integer part, ! 9 exponent decimal point, ! 10 exponent fraction part, ! 11 blanks, ! 12 final comma or semicolon. ! ! with most quantities optional. ! ! Examples: ! ! S R ! ! '1' 1.0 ! ' 1 ' 1.0 ! '1A' 1.0 ! '12,34,56' 12.0 ! ' 34 7' 34.0 ! '-1E2ABCD' -100.0 ! '-1X2ABCD' -1.0 ! ' 2E-1' 0.2 ! '23.45' 23.45 ! '-4.2E+2' -420.0 ! '17d2' 1700.0 ! '-14e-2' -0.14 ! 'e2' 100.0 ! '-12.73e-9.23' -12.73 * 10.0**(-9.23) ! ! Modified: ! ! 12 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate at the end of the string, or when no more ! characters can be read to form a legal real. Blanks, ! commas, or other nonnumeric data will, in particular, ! cause the conversion to halt. ! ! Output, real R, the real value that was read from the string. ! ! Output, integer IERROR, error flag. ! ! 0, no errors occurred. ! ! 1, 2, 6 or 7, the input number was garbled. The ! value of IERROR is the last type of input successfully ! read. For instance, 1 means initial blanks, 2 means ! a plus or minus sign, and so on. ! ! Output, integer LCHAR, the number of characters read from ! the string to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none ! logical ch_eqi character c integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer lchar integer nchar integer ndig real r real rbot real rexp real rtop character ( len = * ) s character, parameter :: TAB = char ( 9 ) ! nchar = len_trim ( s ) ierror = 0 r = 0.0E+00 lchar = - 1 isgn = 1 rtop = 0.0E+00 rbot = 1.0E+00 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do lchar = lchar + 1 c = s(lchar+1:lchar+1) ! ! Blank or TAB character. ! if ( c == ' ' .or. c == TAB ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( ihave > 1 ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 lchar = lchar + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 isgn = - 1 else if ( ihave == 6 ) then ihave = 7 jsgn = - 1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 else if ( ihave == 6 ) then ihave = 7 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else if ( ihave >= 6 .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! Exponent marker. ! else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then if ( ihave < 6 ) then ihave = 6 else iterm = 1 end if ! ! Digit. ! else if ( ihave < 11 .and. lge ( c, '0' ) .and. lle ( c, '9' ) ) then if ( ihave <= 2 ) then ihave = 3 else if ( ihave == 4 ) then ihave = 5 else if ( ihave == 6 .or. ihave == 7 ) then ihave = 8 else if ( ihave == 9 ) then ihave = 10 end if call ch_to_digit ( c, ndig ) if ( ihave == 3 ) then rtop = 10.0E+00 * rtop + real ( ndig ) else if ( ihave == 5 ) then rtop = 10.0E+00 * rtop + real ( ndig ) rbot = 10.0E+00 * rbot else if ( ihave == 8 ) then jtop = 10 * jtop + ndig else if ( ihave == 10 ) then jtop = 10 * jtop + ndig jbot = 10 * jbot end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if ! ! If we haven't seen a terminator, and we haven't examined the ! entire string, go get the next character. ! if ( iterm == 1 .or. lchar+1 >= nchar ) then exit end if end do ! ! If we haven't seen a terminator, and we have examined the ! entire string, then we're done, and LCHAR is equal to NCHAR. ! if ( iterm /= 1 .and. lchar+1 == nchar ) then lchar = nchar end if ! ! Number seems to have terminated. Have we got a legal number? ! Not if we terminated in states 1, 2, 6 or 7! ! if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then ierror = ihave return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0E+00 else if ( jbot == 1 ) then rexp = 10.0E+00**( jsgn * jtop ) else rexp = jsgn * jtop rexp = rexp / jbot rexp = 10.0E+00**rexp end if end if r = isgn * rexp * rtop / rbot return end subroutine sgb_check ( lda, m, n, ml, mu, ierror ) ! !******************************************************************************* ! !! SGB_CHECK checks the dimensions of a general band matrix. ! ! ! Modified: ! ! 18 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least 2 * ML + MU + 1. ! ! Input, integer M, the number of rows of the matrix. ! M must be positive. ! ! Input, integer N, the number of columns of the matrix. ! N must be positive. ! ! Input, integer ML, MU, the lower and upper bandwidths. ! ML and MU must be nonnegative, and no greater than min(M,N)-1. ! ! Output, integer IERROR, reports whether any errors were detected. ! IERROR is set to 0 before the checks are made, and then: ! IERROR = IERROR + 1 if LDA is illegal; ! IERROR = IERROR + 2 if M is illegal; ! IERROR = IERROR + 4 if ML is illegal; ! IERROR = IERROR + 8 if MU is illegal; ! IERROR = IERROR + 16 if N is illegal. ! implicit none ! integer ierror integer lda integer m integer ml integer mu integer n ! ierror = 0 if ( lda < 2 * ml + mu + 1 ) then ierror = ierror + 1 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGB_CHECK - Illegal LDA = ', lda end if if ( m < 1 ) then ierror = ierror + 2 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGB_CHECK - Illegal M = ', m end if if ( ml < 0 .or. ml > min ( m, n ) - 1 ) then ierror = ierror + 4 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGB_CHECK - Illegal ML = ', ml end if if ( mu < 0 .or. mu > min ( m, n ) - 1 ) then ierror = ierror + 8 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGB_CHECK - Illegal MU = ', mu end if if ( n < 1 ) then ierror = ierror + 16 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGB_CHECK - Illegal N = ', n end if return end subroutine sgb_fa ( a, lda, n, ml, mu, ipivot, info ) ! !******************************************************************************* ! !! SGB_FA factors a matrix stored in LINPACK general band storage. ! ! ! Discussion: ! ! The matrix is stored in the array using LINPACK general band storage. ! The following program segment will set up the input. ! ! m = ml + mu + 1 ! do j = 1, n ! i1 = max ( 1, j-mu ) ! i2 = min ( n, j+ml ) ! do i = i1, i2 ! k = i - j + m ! a(k,j) = afull(i,j) ! end do ! end do ! ! This uses rows ML+1 through 2*ML+MU+1 of the array A. ! In addition, the first ML rows in the array are used for ! elements generated during the triangularization. ! The total number of rows needed in A is 2*ML+MU+1. ! The ML+MU by ML+MU upper left triangle and the ! ML by ML lower right triangle are not referenced. ! ! Modified: ! ! 04 March 1999 ! ! Parameters: ! ! Input/output, real A(LDA,N), the matrix in band storage. The ! columns of the matrix are stored in the columns of the array, ! and the diagonals of the matrix are stored in rows ML+1 through ! 2*ML+MU+1. On return, A has been overwritten by the LU factors. ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least 2*ML+MU+1. ! ! Input, integer N, the order of the matrix. ! N must be positive. ! ! Input, integer ML, MU, the lower and upper bandwidths. ! ML and MU must be nonnegative, and no greater than N-1. ! ! Output, integer IPIVOT(N), the pivot vector. ! ! Output, integer INFO, singularity flag. ! 0, no singularity detected. ! nonzero, the factorization failed on the INFO-th step. ! implicit none ! integer lda integer n ! real a(lda,n) integer i integer i0 integer ierror integer info integer ipivot(n) 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 t ! ! Check the dimensions. ! call sgb_check ( lda, n, n, ml, mu, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGB_FA - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions.' return end if m = ml + mu + 1 info = 0 ! ! Zero out the initial fill-in columns. ! j0 = mu + 2 j1 = min ( n, m ) - 1 do jz = j0, j1 i0 = m + 1 - jz do i = i0, ml a(i,jz) = 0.0E+00 end do end do jz = j1 ju = 0 do k = 1, n-1 ! ! Zero out the next fill-in column. ! jz = jz + 1 if ( jz <= n ) then do i = 1, ml a(i,jz) = 0.0E+00 end do end if ! ! Find L = pivot index. ! lm = min ( ml, n-k ) l = m do j = m+1, m+lm if ( abs ( a(j,k) ) > abs ( a(l,k) ) ) then l = j end if end do ipivot(k) = l + k - m ! ! Zero pivot implies this column already triangularized. ! if ( a(l,k) == 0.0E+00 ) then info = k write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGB_FA - Fatal error!' write ( *, '(a,i6)' ) ' Zero pivot on step ', info return end if ! ! Interchange if necessary. ! if ( l /= m ) then t = a(l,k) a(l,k) = a(m,k) a(m,k) = t end if ! ! Compute multipliers. ! do j = m+1, m+lm a(j,k) = - a(j,k) / a(m,k) end do ! ! Row elimination with column indexing. ! ju = max ( ju, mu+ipivot(k) ) ju = min ( ju, n ) mm = m do j = k+1, ju l = l - 1 mm = mm - 1 t = a(l,j) if ( l /= mm ) then a(l,j) = a(mm,j) a(mm,j) = t end if do i = 1, lm a(mm+i,j) = a(mm+i,j) + t * a(m+i,k) end do end do end do ipivot(n) = n if ( a(m,n) == 0.0E+00 ) then info = n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGB_FA - Fatal error!' write ( *, '(a,i6)' ) ' Zero pivot on step ', info end if return end subroutine sgb_sl ( a, lda, n, ml, mu, ipivot, b, job ) ! !******************************************************************************* ! !! SGB_SL solves a system factored by SGB_FA. ! ! ! Modified: ! ! 04 March 1999 ! ! Parameters: ! ! Input, real A(LDA,N), the LU factors from SGB_FA. ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least 2*ML+MU+1. ! ! Input, integer N, the order of the matrix. ! N must be positive. ! ! Input, integer ML, MU, the lower and upper bandwidths. ! ML and MU must be nonnegative, and no greater than N-1. ! ! Input, integer IPIVOT(N), the pivot vector from SGB_FA. ! ! Input/output, real B(N). ! On input, the right hand side vector. ! On output, the solution. ! ! Input, integer JOB. ! 0, solve A * x = b. ! nonzero, solve transpose ( A ) * x = b. ! implicit none ! integer lda integer n ! real a(lda,n) real b(n) integer ierror integer ipivot(n) integer j integer job integer k integer l integer la integer lb integer lm integer m integer ml integer mu real t ! ! Check the dimensions. ! call sgb_check ( lda, n, n, ml, mu, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGB_SL - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions!' return end if m = mu + ml + 1 ! ! Solve A * x = b. ! if ( job == 0 ) then ! ! Solve L * Y = B. ! if ( ml >= 1 ) then do k = 1, n-1 lm = min ( ml, n-k ) l = ipivot(k) t = b(l) if ( l /= k ) then b(l) = b(k) b(k) = t end if do j = 1, lm b(k+j) = b(k+j) + t * a(m+j,k) end do end do end if ! ! Solve U * X = Y. ! do k = n, 1, -1 b(k) = b(k) / a(m,k) lm = min ( k, m ) - 1 la = m - lm lb = k - lm t = -b(k) do j = 0, lm-1 b(lb+j) = b(lb+j) + t * a(la+j,k) end do end do ! ! Solve transpose ( A ) * X = B. ! else ! ! Solve transpose(U) * Y = B. ! do k = 1, n lm = min ( k, m ) - 1 la = m - lm lb = k - lm t = 0.0E+00 do j = 0, lm-1 t = t + a(la+j,k) * b(lb+j) end do b(k) = ( b(k) - t ) / a(m,k) end do ! ! Solve transpose(L) * X = Y. ! if ( ml >= 1 ) then do k = n-1, 1, -1 lm = min ( ml, n-k ) t = 0.0E+00 do j = 1, lm t = t + a(m+j,k) * b(k+j) end do b(k) = b(k) + t l = ipivot(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 shape ( dpsi, norder, psi, xi ) ! !******************************************************************************* ! !! SHAPE evaluates the shape functions at the local coordinate XI. ! ! ! Discussion: ! ! The shape functions are defined for a reference element which ! extends from -1 to 1. NORDER equally spaced nodes lie within ! the element, and one shape function is associated with each node. ! Shape function I must be 1 at node I, and zero at the other nodes. ! ! Modified: ! ! 18 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Output, real DPSI(NORDER), the value of the derivative of each ! shape function at the point XI. ! ! Input, integer NORDER, the order of the element. ! ! Output, real PSI(NORDER), the value of each of the shape functions ! at the point XI. ! ! Input, real XI, the point at which the shape functions are to be ! evaluated. XI is assumed to be a coordinate in the "reference ! system, in which the reference element comprises the interval [-1, 1]. ! implicit none ! integer norder ! real a real b real c real d real dpsi(norder) real e real f real psi(norder) real xi ! ! Constant shape function. ! if ( norder == 1 ) then psi(1) = 1.0E+00 dpsi(1) = 0.0E+00 ! ! Linear shape functions. ! Nodes at -1 and +1. ! else if ( norder == 2 ) then psi(1) = 0.5E+00 * ( 1.0E+00 - xi ) dpsi(1) = - 0.5E+00 psi(2) = 0.5E+00 * ( 1.0E+00 + xi ) dpsi(2) = 0.5E+00 ! ! Quadratics. ! Nodes at -1, 0, +1. ! else if ( norder == 3 ) then psi(1) = xi * ( xi - 1.0E+00 ) * 0.5E+00 dpsi(1) = xi - 0.5E+00 psi(2) = ( 1.0E+00 + xi ) * ( 1.0E+00 - xi ) dpsi(2) = - 2.0E+00 * xi psi(3) = xi * ( xi + 1.0E+00 ) * 0.5E+00 dpsi(3) = xi + 0.5E+00 ! ! Cubics. ! Nodes at -1, -1/3, +1/3, +1. ! else if ( norder == 4 ) then psi(1) = - ( 3.0E+00 * xi + 1.0E+00 ) * ( 3.0E+00 * xi - 1.0E+00 ) & * ( xi - 1.0E+00 ) / 16.0E+00 dpsi(1) = - ( 27.0E+00 * xi**2 - 18.0E+00 * xi - 1.0E+00 ) / 16.0E+00 psi(2) = 9.0E+00 * ( xi + 1.0E+00 ) * ( 3.0E+00 * xi - 1.0E+00 ) & * ( xi - 1.0E+00 ) / 16.0E+00 dpsi(2) = ( 81.0E+00 * xi**2 - 18.0E+00 * xi - 27.0E+00 ) / 16.0E+00 psi(3) = - 9.0E+00 * ( xi + 1.0E+00 ) * ( 3.0E+00 * xi + 1.0E+00 ) & * ( xi - 1.0E+00 ) / 16.0E+00 dpsi(3) = ( -81.0E+00 * xi**2 - 18.0E+00 * xi + 27.0E+00 ) / 16.0E+00 psi(4) = ( xi + 1.0E+00 ) * ( 3.0E+00 * xi + 1.0E+00 ) & * ( 3.0E+00 * xi - 1.0E+00 ) / 16.0E+00 dpsi(4) = - ( -27.0E+00 * xi**2 - 18.0E+00 * xi + 1.0E+00 ) / 16.0E+00 ! ! Quartics ! Nodes at -1, -1/2, 0, +1/2, 1. ! else if ( norder == 5 ) then psi(1) = 2.0E+00 * ( xi + 0.5E+00 ) * xi * ( xi - 0.5E+00 ) & * ( xi - 1.0E+00 ) / 3.0E+00 dpsi(1) = ( 16.0E+00 * xi**3 - 12.0E+00 * xi**2 - 2.0E+00 * xi + 1.0E+00 ) & / 6.0E+00 psi(2) = -8.0E+00 * ( xi + 1.0E+00 ) * xi * ( xi - 0.5E+00 ) & * ( xi - 1.0E+00 ) / 3.0E+00 dpsi(2) = -4.0E+00 * ( 8.0E+00 * xi**3 - 3.0E+00 * xi**2 - 4.0E+00 * xi & + 1.0E+00 ) / 3.0E+00 psi(3) = 4.0E+00 * ( xi + 1.0E+00 ) * ( xi + 0.5E+00 ) * ( xi - 0.5E+00 ) & * ( xi - 1.0E+00 ) dpsi(3) = 16.0E+00 * xi**3 - 10.0E+00 * xi psi(4) = -8.0E+00 * ( xi + 1.0E+00 ) * ( xi + 0.5E+00 ) * xi & * ( xi - 1.0E+00 ) / 3.0E+00 dpsi(4) = -4.0E+00 * ( 8.0E+00 * xi**3 + 3.0E+00 * xi**2 - 4.0E+00 * xi & - 1.0E+00 ) / 3.0E+00 psi(5) = 2.0E+00 * ( xi + 1.0E+00 ) * ( xi + 0.5E+00 ) * xi & * ( xi - 0.5E+00 ) / 3.0E+00 dpsi(5) = ( 16.0E+00 * xi**3 + 12.0E+00 * xi**2 - 2.0E+00 * xi - 1.0E+00 ) & / 6.0E+00 ! ! Quintics. ! Nodes at -1, -3/5, -1/5, +1/5, 3/5, 1. ! (Gee, I'll bet there's nine better ways of doing this!) ! else if ( norder == 6 ) then a = - 1.0E+00 b = - 0.6E+00 c = - 0.2E+00 d = + 0.2E+00 e = + 0.6E+00 f = + 1.0E+00 psi(1) = ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & / ( ( a - b ) * ( a - c ) * ( a - d ) * ( a - e ) * ( a - f ) ) psi(2) = ( xi - a ) * ( xi - c ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & / ( ( b - a ) * ( b - c ) * ( b - d ) * ( b - e ) * ( b - f ) ) psi(3) = ( xi - a ) * ( xi - b ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & / ( ( c - a ) * ( c - b ) * ( c - d ) * ( c - e ) * ( c - f ) ) psi(4) = ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - e ) * ( xi - f ) & / ( ( d - a ) * ( d - b ) * ( d - c ) * ( d - e ) * ( d - f ) ) psi(5) = ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - f ) & / ( ( e - a ) * ( e - b ) * ( e - c ) * ( e - d ) * ( e - f ) ) psi(6) = ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - e ) & / ( ( f - a ) * ( f - b ) * ( f - c ) * ( f - d ) * ( f - e ) ) dpsi(1) = ( & ( xi - c ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & + ( xi - b ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & + ( xi - b ) * ( xi - c ) * ( xi - e ) * ( xi - f ) & + ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - f ) & + ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - e ) ) / & ( ( a - b ) * ( a - c ) * ( a - d ) * ( a - e ) * ( a - f ) ) dpsi(2) = ( & ( xi - c ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - c ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - c ) * ( xi - d ) * ( xi - f ) & + ( xi - a ) * ( xi - c ) * ( xi - d ) * ( xi - e ) ) / & ( ( b - a ) * ( b - c ) * ( b - d ) * ( b - e ) * ( b - f ) ) dpsi(3) = ( & ( xi - b ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - d ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - d ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - d ) * ( xi - e ) ) / & ( ( c - a ) * ( c - b ) * ( c - d ) * ( c - e ) * ( c - f ) ) dpsi(4) = ( & ( xi - b ) * ( xi - c ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - c ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - e ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - e ) ) / & ( ( d - a ) * ( d - b ) * ( d - c ) * ( d - e ) * ( d - f ) ) dpsi(5) = ( & ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - f ) & + ( xi - a ) * ( xi - c ) * ( xi - d ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - d ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - f ) & + ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - d ) ) / & ( ( e - a ) * ( e - b ) * ( e - c ) * ( e - d ) * ( e - f ) ) dpsi(6) = ( & ( xi - b ) * ( xi - c ) * ( xi - d ) * ( xi - e ) & + ( xi - a ) * ( xi - c ) * ( xi - d ) * ( xi - e ) & + ( xi - a ) * ( xi - b ) * ( xi - d ) * ( xi - e ) & + ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - e ) & + ( xi - a ) * ( xi - b ) * ( xi - c ) * ( xi - d ) ) / & ( ( f - a ) * ( f - b ) * ( f - c ) * ( f - d ) * ( f - e ) ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SHAPE - Fatal error!' write ( *, '(a,i6)' ) ' Illegal order = ', norder stop end if return end subroutine solution_evaluate ( dpsi, dudx, gu, kind, maxelm, maxnod, maxnpe, & nel, nodes, psi, u, x, xi ) ! !******************************************************************************* ! !! SOLUTION_EVALUATE evaluates the solution U(X) and its derivative dUdX. ! ! ! Discussion: ! ! The solution is to be evaluated within element NEL, at the point ! with reference coordinate XI. ! ! Modified: ! ! 17 May 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Parameters: ! ! Output, real DPSI(MAXNPE), the value of the derivatives of the ! shape functions at the evaluation point. ! ! Output, real DUDX, the value of the derivative of the computed solution ! at the evaluation point. ! ! Input, real GU(MAXNOD), the coefficient vector defining the computed ! solution. ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer NEL, the index of the element within which the ! solution is to be evaluated. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Output, real PSI(MAXNPE), the value of the shape functions at the ! evaluation point. ! ! Output, real U, the value of the computed solution at the evaluation point. ! ! Input, real X(MAXNOD), the physical coordinates of the nodes. ! ! Input, real XI, the reference coordinate, within element NEL, of the ! evaluation point. ! implicit none ! integer maxelm integer maxnod integer maxnpe ! real dpsi(maxnpe) real dudx real dxidx real gu(maxnod) integer ie integer ig integer kind(maxelm) integer nel integer nodes(maxnpe,maxelm) integer norder real psi(maxnpe) real u real x(maxnod) real xi ! ! At the point XI, evaluate the shape functions PSI(XI) and their ! derivatives DPSI(XI). ! norder = kind(nel) + 1 call shape ( dpsi, norder, psi, xi ) ! ! Compute the derivative correction factor, d XI / d X. ! The values of DPSI are originally calculated on an interval [-1,1], ! and need to be adjusted to the interval [ X(nodes(1,nel)), X(nodes(n,nel)) ]. ! dxidx = 2.0E+00 / ( x(nodes(norder,nel)) - x(nodes(1,nel)) ) ! ! Now compute: ! ! U = sum G(IG) * PSI(IE)(XI) ! dU/dX = sum G(IG) * dPSI(IE)(XI)/dXI * dXI/dX ! ! where "IE" is a local element index, and IG is a global index. ! u = 0.0E+00 dudx = 0.0E+00 do ie = 1, norder ig = nodes(ie,nel) u = u + gu(ig) * psi(ie) dudx = dudx + gu(ig) * dpsi(ie) * dxidx end do return end subroutine solution_write ( dpsi, gu, iexact, infix, irpn, kind, maxelm, & maxfix, maxmat, maxnod, maxnpe, maxrpn, nelem, nodes, psi, title, x ) ! !******************************************************************************* ! !! SOLUTION_WRITE prints out the computed solution. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Workspace, real DPSI(MAXNPE), derivatives of the shape functions. ! ! Input, real GU(MAXNOD), the coefficient vector defining the computed ! solution. ! ! Input, integer IEXACT, records what is known: ! 0, no formulas for the exact solution are known; ! 1, a formula for the exact value of U(X) is known; ! 2, formulas for the exact values of U(X) and dU/dX are known. ! ! Input, integer KIND(MAXELM), the degree of each element. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer MAXFIX, the maximum number of characters in INFIX, ! used to contain a raw formula. 80 should be enough. ! ! Input, integer MAXMAT, the maximum number of materials. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXNPE, the maximum number of nodes per element. ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODES(MAXNPE,MAXELM), the indices of the nodes ! belonging to each element. ! ! Workspace, real PSI(MAXNPE), shape function values. ! ! Input, character ( len = * ) TITLE, the title of the problem. ! ! Input, real X(MAXNOD), the physical coordinates of the nodes. ! implicit none ! integer maxelm integer maxfix integer maxnod integer maxnpe integer maxrpn ! character com real dpsi(maxnpe) real dudx real duhdx real gu(maxnod) integer i integer i1 integer i2 integer ierror integer iexact integer ifrm character ( len = * ) infix integer irpn(maxrpn) integer j integer jhi integer kind(maxelm) integer maxmat character ( len = 10 ) namvar integer nel integer nelem integer nodes(maxnpe,maxelm) real psi(maxnpe) character ( len = 80 ) title real uh real ux real value real x(maxnod) real x1 real x2 real xi real xx ! namvar = ' ' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'The solution has been computed.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ' if ( iexact == 0 ) then write ( *, '(a)' ) ' X U dUdx ' else if ( iexact==1 ) then write ( *, '(a)' ) ' X U computed U exact dUdx computed' else if ( iexact == 2 ) then write ( *, '(a)' ) ' X U computed U exact dUdx computed' // & ' dUdx exact' end if write ( *, '(a)' ) ' ' ifrm = 1 do i = 1, nelem nel = i jhi = kind(i) if ( i == nelem ) then jhi = kind(i) + 1 end if i1 = nodes(1,i) x1 = x(i1) i2 = nodes(kind(i)+1,i) x2 = x(i2) do j = 1, jhi xi = - 1.0E+00 + 2.0E+00 * real ( j - 1 ) / real ( kind(i) ) xx = x1 + 0.5E+00 * ( 1.0E+00 + xi ) * ( x2 - x1 ) ! ! Set the current X value. ! if ( iexact >= 1 ) then namvar = 'X' com = 'V' value = xx call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) end if ! ! Evaluate the exact solution UX. ! if ( iexact >= 1 ) then ifrm = maxmat*4+1 com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) ux = value end if ! ! Evaluate the exact derivative DUDX. ! if ( iexact >= 2 ) then ifrm = maxmat*4+2 com = 'E' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, value ) dudx = value end if ! ! Evaluate the computed solution, UH, DUHDX. ! call solution_evaluate ( dpsi, duhdx, gu, kind, maxelm, maxnod, maxnpe, & nel, nodes, psi, uh, x, xi ) if ( iexact == 0 ) then write ( *, '(5g14.6)' ) xx, uh, duhdx else if ( iexact == 1 ) then write ( *, '(5g14.6)' ) xx, uh, ux, duhdx else if ( iexact == 2 ) then write ( *, '(5g14.6)' ) xx, uh, ux, duhdx, dudx end if end do end do return end subroutine symadd ( ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvr, nsym, nsymp, symbol, valsym ) ! !*********************************************************************** ! !! SYMADD adds a symbol name to the list of symbolic names. ! ! ! Modified: ! ! 01 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred, variable was added to list. ! 1, error occurred, variable was not added to list. ! ! Input/output, integer IFREE, the index of the next free ! memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed ! in VALSYM. ! ! Input, character ( len = MAXLEN ) NAMVR, the name of the variable. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUMDIM(2,MAXSYM). ! For each symbol I, NUMDIM(1,I) is the number of rows in its ! value, and NUMDIM(2,I) is the number of columns. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! implicit none ! integer, parameter :: maxlen = 20 ! integer maxsym integer maxval ! integer i integer ierror integer ifree integer indx integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer lennam character ( len = maxlen ) namvr integer nsym integer nsymp logical s_eqi character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! ! Get the length of the variable name. ! ierror = 0 lennam = len_trim ( namvr ) if ( lennam <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMADD - Error!' write ( *, '(a)' ) ' This variable name has zero length!' ierror = 1 return end if ! ! Check if the name is already in use. ! do i = 1, nsymp if ( s_eqi ( namvr, symbol(i) ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMADD - Error!' write ( *, '(a)' ) ' The name ', trim ( namvr ), ' is already in use.' ierror = 1 return end if end do ! ! Check for user overriding earlier use of same name. ! do i = nsymp+1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then symbol(i) = 'VOID' end if end do ! ! Is there enough space in SYMBOL for another symbol name? ! if ( nsym >= maxsym ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMADD - Error!' write ( *, '(a)' ) ' There is not enough memory for ' // trim ( namvr ) write ( *, '(a)' ) ' Perhaps the KILL or INIT commands would help.' ierror = 1 return end if ! ! Is there enough space in VALSYM for the symbol's value? ! if ( ifree > maxval ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMADD - Error!' write ( *, '(a)' ) ' There is not enough memory for ' // trim ( namvr ) write ( *, '(a)' ) ' Perhaps the KILL or INIT commands would help.' return end if ! ! Insert information about the symbol into various tables. ! nsym = nsym + 1 symbol(nsym) = namvr iopsym(nsym) = 0 iprsym(nsym) = 10 ipval(nsym) = ifree ifree = ifree + 1 indx = ipval(nsym) valsym(indx) = 0.0E+00 return end subroutine symlst ( ifree, intsym, iopsym, iprsym, ipval, irpn, maxfrm, & maxrpn, maxsym, maxval, namvr, nints, nsym, nsymp, nsyms, symbol, valsym ) ! !******************************************************************************* ! !! SYMLST prints out information about a variable. ! ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXRPN, the maximum number of entries in IRPN. ! ! Input, CHARACTER*10 NAMVR, the name of the variable ! for which information is desired. ! ! There are a few "special" values for NAMVR: ! ! blank, prints values of all user defined symbols. ! 'ALL', same as blank, plus values of all built in symbols. ! 'DEBUG', same as 'ALL', plus a table of built in symbols. ! implicit none ! integer, parameter :: maxlen = 10 ! integer maxrpn integer maxsym integer maxval ! integer i integer ierror integer ifree integer ifrm integer ihi integer ilo integer inc integer intsym(maxsym) integer iopsym(maxsym) integer ipoint integer iprsym(maxsym) integer ipval(maxsym) integer irpn(maxrpn) integer itemp integer j integer jhi integer jlo integer k integer khi integer klo integer lens logical s_eqi integer maxfrm character ( len = maxlen ) namtmp character ( len = maxlen ) namvr integer ncol integer ndim integer nints integer nrow integer nsym integer nsymp integer nsyms character ( len = 100 ) output character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! output = ' ' if ( s_eqi ( namvr,'debug') ) then nsyms = max(nsym,nsyms) write ( *, '(a)' ) ' Symbol Loc Prior Args Nrow Ncol' write ( *, '(a)' ) ' ' do i = 1, nsyms write ( *, '(a10,1x,3i6)' ) symbol(i),ipval(i),iprsym(i), & iopsym(i) end do write ( *, '(a,i6)' ) 'Next free memory location at ', ifree do i = 1, nsyms ndim = 1 jlo = ipval(i) jhi = jlo+ndim-1 do klo = jlo, jhi, 5 khi = min(klo+4,jhi) write ( *,'(1x,5g14.6)')(valsym(k),k = klo,khi) end do end do end if if ( s_eqi ( namvr,'debug').or. s_eqi ( namvr,'all') ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Predefined symbols:' write ( *, '(a)' ) ' ' inc = 6 ihi = ((nsymp-1)/inc)+1 do i = 1, ihi jhi = i+(inc-1)*ihi jhi = min(jhi,nsymp) write ( *,'(1x,7(a10,1x))') (symbol(j),j = i,jhi,ihi) end do end if if ( s_eqi ( namvr,'debug').or.s_eqi ( namvr,'all').or. namvr == ' ' ) then if ( nsym > nsymp ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'User-defined variables:' do i = nsymp+1, nsym nrow = 1 ncol = 1 ipoint = ipval(i) namtmp = symbol(i) call matwrt ( ierror, namtmp, ncol, nrow, valsym(ipoint) ) end do else write ( *, '(a)' ) 'There are NO user defined variables now.' end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Formulas:' write ( *, '(a)' ) ' ' do ifrm = 1, maxfrm+1 if ( ifrm == maxfrm+1 ) then if ( nints <= 0 ) go to 100 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Most recent infix formula.' write ( *, '(a)' ) ' ' ilo = 1 ihi = nints else ilo = (ifrm-1)*80+1 ihi = ilo+79 end if jlo = 0 write ( *, '(a)' ) ' ' do i = ilo, ihi if ( ifrm == maxfrm+1 ) then itemp = intsym(i) else itemp = irpn(i) end if if ( itemp<=0.or.itemp>maxsym)go to 80 sym = symbol(itemp) if ( (jlo+lens)>80 ) then write ( *, '(1x,a)' ) output output = ' ' jlo = 0 end if jhi = jlo+lens jlo = jlo+1 output(jlo:jhi) = sym(1:lens) jlo = jhi+1 if ( sym == '$')go to 90 80 continue end do 90 continue write ( *, '(1x,a)' ) output 100 continue end do write ( *, '(a)' ) ' ' return end if ! ! List a single variable. ! do i = 1, nsym if ( s_eqi ( namvr,symbol(i)) ) then nrow = 1 ncol = 1 ipoint = ipval(i) namtmp = namvr call matwrt ( ierror, namtmp, ncol, nrow, valsym(ipoint) ) return end if end do write ( *, '(a)' ) 'Warning! There is no variable named ' & // trim ( namvr ) return end subroutine symbol_value ( com, ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvar, ncol, nrow, nsym, nsymp, symbol, valsym, value ) ! !******************************************************************************* ! !! SYMBOL_VALUE evaluates, sets or deletes a variable. ! ! ! Discussion: ! ! SYMVAL accepts the name of a symbol in NAMVAR, and does one ! of three things, depending on the value of COM: ! ! 'R' - "reads" the value of the variable, returning it in VALUE. ! 'V' - "sets" the value of the variable from the array VALUE. ! 'D' - "deletes" the variable from memory. ! ! Modified: ! ! 15 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character COM, determines what task is to be done: ! 'R', "reads" the variable, returning its value and dimensions. ! 'V', "sets" the value and dimensions of the variable. ! 'D', "deletes" the variable from memory. ! ! Output, integer ierror, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input/output, integer IFREE, the index of the next free ! memory location in VALSYM. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(i) is 0. ! If a symbol represents a unary operator such as abs, ! IOPSYM(i) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the function priorities. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed. ! ! Input, character ( len = maxlen ) namvar, the name of the variable. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, character ( len = maxlen ) SYMBOL(MAXSYM), symbolic names. ! ! Input, real VALSYM(MAXVAL), the values of the symbolic variables. ! ! Input/output, real VALUE, holds the value of the variable. ! implicit none ! integer, parameter :: maxlen = 10 ! integer maxsym integer maxval ! character com integer i integer ierror integer ifree integer ihi integer ilo integer indx integer iopsym(maxsym) integer iprsym(maxsym) integer ipval(maxsym) integer j integer lennam logical s_eqi integer match character ( len = maxlen ) namvar integer ncol integer nrow integer nsym integer nsymp character ( len = maxlen ) symbol(maxsym) real valsym(maxval) real value ! ierror = 0 if ( len_trim ( namvar ) <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' Bad variable name ' // trim ( namvar ) return end if ! ! Search for a match between NAMVAR and any predefined symbol. ! match = 0 do i = 1, nsym if ( s_eqi ( namvar, symbol(i) ) ) then match = i exit end if end do ! ! If no such variable seems to exist, but we have been asked ! to DELETE or READ such a variable, then exit. ! if ( match == 0 ) then if ( s_eqi ( com, 'D' ) .or. s_eqi ( com, 'R' ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' No variable named ' // trim ( namvar ) if ( s_eqi ( com, 'R' ) ) then ierror = 1 return end if end if ! ! Add the name of the new variable to the internal list. ! call symadd ( ierror, ifree, iopsym, iprsym, ipval, & maxsym, maxval, namvar, nsym, nsymp, symbol, valsym ) if ( ierror /= 0 ) then return end if match = nsym ! ! The matched symbol has index MATCH. ! else ! ! If the symbol is not the name of a variable, but ! rather the name of an operator, then this is an error. ! if ( iopsym(match) /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' You are trying to assign a value to' write ( *, '(a)' ) ' the name of a function or operator:' write ( *, '(a)' ) symbol(match) return end if ! ! Some variables may not be revalued. ! if ( s_eqi ( com, 'V' ) ) then if ( s_eqi ( symbol(match), 'E' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'EPS' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'PI' ) ) ierror = 1 if ( symbol(match)(1:1) == '"' ) ierror = 1 end if ! ! Some variables may not be deleted. ! if ( s_eqi ( com, 'D' ) ) then if ( s_eqi ( symbol(match), 'E' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'EPS' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'PI' ) ) ierror = 1 if ( s_eqi ( symbol(match), 'SEED' ) ) ierror = 1 if ( symbol(match)(1:1) == '"' ) ierror = 1 end if if ( ierror /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE - Error!' write ( *, '(a)' ) ' Try to change special variable:' write ( *, '(a)' ) symbol(match) return end if nrow = 1 ncol = 1 end if indx = ipval(match) if ( s_eqi ( com, 'V' ) ) then valsym(indx) = value else if ( s_eqi ( com, 'R' ) ) then value = valsym(indx) else if ( s_eqi ( com, 'D' ) ) then ilo = indx ihi = ifree-1 do i = 1, ihi+1-ilo valsym(ilo-nrow*ncol+i-1) = valsym(ilo+i-1) end do ifree = ifree-nrow*ncol do i = match+1, nsym iopsym(i-1) = iopsym(i) iprsym(i-1) = iprsym(i) ipval(i-1) = ipval(i)-nrow*ncol symbol(i-1) = symbol(i) end do nsym = nsym-1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SYMBOL_VALUE:' write ( *, '(a)' ) ' Removed the variable ' // trim ( namvar ) write ( *, '(a,i6)' ) ' Free memory now: ', maxval+1-ifree end if return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date 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 character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, 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 ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine tokens ( ierror, ifinis, ifree, implic, indx1, indx2, ineg, infix, & intsym, iopsym, iprsym, ipval, maxsym, maxval, nints, nsym, nsymp, & nuvar, symbol, valsym ) ! !******************************************************************************* ! !! TOKENS parses a character string into recognized symbols and constants. ! ! ! Discussion: ! ! The routine produces the output array of integers INTSYM, containing ! whose entries stand for recognized symbols, variable names or ! constants. ! ! Modified: ! ! 02 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IERROR, error flag. ! 0, no error occurred. ! nonzero, an error occurred. ! ! Input, integer IFINIS, the index in SYMBOL of the symbol ! '$', meaning the end of the formula. ! ! Input, integer IFREE, the address of the next free ! memory location in VALSYM. ! ! Output, integer IMPLIC, implicit variable flag. ! 0, the assignment variable was not implicit. ! nonzero, the assignment variable was implicit, and has been ! added to the symbol table, and assigned the symbol index of IMPLIC. ! ! Input, integer INDX1, the index of the symbol 'INDEX1'. ! ! Input, integer INDX2, the index of the symbol 'INDEX2'. ! ! Input, integer INEG, the index of the symbol 'NEG'. ! ! Input, character ( len = * ) INFIX, a sequence of characters representing ! an infix formula. ! ! Output, integer INTSYM(80), a sequence of integers which are ! the indices of the symbols used in INFIX. ! ! Input, integer IOPSYM(MAXSYM), contains, for each ! symbol, the number of operands. In particular, if ! a symbol represents a constant, IOPSYM(I) is 0. ! If a symbol represents a unary operator such as ABS, ! IOPSYM(I) is 1. IOPSYM may be as large as 3. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Input, integer IPVAL(MAXSYM), contains, for each symbol, ! the address in VALSYM where the value of the symbol ! is stored, if it is a scalar. If the symbol represents ! a vector or matrix, IPVAL(IARG) points to location of ! the first entry. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed in VALSYM. ! ! Output, integer NINTS, the number of useful entries in INTSYM. ! ! Input, integer NSYM, the total number of symbols, including temporaries. ! ! Input, integer NSYMP, the number of permanent symbols. ! ! Input, integer NUVAR. ! If NUVAR is 0, then the formula may refer to variables which ! have not yet been declared. ! Otherwise, all variables on the right hand side of the "=" ! sign must have been previously declared. ! ! Input, character ( len = MAXLEN ) SYMBOL(MAXSYM), the symbolic names ! of the internally defined constants and functions. ! ! Input, real VALSYM(MAXVAL), contains the values of all the ! symbolic variables. ! implicit none ! integer, parameter :: maxlen = 20 ! integer maxsym ! character ctemp integer i integer iback integer icom1 integer icom2 integer ierror integer ifinis integer ifree integer ilpr integer implic integer indx1 integer indx2 integer ineg character ( len = * ) infix integer intsym(80) integer iopsym(maxsym) integer ipoint integer ipos integer iprsym(maxsym) integer ipval(maxsym) integer irpr integer j integer jback integer lchar integer lenfix integer lenmat integer lens integer loc integer loc1 integer loc2 integer match integer matcho integer maxval character ( len = maxlen ) namvr integer nints logical notnum integer nsym integer nsymp integer nuvar real rval logical s_eqi logical s_is_alpha character sym character ( len = maxlen ) sym1 character ( len = maxlen ) sym2 character ( len = maxlen ) sym3 character ( len = maxlen ) symbol(maxsym) real valsym(maxval) ! implic = 0 ierror = 0 nints = 0 match = 0 ipos = 1 lenfix = len_trim ( infix ) ! ! 1: Does this formula begin with a variable name followed ! by an equals sign? ! sym = '=' loc1 = index ( infix, sym ) if ( loc1 <= 1 ) then go to 20 end if ! ! The formula begins with a variable name followed by an equals sign. ! ! Pick off the name of the variable, but first make sure ! you stop before any left hand parenthesis, in case the ! formula is something like "alfred(7)=5". ! sym = '(' loc2 = index ( infix, sym ) if ( loc2 <= 0 ) then loc = loc1 else loc = min ( loc1, loc2 ) end if if ( loc <=1 .or. loc > maxlen+1 ) then go to 20 end if loc = loc - 1 namvr = infix(1:loc) ! ! Check to see whether the name of the assignment variable has ! already been defined. ! do i = 1, nsym if ( s_eqi ( namvr, symbol(i) ) ) then go to 20 end if end do ! ! The assignment variable has not previously been defined. ! Prepare to add it to the table of symbols. ! ipos = loc + 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS:' write ( *, '(a)' ) ' Implicit definition of ' // trim ( namvr ) call symadd ( ierror, ifree, iopsym, iprsym, ipval, maxsym, maxval, namvr, & nsym, nsymp, symbol, valsym ) if ( ierror /= 0 ) then return end if nints = nints + 1 intsym(nints) = nsym match = nsym implic = nsym ! ! Process the next chunk of the formula. ! 20 continue ! ! If we have reached the end of the input string, we're done. ! Append an "End-of-formula" marker "$" and return. ! if ( ipos > lenfix ) then nints = nints + 1 intsym(nints) = ifinis return end if ! ! We're not done yet. Remember the last match in MATCHO, and ! look for the next match. Go through all the symbols, and ! try to find the longest match possible. ! matcho = match match = 0 lenmat = 0 do i = 1, nsym lens = len_trim ( symbol(i) ) if ( lens > lenmat .and. ipos + lens - 1 <= lenfix ) then if ( s_eqi ( infix(ipos:ipos+lens-1), symbol(i) ) ) then lenmat = lens match = i end if end if end do if ( match == 0 ) then go to 110 end if ! ! We found a match. ! ! But watch out! The user might be implicitly defining a name, ! but because the first part matches an earlier symbol, you don't ! realize it! For example, in a formula involving SINGLE, you ! might think you see "SIN". So make sure that the match you've ! got is not followed by an alphabetic character. ! if ( iopsym(match) /= 0 ) then go to 40 end if if ( ipos+lenmat-1 >= lenfix ) then go to 40 end if ctemp = infix(ipos+lenmat:ipos+lenmat) if ( s_is_alpha ( ctemp ) ) then go to 140 end if ! ! We have matched an old symbol. ! 40 continue sym1 = symbol(match) sym2 = ' ' if ( matcho /= 0 ) then sym2 = symbol(matcho) end if ! ! Check for unary minus or plus. ! if ( sym1 == '-' ) then if ( sym2 == '**' .or. sym2 == ',' .or. & sym2 == '(' .or. sym2 == '=' .or. & ipos == 1 ) then sym = infix(ipos+1:ipos+1) if ( lge ( sym, '0' ) .and. lle ( sym, '9' ) ) then go to 110 end if if ( sym == '.' ) then go to 110 end if end if end if if ( sym1 == '-' ) then if ( ipos == 1 .or. sym2 == '(' .or. sym2 == ',' .or. sym2 == '=' ) then nints = nints + 1 intsym(nints) = ineg ipos = ipos + 1 go to 20 end if end if if ( sym1 == '+' ) then if ( ipos == 1 .or. sym2 == '(' .or. sym2 == '=' ) then ipos = ipos + 1 go to 20 end if end if ! ! If we've hit a right parenthesis, then we assume we're dealing with: ! ! a vector index: A(3) ! a matrix index: A(1,2) ! a two place operator: MAX(A,B), or MIN or POLY ! a three place operator: KRYLOV(A,X,N) ! ! We need to find the matching left parenthesis, count the ! arguments by assuming that commas separate them, and read ! off the name preceding the left parenthesis. ! if ( sym1 == ')' ) then icom1 = 0 icom2 = 0 ilpr = 0 irpr = 1 do iback = 1, nints i = nints + 1 - iback if ( intsym(i) <= 0 ) then go to 90 end if sym3 = symbol ( intsym(i) ) if ( sym3 == ',' .and. irpr-ilpr == 1 ) then if ( icom1 == 0 ) then icom1 = i else if ( icom2 == 0 ) then icom2 = i else ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' Too many commas in formula.' return end if else if ( sym3 == '(' ) then ilpr = ilpr + 1 else if ( sym3 == ')' ) then irpr = irpr + 1 end if ! ! If ILPR=IRPR, we've reached the matching left parenthesis. ! if ( ilpr == irpr ) then if ( i <= 1 ) then go to 100 end if if ( intsym(i-1) <= 0 ) then go to 100 end if sym3 = symbol ( intsym(i-1) ) if ( iopsym(intsym(i-1)) /= 0 .and. & .not. s_eqi ( sym3, 'MAX' ) .and. & .not. s_eqi ( sym3, 'MIN' ) .and. & .not. s_eqi ( sym3, 'ATAN2' ) .and. & .not. s_eqi ( sym3, 'POLY' ) .and. & .not. s_eqi ( sym3, 'KRYLOV' ) ) then go to 100 end if ! ! Copy ) into INTSYM. ! nints = nints + 1 intsym(nints) = match ! ! INDX1 is our name for a vector reference. ! if ( icom1 == 0 ) then match = indx1 ! ! INDX2 is our name for a matrix reference. ! else if ( icom2 == 0 ) then intsym(icom1) = match do jback = 1, nints-icom1 j = nints + 1 - jback intsym(j+1) = intsym(j) end do nints = nints + 1 intsym(icom1+1) = match - 1 if ( iopsym(intsym(i-1)) /= 0 ) then ipos = ipos + 1 go to 20 end if match = indx2 ! ! Insert )( to divide a trio of arguments. ! else intsym(icom1) = match do jback = 1, nints-icom1 j = nints + 1 - jback intsym(j+1) = intsym(j) end do intsym(icom1+1) = match - 1 nints = nints + 1 intsym(icom2) = match do jback = 1, nints-icom2 j = nints + 1 - jback intsym(j+1) = intsym(j) end do intsym(icom2+1) = match - 1 nints = nints + 1 if ( iopsym(intsym(i-1)) /= 0 ) then ipos = ipos + 1 go to 20 end if match = indx2 end if ipos = ipos + 1 nints = nints + 1 intsym(nints) = match go to 20 end if 90 continue end do end if 100 continue lens = len_trim ( symbol(match) ) nints = nints + 1 intsym(nints) = match ipos = ipos + lens go to 20 ! ! Check for a constant. ! 110 continue call s_to_r ( infix(ipos:), rval, ierror, lchar ) if ( lchar <= 0 ) then go to 140 end if if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' Illegal real number.' return end if ! ! Look for a constant already defined with the same value. ! do i = nsymp+1, nsym if ( iprsym(i) == 10 ) then if ( symbol(i)(1:1) == '"' ) then ipoint = ipval(i) if ( valsym(ipoint) == rval ) then match = i go to 130 end if end if end if end do ! ! This is a new constant. Make up a name for it. ! if ( real ( int ( rval ) ) == rval ) then write ( namvr, '(''"'',i19)' ) int ( rval ) else write ( namvr, '(''"'', g19.10)' ) rval end if call s_blank_delete ( namvr ) ! ! Add the constant to the symbol list. ! call symadd ( ierror, ifree, iopsym, iprsym, ipval, maxsym, maxval, namvr, & nsym, nsymp, symbol, valsym ) if ( ierror /= 0 ) then return end if ipoint = ipval(nsym) valsym(ipoint) = rval match = nsym 130 continue ! ! Advance LCHAR characters, except that CHRCTR will read in a ! comma as the end of the number, and we have to back up in ! such a case. ! if ( infix(ipos+lchar-1:ipos+lchar-1) /= ',' ) then ipos = ipos + lchar else ipos = ipos + lchar - 1 end if nints = nints + 1 intsym(nints) = match go to 20 ! ! Consider implicit variable declaration on right hand side. ! 140 continue if ( nuvar == 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' Undeclared variable in formula.' return end if namvr = ' ' do i = 1, maxlen sym = infix(ipos-1+i:ipos-1+i) notnum = llt ( sym, '0' ) .or. lgt ( sym, '0' ) if ( ( .not. s_is_alpha(sym) ) .and. notnum ) then exit end if namvr(i:i) = sym end do lchar = len_trim ( namvr ) if ( lchar <= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS - Error!' write ( *, '(a)' ) ' The string starting at location ', ipos, & ' is unreadable.' return end if ipos = ipos + lchar call symadd ( ierror, ifree, iopsym, iprsym, ipval, maxsym, maxval, namvr, & nsym, nsymp, symbol, valsym ) nints = nints + 1 intsym(nints) = nsym match = nsym write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TOKENS:' write ( *, '(a)' ) ' The undeclared variable ' // trim ( namvr ) //& ' will be assumed to be a scalar.' go to 20 end