program femod2 ! !******************************************************************************* ! !! FEMOD2 is a two dimensional interactive or batch finite element code. ! ! ! Discussion: ! ! FEMOD2 can solve boundary value problems on a two dimensional region using ! the finite element method. The problems this code can handle are of the ! following form- ! ! 1. Differential equation inside a region ! ! -DEL DOT ( K(X,Y) * DEL(U)) + B(X,Y) * U = F(X,Y) + P(X,Y) ! ! DEL is the 2-dimensional vector operator defined by ! DEL ( U(X,Y) ) = (DU/DX, DU/DY). ! DEL DOT DEL ( U(X,Y)) = D2 U/DX2 + D2 U/DY2 ! ! P(X,Y) stands for a point load function. ! ! 2. Boundary conditions ! ! If N is a node on the boundary, an essential boundary condition ! of the form U(N) = VALUE(X,Y) may be applied. The program allows you ! to specify such conditions for groups of nodes. You may ! specify a list of nodes, a side of the region, or the entire ! boundary. Then you specify the value to be assigned to ! these nodes as a function of X and Y. ! ! If E is an element with side s on the boundary, then the natural ! boundary condition -K*DU/DN = VALUE1(X,Y)*U - VALUE2(X,Y) may be ! applied on side S of E by specifying the element number E, ! the side S (E, N, W or S), and formulas for the quantities ! VALUE1(X,Y), VALUE2(X,Y). DU/DN is the outward normal ! derivative. Note that if you apply natural boundary conditions ! on all sides of the region, the solution is only determined ! up to a constant. The program resolves this problem by forcing ! the solution at node 1 to be 0. ! ! The point load function P(X,Y) may be defined to consist of ! several separate point loads, each with a distinct location ! and value. ! ! The program as now written will set up a mesh of NX by NY points, ! smoothly spaced in a quadrilateral whose four corners are defined ! by the user. ! ! 3. Program limitations - ! ! FEMOD2 cannot handle time dependent problems. ! ! The element types available are ! 3 node linear triangles, ! 4 node linear quadrilaterals, ! 6 node quadratic triangles, ! 8 node serendipity quadrilaterals, ! 9 node quadratic quadrilaterals. ! ! Integration rules: ! Quadrilaterals: 4 and 9 point schemes ! Triangles: 3, 4 and 7 point schemes. ! ! Up to 5 material types. ! Up to 400 elements. ! Up to 800 nodes. ! Up to 20 point loads. ! Up to 200 essential boundary conditions. ! Up to 200 natural boundary conditions. ! ! The region must be a quadrilateral, with evenly spaced nodes in ! the X and Y directions, although there may be more in ! one direction than the other. ! ! Reference ! ! Eric Becker, Graham Carey, J Tinsley Oden, ! Finite Elements, An Introduction, Volume I. ! Prentice-Hall, Englewood Cliffs, New Jersey, 1981. ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer, parameter :: maxban = 10 integer, parameter :: max_ebc = 50 integer, parameter :: max_nbc = 50 integer, parameter :: maxelm = 100 integer, parameter :: maxfix = 80 integer, parameter :: maxgauss = 9 integer, parameter :: maxmat = 2 integer, parameter :: maxnod = 100 integer, parameter :: maxord = 9 integer, parameter :: maxpt = 20 integer, parameter :: maxrpn = 1600 ! character com character ( len = 10 ) command real dpsi(2,maxord) real ef(maxord) real ek(maxord,maxord) real gf(maxnod) real gk(maxban,maxnod) real gu(maxnod) integer ierror integer ifrm character ( len = 80 ) infix integer ipivot(maxnod) integer irhs integer irpn(maxrpn) character isay integer itrue integer mat(maxelm) integer ml integer mu character ( len = 10 ) namvar integer num_ebc integer num_nbc integer nel integer nel_nbc(max_nbc) integer nelem integer ngauss integer nmat integer nnode integer node_ebc(max_ebc) integer nodes(maxord,maxelm) integer norder integer npoint integer npt(maxpt) integer nrow integer nx integer ny real psi(maxord) logical s_eqi integer side_nbc(max_nbc) character ( len = 70 ) title real value real val_ebc(max_ebc) real val_nbc(2,max_nbc) real vpt(maxpt) real wgauss(maxgauss) character ( len = 20 ) what real x(2,maxnod) real xgauss(2,maxgauss) ! ierror = 0 ifrm = 0 nmat = 0 call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD2' write ( *, '(a)' ) ' Last modified on 19 June 2000' ! ! Initialize the compiler. ! com = 'i' namvar = ' ' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) ! ! Define the variable names X and Y. ! com = 'a' namvar = 'x' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) com = 'a' namvar = 'y' call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) ! ! Initialize data. ! itrue = 0 ml = 0 mu = 0 nrow = 0 do write ( *, '(a)' ) 'Enter command' read ( *, '(a)' ) command ! ! PREP ! Call the preprocessor. ! if ( s_eqi ( command, 'PREP' ) ) then call prep ( infix, irpn, mat, max_ebc, max_nbc, maxelm, maxfix, & maxgauss, maxmat, maxnod, maxord, maxpt, maxrpn, num_ebc, num_nbc, & nel_nbc, nelem, ngauss, nmat, nnode, node_ebc, nodes, norder, npoint, & npt, side_nbc, nx, ny, title, val_ebc, val_nbc, vpt, x, wgauss, xgauss ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD2 - Fatal error!' write ( *, '(a)' ) ' Error during preprocessor!' exit end if ! ! PRINT ! EBC ! EKF ! ELEM ! GKF ! NBC ! NODE ! POINT ! else if ( command == 'PRINT' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Print WHAT?' read ( *, '(a)' ) what if ( what == 'EBC' ) then call ebc_write ( max_ebc, num_ebc, node_ebc, val_ebc ) else if ( what == 'EKF' ) then write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Enter an element number between 1 and ', nelem read ( *, * ) nel ! ! NEED EKF_SET call here. ! call ekf_write ( ef, ek, maxelm, maxord, nel, nodes, norder ) else if ( what == 'ELEM' ) then call element_write ( mat, maxelm, maxord, nelem, nodes, norder ) else if ( what == 'GKF' ) then irhs = 1 call gkf_write ( gk, ml, mu, nnode, nrow, irhs, gf ) else if ( what == 'NBC' ) then call nbc_write ( max_nbc, num_nbc, nel_nbc, norder, side_nbc, val_nbc ) else if ( what == 'NODE' ) then call node_write ( maxnod, nnode, x ) else if ( what == 'POINT' ) then call point_write ( maxpt, npoint, npt, vpt ) end if ! ! PROC ! Call the processor. ! else if ( command == 'PROC' ) then call proc ( dpsi, ef, ek, gf, gk, gu, ierror, infix, ipivot, irpn, mat, & max_ebc, max_nbc, maxelm, maxfix, maxgauss, maxnod, maxord, maxpt, & maxban, maxrpn, ml, mu, num_ebc, num_nbc, nel_nbc, nelem, ngauss, & nnode, node_ebc, nodes, norder, npoint, npt, nrow, side_nbc, psi, & val_ebc, val_nbc, vpt, wgauss, x, xgauss ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD2 - Fatal error!' write ( *, '(a)' ) ' Error during processing!' exit end if ! ! POST ! Call the postprocessor. ! else if ( command == 'POST' ) then call post ( dpsi, gu, infix, irpn, itrue, mat, maxelm, maxfix, & maxgauss, maxmat, maxnod, maxord, maxrpn, nelem, ngauss, nodes, & norder, nx, ny, psi, title, wgauss, x, xgauss ) ! ! QUIT ! else if ( command(1:1) == 'Q' ) then exit ! ! SOLUTION ! Enter formula for the solution. ! else if ( command == 'SOLUTION' ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter Y if true solution is known.' read ( *, '(a)' ) isay if ( s_eqi ( isay, 'Y' ) ) then call exact_read ( infix, irpn, itrue, maxfix, maxmat, maxrpn ) end if end if end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEMOD2:' write ( *, '(a)' ) ' Normal end of execution.' stop end subroutine addmat ( a, ml, mu, neqn, nrow, entry, i, j ) ! !******************************************************************************* ! !! ADDMAT adds a given value to the band matrix A in position (I,J). ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(NROW,NEQN), the matrix to be modified. ! On output, A(I,J) has been increased by ENTRY. ! ! Input, integer ML, MU, the lower and upper bandwidths of the matrix. ! ! Input, integer NEQN, the number of equations, and rows of the matrix. ! ! Input, integer NROW, the leading dimension of the matrix, which must ! be at least 2*ML+MU+1. ! ! Input, real ENTRY, the value to be added to the matrix entry A(I,J). ! ! Input, integer I, J, the row and column of the entry to be incremented. ! implicit none ! integer neqn integer nrow ! real a(nrow,neqn) real entry integer i integer j integer ml integer mu ! if ( entry == 0.0E+00 ) then return end if ! ! Check range of I and J. ! if ( i <= 0 .or. i > neqn .or. j <= 0 .or. j > neqn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ADDMAT - Fatal error!' write ( *, '(a)' ) ' Illegal indices:' write ( *, '(a,i6)' ) ' I = ', i write ( *, '(a,i6)' ) ' J = ', j stop end if if ( (j-i) <= mu .and. (i-j) <= ml ) then a(i-j+ml+mu+1,j) = a(i-j+ml+mu+1,j) + entry else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ADDMAT - Fatal error!' write ( *, '(a)' ) ' Illegal indices:' write ( *, '(a,i6)' ) ' I = ', i write ( *, '(a,i6)' ) ' J = ', j write ( *, '(a,i6)' ) ' ML = ', ml write ( *, '(a,i6)' ) ' MU = ', mu write ( *, '(a,g14.6)' ) ' ENTRY = ', entry stop end if 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, ncol, nrow, 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. ! ! 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. ! ! GCD(I,J) Greatest common divisor 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. ! Here "small" means less than RTOL in absolute ! value. ! ! RTOL The tolerance for the ROUND command, initially set ! to EPS. To change the value of RTOL, simply ! reassign it with a statement like "RTOL = 0.1E-6". ! ! 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. ! ! ZEROS(I) Matrix zero of order I. ! ZEROS(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. ! !******************************************************************************* ! ! 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, variables already 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 COMRPN's internal ! 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. ! nonzero means some kind of error. These errors are usually ! signaled by a printed message. ! ! 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 user's formula to be compiled and evaluated. ! ! Workspace, 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 ( len = * ) 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 in the variable. ! ! Input/output, integer NROW, the number of rows in the variable. ! ! Workspace, character ( len = 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. ! 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 IPRSYM(MAXSYM), contains, for each ! symbol, the operator "priority" of that symbol. This ! is simply to help COMRPN deal with ambiguous formulas like A*B+C. ! ! 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 MAXIW, controls the size of the integer ! work vector IWORK, which is needed by EIGEN, SGECO, and SGEFA. ! ! 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 = 2500 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, save, dimension ( 80 ) :: intsym integer, save, dimension ( maxsym ) :: iopsym integer, save, dimension ( maxsym ) :: iprsym integer, save, dimension ( maxsym ) :: ipval integer, save :: irad integer irpn(maxrpn) integer, save :: irtol integer, save, dimension ( maxsym ) :: 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 ) :: symbol real value real, dimension ( maxval ), save :: valsym ! ierror = 0 implic = 0 namvr = namvar maxfix = len ( infix ) if ( maxfix > 0 ) then maxfrm = ( maxrpn / maxfix ) else maxfrm = 0 end if nrow = max ( nrow, 1 ) ncol = max ( ncol, 1 ) call s_blank_delete ( namvr ) 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, irtol, 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 ) if ( ierror == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Note:' write ( *, '(a)' ) ' Added the variable ' // trim ( namvr ) end if 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 ) ! ! Reject the formula if it has no information in it. ! 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 did not pass parentheses checks.' return end if ! ! Convert INFIX, the string of characters, into INTSYM, an 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' end if return end if ! ! Convert INTSYM, which is an infix formula, into IRPN, which ! is an RPN formula. ! imin = ( ifrm - 1 ) * 80 + 1 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 ! ! Check that operators and arguments are balanced. ! ihi = imin + nrpn - 1 call rpnchk ( ierror, ihi, ilo, iopsym, irpn, maxrpn, maxsym ) if ( ierror /= 0 .or. ilo /= imin ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' Formula did not pass RPN checks.' 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 go to 10 end if else if ( s_eqi ( com, 'E' ) ) then go to 10 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMRPN - Fatal error!' write ( *, '(a)' ) ' Unknown command = ' // trim ( com ) ierror = 1 end if return ! ! COM=E, Evaluate formula. ! 10 continue imin = ( ifrm - 1 ) * 80 + 1 nrpn = 80 ifreesv = ifree call rpnval ( ierror, ifree, imin, iopsym, iprsym, ipval, & irad, irpn, irtol, istack, maxrpn, maxsym, maxval, & nrpn, nsym, nsyms, symbol, valsym, value ) ifree = ifreesv if ( ierror /= 0 ) then nrow = 1 ncol = 1 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' 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 ebc_apply ( gf, gk, max_ebc, maxnod, ml, mu, num_ebc, nnode, & node_ebc, nrow, val_ebc ) ! !******************************************************************************* ! !! EBC_APPLY enforces essential boundary conditions by modifying the system. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real GF(MAXNOD), the global force vector. ! ! Input/output, real GK(NROW,NNODE), the global stiffness matrix. ! ! Input, integer MAX_EBC, the maximum number of essential boundary ! conditions. ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer ML, MU, the lower and upper bandwidths of the matrix. ! ! Input, integer NUM_EBC, the number of essential boundary conditions. ! ! Input, integer NNODE, the number of nodes. ! ! Input, integer NODE_EBC(MAX_EBC), the nodes associated with essential ! boundary conditions. ! ! Input, integer NROW, the leading dimension of GK. ! ! Input, real VAL_EBC(MAX_EBC), the values associated with essential ! boundary conditions. ! implicit none ! integer max_ebc integer maxnod integer nnode integer nrow ! real entry real gf(maxnod) real gk(nrow,nnode) integer i integer irow integer j integer jcol integer ml integer mu integer num_ebc integer node_ebc(max_ebc) real val_ebc(max_ebc) ! ! Apply the essential boundary conditions. ! ! Replace the original equation for the variable associated with ! node NODE_EBC(I) by the equation ! ! V(NODE) = value. ! ! To preserve the symmetry of the stiffness matrix, zero out columns ! corresponding to the essential boundary conditions. Note that we ! must do this first, in case compressed storage is used. ! if ( num_ebc > 0 ) then do i = 1, num_ebc jcol = node_ebc(i) do j = 1, nnode irow = j if ( irow /= jcol ) then call getmat ( gk, ml, mu, nnode, nrow, entry, irow, jcol ) gf(irow) = gf(irow) - entry * val_ebc(i) entry = 0.0E+00 call setmat(gk, ml,mu,nnode,nrow,entry,irow,jcol) end if end do end do do i = 1, num_ebc irow = node_ebc(i) gf(irow) = val_ebc(i) entry = 0.0E+00 do j = 1, nnode jcol = j call setmat ( gk, ml,mu,nnode,nrow,entry,irow,jcol) end do entry = 1.0E+00 jcol = irow call setmat ( gk, ml,mu,nnode,nrow,entry,irow,jcol) end do else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Artificial boundary condition applied.' write ( *, '(a)' ) 'U(Node 1) = 0' irow = 1 entry = 0.0E+00 do j = 1, nnode jcol = j call setmat ( gk, ml, mu, nnode, nrow, entry, irow, jcol ) call setmat ( gk, ml, mu, nnode, nrow, entry, jcol, irow ) end do entry = 1.0E+00 jcol = irow call setmat ( gk, ml,mu,nnode,nrow,entry,irow,jcol) gf(irow) = 0.0E+00 end if return end subroutine ebc_read ( infix, irpn, max_ebc, maxfix, maxmat, maxnod, maxrpn, & num_ebc, nnode, node_ebc, norder, nx, ny, val_ebc, x ) ! !******************************************************************************* ! !! EBC_READ reads the essential boundary conditions. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_EBC, the maximum number of essential boundary ! conditions. ! ! Output, integer NUM_EBC, the number of essential boundary conditions. ! implicit none ! integer max_ebc integer maxfix integer maxnod integer maxrpn ! real bcval character com integer i integer ierror integer ifirst integer ifrm integer ihi integer ilast integer ilo integer inc character ( len = * ) infix integer inode integer irpn(maxrpn) character isay integer maxmat character ( len = 10 ) namvar integer num_ebc integer nnode integer nod integer node_ebc(max_ebc) integer norder integer nx integer ny integer nyp1h logical s_eqi real value real val_ebc(max_ebc) real x(2,maxnod) real xval real yval ! namvar = ' ' ! ! Read essential boundary condition data ! num_ebc = 0 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Specify essential boundary conditions ' write ( *, '(a)' ) 'of the form u(node) = value1' 10 continue write ( *, '(a)' ) '0=done, L=you will list individual nodes' write ( *, '(a)' ) 'N, S, E, W = specify nodes on given side' write ( *, '(a)' ) 'A=all boundary nodes have same condition.' write ( *, '(a)' ) ' Enter 0, L, S, N, E, W or A' read ( *, '(a)' ) isay if ( isay == '0' ) then return end if if ( s_eqi ( isay, 'L' ) ) go to 20 if ( s_eqi ( isay, 'S' ) ) go to 40 if ( s_eqi ( isay, 'N' ) ) go to 40 if ( s_eqi ( isay, 'E' ) ) go to 40 if ( s_eqi ( isay, 'W' ) ) go to 40 if ( s_eqi ( isay, 'A' ) ) go to 60 go to 10 20 continue ifirst = num_ebc + 1 ilast = num_ebc do write ( *, '(a)' ) 'Enter node number, 0 when done.' read ( *, * ) nod if ( nod <= 0 .and. ilast < ifirst)go to 10 if ( nod <= 0 .and. ilast >= ifirst)go to 90 ilast = ilast + 1 num_ebc = num_ebc + 1 node_ebc(num_ebc) = nod end do 40 continue if ( s_eqi ( isay, 'W' ) ) then ilo = 1 ihi = ny inc = 1 else if ( s_eqi ( isay, 'E' ) ) then ilo = nnode-ny+1 ihi = nnode inc = 1 else if ( s_eqi ( isay, 'S' ) ) then ilo = 1 ihi = nnode-ny+1 inc = ny else if ( s_eqi ( isay, 'N' ) ) then ilo = ny ihi = nnode inc = ny end if if ( norder == 8 .and. s_eqi ( isay, 'N' ) )go to 44 if ( norder == 8 .and. s_eqi ( isay, 'S' ) )go to 44 ifirst = num_ebc + 1 do i = ilo, ihi, inc num_ebc = num_ebc + 1 node_ebc(num_ebc) = i end do ilast = num_ebc go to 90 ! ! Handle north and south walls of 8-node quadrilaterals. ! East and west can be handled as above without change. ! 44 continue ifirst = num_ebc + 1 nyp1h = ( ny + 1 ) / 2 if ( s_eqi ( isay, 'S' ) ) then inode = 1 - nyp1h else inode = 0 end if do i = 1, nx if ( ( mod(i,2) == 1 .and. s_eqi ( isay, 'N' ) ) .or. & ( mod(i,2) == 0 .and. s_eqi ( isay, 'S' ) ) ) then inode = inode+ny else inode = inode+nyp1h end if num_ebc = num_ebc+1 node_ebc(num_ebc) = inode end do ilast = num_ebc go to 90 60 continue num_ebc = 0 if ( norder /= 8 ) then ifirst = num_ebc+1 do inode = 1,nnode if ( inode <= ny .or. inode > nnode-ny .or. mod ( inode, ny ) == 0 .or. & mod(inode,ny) == 1 ) then num_ebc = num_ebc+1 node_ebc(num_ebc) = inode end if end do ilast = num_ebc ! ! Special treatment for 8-node quadrilaterals ! else ifirst = num_ebc+1 do inode = 1,nnode if ( inode <= ny .or. inode > nnode-ny ) then num_ebc = num_ebc+1 node_ebc(num_ebc) = inode end if end do nyp1h = (ny+1)/2 inode = ny ihi = nx-1 do i = 2,ihi inode = inode+1 num_ebc = num_ebc+1 node_ebc(num_ebc) = inode if ( mod(i,2) == 1 ) then inode = inode+ny-1 else inode = inode+nyp1h-1 end if num_ebc = num_ebc+1 node_ebc(num_ebc) = inode end do ilast = num_ebc end if ! ! We now have a list of all the nodes to be assigned ! so now let's get and evaluate the boundary condition ! 90 continue write ( *, '(a)' ) 'Enter formula for boundary value' read ( *, '(a)' ) infix com = 'f' ifrm = maxmat*3+4 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) if ( ierror /= 0 ) then ierror = 0 write ( *, '(a)' ) 'The formula did not compile.' go to 90 end if do i = ifirst, ilast nod = node_ebc(i) xval = x(1,nod) com = 'V' namvar = 'X' value = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) yval = x(2,nod) com = 'V' namvar = 'Y' value = yval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) com = 'E' ifrm = maxmat*3+4 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) bcval = value val_ebc(i) = bcval end do if ( s_eqi ( isay, 'A' ) ) then return end if go to 10 end subroutine ebc_write ( max_ebc, num_ebc, node_ebc, val_ebc ) ! !******************************************************************************* ! !! EBC_WRITE prints out the essential boundary conditions. ! ! ! Modified: ! ! 25 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_EBC, the maximum number of essential boundary ! conditions. ! ! Input, integer NUM_EBC, the number of essential boundary conditions. ! ! Input, integer NODE_EBC(MAX_EBC), a list of the nodes at which essential ! boundary conditions are applied. ! ! Input, real VAL_EBC(MAX_EBC), the value of the essential boundary ! conditions. ! implicit none ! integer max_ebc ! integer i integer num_ebc integer node_ebc(max_ebc) real val_ebc(max_ebc) ! if ( num_ebc <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'There are no essential boundary conditions.' write ( *, '(a)' ) ' ' return end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Essential boundary conditions:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Index Node Assigned value' do i = 1, num_ebc write ( *, '(i4,4x,i4,4x,g14.6)' ) i, node_ebc(i), val_ebc(i) end do return end subroutine ekf_set ( ef, ek, ierror, infix, irpn, mater, maxfix, maxgauss, & maxrpn, norder, ngauss, wgauss, x, xgauss ) ! !******************************************************************************* ! !! EKF_SET sets up an element stiffness matrix and right hand side. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real EF(NORDER), the elemental right hand side vector. ! ! Output, real EK(NORDER,NORDER), the elemental stiffness matrix. ! ! Output, integer IERROR, error flag. ! 0, no errors were detected. ! nonzero, an error occurred. ! ! Input, character ( len = * ) INFIX, a string used for formulas. ! ! Input, integer IRPN(MAXRPN), contains "compiled" formulas. ! ! Input, integer MATER, the index of the element material. ! ! Input, integer MAXFIX, the maximum length of INFIX. ! ! Input, integer MAXGAUSS, the maximum number of Gauss integration points. ! ! Input, integer MAXRPN, the dimension of IRPN. ! ! Input, integer NORDER, the order of the element. ! ! Input, integer NGAUSS, the number of Gauss integration points. ! ! Input, real WGAUSS(NGAUSS), the Gauss integration weights. ! ! Input, real X(2,NORDER), the coordinates of the element nodes. ! ! Input, real XGAUSS(NGAUSS), the Gauss integration points. ! implicit none ! integer maxfix integer maxrpn integer maxgauss integer norder ! character com real detj real dpsi(2,9) real dpsix(9) real dpsiy(9) real dsdx(2,2) real dxds(2,2) real ef(norder) real ek(norder,norder) real fac integer i integer ierror integer ifrm character ( len = * ) infix integer irpn(maxrpn) integer j integer k integer l integer mater character ( len = 10 ) namvar integer ngauss real psi(9) real tol real value real wgauss(maxgauss) real x(2,norder) real xb real xf real xgauss(2,maxgauss) real xk real xval real yval ! ierror = 0 ifrm = 0 ! ! Initialize element arrays ! ef(1:norder) = 0.0E+00 ek(1:norder,1:norder) = 0.0E+00 ! ! Begin integration point loop ! do l = 1, ngauss call shape ( dpsi, norder, psi, xgauss(1,l), xgauss(2,l) ) call xsitox ( norder, x, xgauss(1,l), xval, yval ) ! ! Store values of X and Y. ! com = 'V' namvar = 'X' value = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) com = 'V' namvar = 'Y' value = yval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) ! ! Compute K(X,Y). ! com = 'E' ifrm = (mater-1)*3+1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) xk = value ! ! Compute B(X,Y). ! com = 'e' ifrm = (mater-1)*3+2 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) xb = value ! ! Compute F(X,Y). ! com = 'e' ifrm = (mater-1)*3+3 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) xf = value ! ! Calculate DXDS...Equation (5.3.6) ! do i = 1, 2 do j = 1, 2 dxds(i,j) = 0.0E+00 do k = 1, norder dxds(i,j) = dxds(i,j) + dpsi(j,k) * x(i,k) end do end do end do ! ! Calculate DSDX...Equation (5.2.5) ! detj = dxds(1,1) * dxds(2,2) - dxds(1,2) * dxds(2,1) detj = abs ( detj ) tol = 0.0001E+00 if ( abs ( detj ) <= tol ) then write ( *, 1020 ) detj go to 100 end if dsdx(1,1) = dxds(2,2) / detj dsdx(2,2) = dxds(1,1) / detj dsdx(1,2) = - dxds(1,2) / detj dsdx(2,1) = - dxds(2,1) / detj ! ! Calculate D(PSI)/DX...Equation (5.3.5) ! do i = 1, norder dpsix(i) = dpsi(1,i) * dsdx(1,1) + dpsi(2,i) * dsdx(2,1) dpsiy(i) = dpsi(1,i) * dsdx(1,2) + dpsi(2,i) * dsdx(2,2) end do ! ! Accumulate integration point value of integrals ! fac = detj * wgauss(l) do i = 1, norder ef(i) = ef(i) + fac * xf * psi(i) do j = 1, norder ek(i,j) = ek(i,j) + fac * ( xk * ( dpsix(i) * dpsix(j) + dpsiy(i) & * dpsiy(j) ) + xb * psi(i) * psi(j) ) end do end do end do return 100 continue ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'EKF_SET - Fatal error!' write ( *, '(a)' ) ' Bad element matrix:' write ( *, 1050 ) x(1,1:norder) write ( *, 1060 ) x(2,1:norder) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'dXdS:' write ( *, '(a)' ) ' ' do i = 1, 2 write ( *, '(2g14.6)' ) dxds(i,1:2) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DPSI:' write ( *, '(a)' ) ' ' do i = 1, 2 write ( *, * ) dpsi(i,1:norder) end do return 1020 format('fatal error - small determinant = ',g14.6) 1050 format('x(1,i) = ',9g14.6) 1060 format('x(2,i) = ',9g14.6) end subroutine ekf_write ( ef, ek, maxelm, maxord, nel, nodes, norder ) ! !******************************************************************************* ! !! EKF_WRITE writes an element matrix and right hand side. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxord integer norder ! real ef(norder) real ek(norder,norder) integer i integer j integer nel integer nodes(maxord,maxelm) ! write ( *, '(a)' ) ' ' write ( *, 1000 ) nel write ( *, 1010 ) (nodes(j,nel),j = 1,norder) write ( *, '(a)' ) ' ' do i = 1, norder write ( *, 1020 ) nodes(i,nel), (ek(i,j),j = 1,norder), ef(i) end do return 1000 format('element matrix ek, right hand side ef, element ',i6) 1010 format(5x,9(4x,i6,4x)) 1020 format(i6,10g12.4) end subroutine element_read ( ngauss, norder ) ! !******************************************************************************* ! !! ELEMENT_READ reads the element order. ! ! ! Modified: ! ! 18 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer NGAUSS, the number of Gauss quadrature points. ! ! Output, integer NORDER, the element order. ! implicit none ! integer ngauss integer norder ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Define the element order' do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Choose an element type by number of nodes.' write ( *, '(a)' ) ' Legal choices are ( 3, 4, 6, 8 or 9 ):' read ( *, * ) norder if ( norder == 3 .or. norder == 6 ) then ngauss = 4 exit end if if ( norder == 4 .or. norder == 8 .or. norder == 9 ) then ngauss = 9 exit end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Illegal choice for element order.' end do return end subroutine element_write ( mat, maxelm, maxord, nelem, nodes, norder ) ! !******************************************************************************* ! !! ELEMENT_WRITE prints the element definitions. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxord ! integer i integer mat(maxelm) integer n integer nelem integer nodes(maxord,maxelm) integer norder ! write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Number of nodes used per element is ', norder write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Element Material Node numbers' write ( *, '(a)' ) ' ' do n = 1, nelem write ( *, 1020 ) n, mat(n), (nodes(i,n),i = 1,norder) end do return 1020 format(i7,2x,i8,4x,9i5) end subroutine enorms ( dpsi, gu, infix, irpn, mat, maxelm, maxfix, maxgauss, & maxmat, maxnod, maxord, maxrpn, nelem, ngauss, nodes, norder, psi, wgauss, & x, xgauss ) ! !******************************************************************************* ! !! ENORMS computes the RMS and energy norm error in the approximate solution. ! ! ! Discussion: ! ! The calculation requires that the user supply formulas for the ! exact solution U and its partial derivatives dUdX and dUdY. ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real GU(MAXNOD), the global solution vector. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxord integer maxrpn ! character com real delu real deluh real deluuh real detj real dpsi(2,maxord) real dsdx(2,2) real dudx real dudy real duhdx real duhdy real dxds(2,2) real eu real eue real euh real euhe real euuh real euuhe real gu(maxnod) integer i integer ierror integer ifrm integer ii character ( len = * ) infix integer irpn(maxrpn) integer j integer k integer kk integer mat(maxelm) integer mater integer maxmat character ( len = 10 ) namvar integer nelem integer ngauss integer nodes(maxord,maxelm) integer norder real psi(maxord) real squ real sque real squh real squhe real squuh real squuhe real uh real ux real value real wgauss(maxgauss) real x(2,maxnod) real xb real xgauss(2,maxgauss) real xk real xval real xyelem(2,9) real yval ! ifrm = 0 ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Comparison of computed solution UH to true solution U' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'element l2(u) l2(uh) l2(error) energy(u) ' // & 'energy(uh) energy(error)' write ( *, '(a)' ) ' ' squ = 0.0E+00 squh = 0.0E+00 squuh = 0.0E+00 eu = 0.0E+00 euh = 0.0E+00 euuh = 0.0E+00 do i = 1, nelem mater = mat(i) do j = 1,norder k = nodes(j,i) xyelem(1,j) = x(1,k) xyelem(2,j) = x(2,k) end do sque = 0.0E+00 squhe = 0.0E+00 squuhe = 0.0E+00 eue = 0.0E+00 euhe = 0.0E+00 euuhe = 0.0E+00 ! ! For the K-th Gauss point... ! do k = 1, ngauss ! ! Evaluate the shape functions. ! call shape ( dpsi, norder, psi, xgauss(1,k), xgauss(2,k) ) ! ! Map the Gauss point to an X, Y point. ! call xsitox ( norder, xyelem, xgauss(1,k), xval, yval ) com = 'v' namvar = 'x' value = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) namvar = 'y' value = yval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) ! ! Compute K(X,Y). ! com = 'e' ifrm = (mater-1)*3+1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) xk = value ! ! Compute B(X,Y). ! com = 'e' ifrm = (mater-1)*3+2 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) xb = value ! ! Compute Uexact(X,Y). ! com = 'e' ifrm = maxmat*3+1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) ux = value ! ! Compute d Uexact/dX ! com = 'e' ifrm = maxmat*3+2 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) dudx = value ! ! Compute d Uexact/dY ! com = 'e' ifrm = maxmat*3+3 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) dudy = value ! ! Compute PSI then UH, DUHDX, DUHDY ! do ii = 1,2 do j = 1,2 dxds(ii,j) = 0.0E+00 do kk = 1,norder dxds(ii,j) = dxds(ii,j)+dpsi(j,kk)*xyelem(ii,kk) end do end do end do detj = abs ( dxds(1,1) * dxds(2,2) - dxds(1,2) * dxds(2,1) ) dsdx(1,1) = dxds(2,2) / detj dsdx(2,2) = dxds(1,1) / detj dsdx(1,2) = - dxds(1,2) / detj dsdx(2,1) = - dxds(2,1) / detj uh = 0.0E+00 duhdx = 0.0E+00 duhdy = 0.0E+00 do ii = 1,norder kk = nodes(ii,i) uh = uh + gu(kk) * psi(ii) duhdx = duhdx + gu(kk) * ( dpsi(1,ii) * dsdx(1,1) + dpsi(2,ii) & * dsdx(2,1) ) duhdy = duhdy + gu(kk) * ( dpsi(1,ii) * dsdx(1,2) + dpsi(2,ii) & * dsdx(2,2) ) end do ! ! Compute error terms ! delu = dudx**2 + dudy**2 deluh = duhdx**2 + duhdy**2 deluuh = (dudx-duhdx)**2 +(dudy-duhdy)**2 eue = eue + detj*wgauss(k)*(-xk*delu+xb*ux**2) euhe = euhe + detj*wgauss(k)*(-xk*deluh+xb*uh**2) euuhe = euuhe + detj*wgauss(k)*(-xk*deluuh+xb*(ux-uh)**2) sque = sque + detj*wgauss(k)*ux**2 squhe = squhe + detj*wgauss(k)*uh**2 squuhe = squuhe + detj*wgauss(k)*(ux-uh)**2 end do eu = eu+eue euh = euh+euhe euuh = euuh+euuhe squ = squ+sque squh = squh+squhe squuh = squuh+squuhe eue = sqrt(eue) euhe = sqrt(euhe) euuhe = sqrt(euuhe) squuhe = sqrt(squuhe) squhe = sqrt(squhe) sque = sqrt(sque) write ( *, 1020 ) i, sque, squhe, squuhe, eue, euhe, euuhe end do squ = sqrt(squ) squh = sqrt(squh) squuh = sqrt(squuh) eu = sqrt(eu) euh = sqrt(euh) euuh = sqrt(euuh) write ( *, '(a)' ) ' ' write ( *, * ) 'Total', squ, squh, squuh, eu, euh, euuh return 1020 format(i5,6g11.3) end subroutine exact_read ( infix, irpn, itrue, maxfix, maxmat, maxrpn ) ! !******************************************************************************* ! !! EXACT_READ reads a formula for the exact solution. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! implicit none ! integer maxfix integer maxrpn ! character com integer ierror integer ifrm character ( len = * ) infix integer irpn(maxrpn) integer itrue integer maxmat character ( len = 10 ) namvar real value ! namvar = ' ' write ( *, '(a)' ) 'Enter formula for u(x,y)' read ( *, '(a)' ) infix com = 'f' ifrm = maxmat*3+1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) write ( *, '(a)' ) 'Enter formula for du/dx' read ( *, '(a)' ) infix com = 'f' ifrm = maxmat*3+2 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) write ( *, '(a)' ) 'Enter formula for du/dy' read ( *, '(a)' ) infix com = 'f' ifrm = maxmat*3+3 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) itrue = 1 return end subroutine flux ( dpsi, gu, ielem, infix, irpn, mater, maxelm, maxfix, & maxgauss, maxnod, maxord, maxrpn, ngauss, nodes, norder, psi, x, xgauss ) ! !******************************************************************************* ! !! FLUX computes flux terms. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxord integer maxrpn ! character com real detj real dpsi(2,maxord) real dsdx(2,2) real dudx real dudy real dxds(2,2) real gu(maxnod) integer i integer ielem integer ierror integer ifrm integer ik character ( len = * ) infix integer irpn(maxrpn) integer j integer j1 integer k integer kk integer l integer mater character ( len = 10 ) namvar integer ngauss integer nodes(maxord,maxelm) integer norder real psi(maxord) real sigxh real sigyh real tol real value real x(2,maxnod) real xgauss(2,maxgauss) real xk real xval real xx(2) real yval ! ifrm = 0 xx(1) = 0.0E+00 xx(2) = 0.0E+00 ! ! Store values of X and Y ! xval = xx(1) com = 'v' namvar = 'x' value = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) yval = xx(2) com = 'v' namvar = 'y' value = yval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) ! ! Evaluate material property K(X,Y). ! com = 'e' ifrm = (mater-1)*3+1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) xk = value do l = 1, ngauss call shape ( dpsi, norder, psi, xgauss(1,l), xgauss(2,l) ) do i = 1, 2 xx(i) = 0.0E+00 do ik = 1,norder j1 = nodes(ik,ielem) xx(i) = xx(i) + psi(ik) * x(i,j1) end do do j = 1, 2 dxds(i,j) = 0.0E+00 do k = 1,norder j1 = nodes(k,ielem) dxds(i,j) = dxds(i,j) + dpsi(j,k) * x(i,j1) end do end do end do detj = dxds(1,1) * dxds(2,2) - dxds(1,2) * dxds(2,1) tol = 0.0001E+00 if ( abs(detj) <= tol ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FLUX - Fatal error!' write ( *, * ) ' Determinant too small, = ', detj return end if dsdx(1,1) = dxds(2,2) / detj dsdx(2,2) = dxds(1,1) / detj dsdx(1,2) = - dxds(1,2) / detj dsdx(2,1) = - dxds(2,1) / detj ! ! Calculate DUH/DX and DUH/DY ! dudx = 0.0E+00 dudy = 0.0E+00 do i = 1, norder kk = nodes(i,ielem) dudx = dudx + ( dpsi(1,i) * dsdx(1,1) + dpsi(2,i) * dsdx(2,1) ) * gu(kk) dudy = dudy + ( dpsi(1,i) * dsdx(1,2) + dpsi(2,i) * dsdx(2,2) ) * gu(kk) end do ! ! Calculate the stress SIGMAH ! sigxh = - xk * dudx sigyh = - xk * dudy write ( *, '(4(g14.6,5x))' ) xx(1), xx(2), sigxh, sigyh end do return end subroutine funscl ( angle_to_radian, arg1, arg2, ierror, result, rtol, 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, real RTOL, the rounding tolerance. ! ! 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 real rtol 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 ( *, * ) ' ', 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 ( *, * ) ' 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 ( *, * ) ' 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 ( *, * ) ' 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 ) < rtol ) 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 ( *, * ) ' 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 ( *, * ) ' Unrecognized function name:', trim ( sym ) end if return end subroutine funval ( iarg1, iarg2, iarg3, ierror, ifree, iopsym, iprsym, & ipset, ipval, irad, irtol, itemp, maxsym, maxval, nsyms, & sym, symbol, valsym, value ) ! !******************************************************************************* ! !! FUNVAL evaluates a scalar, vector or matrix valued function. ! ! ! Modified: ! ! 17 September 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 IFREE, the address 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 IOUNIT(4). ! IOUNIT(1) is the FORTRAN input unit. ! IOUNIT(2) is the standard output unit, while IOUNIT(3) and ! IOUNIT(4), if nonzero, are auxilliary output units. ! ! Input, integer IPRSYM(MAXSYM), the relative priorities ! of the functions. ! ! Output, integer IPSET. If the function was an 'INDX1' or ! 'INDX2', then IPSET is the relative offset of the value. ! ! 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 ITEMP, the index of the constant. ! ! Workspace, integer IWORK(MAXIW). ! ! Input, integer MAXIW, the amount of space in IWORK, which ! should be at least MAXROW. ! ! Input, integer MAXROW, the maximum number of rows and columns ! for matrices. ! ! Input, integer MAXSYM, the maximum number of symbols allowed. ! ! Input, integer MAXVAL, the maximum number of values allowed ! in VALSYM. ! ! 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 = MAXLEN ) SYM, the symbolic name of the function. ! ! 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. ! ! Workspace, real VALUE(MAXROW,MAXROW), space used to hold ! the value of a variable. ! implicit none ! integer, parameter :: maxlen = 20 ! integer maxsym integer maxval ! real adet real angle_to_radian real arg1 real arg2 real coef character ( len = 3 ) ctemp integer i integer iarg1 integer iarg2 integer iarg3 integer ierror integer ifree integer index1 integer index2 integer index3 integer index4 integer index5 integer index6 integer index7 integer indx integer info integer iopsym(maxsym) integer iprsym(maxsym) integer ipset integer ipval(maxsym) integer irad integer irtol integer itemp integer j integer job integer jseed integer jtemp logical s_eqi integer matz integer na integer nb integer ncoef integer ncol integer ncol1 integer ncol2 integer ndim integer ndim1 integer ndim2 integer npow integer nrow integer nrow1 integer nrow2 integer nsyms integer ntemp character ( len = 100 ) output real rcond real result real rmat_trace real rmat_norme real rmat_normf real rtol real sdot character ( len = maxlen ) sym character ( len = maxlen ) symbol(maxsym) real temp real temp1 real tnorm real v1 real v2 real valsym(maxval) real value ! 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 ( *, * ) 'Error! You may not change the value of '// symbol(iarg1) return end if end if ! ! Set dimensions for HILBERT, HILBINV, ID and ZEROS 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, 'ZEROS' ) ) then nrow1 = int ( valsym(ipval(iarg1)) ) ncol1 = nrow1 else nrow1 = 1 ncol1 = 1 end if ndim1 = nrow1*ncol1 index1 = ipval(iarg1) nrow2 = 0 ncol2 = 0 ndim2 = 0 index2 = 0 if ( iarg2 /= 0 ) then nrow2 = 1 ncol2 = 1 ndim2 = nrow2*ncol2 index2 = ipval(iarg2) end if index3 = 0 if ( iarg3 /= 0 ) then index3 = ipval(iarg3) end if 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) if ( iarg2 == 0 ) then arg2 = 0.0E+00 else arg2 = valsym(ipval(iarg2)) end if angle_to_radian = valsym(irad) rtol = valsym(irtol) call funscl ( angle_to_radian, arg1, arg2, ierror, result, rtol, 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 end if ! ! ILLEGAL MATRIX OR VECTOR OPERATIONS. ! NOT SUPPORTED IN THIS VERSION. ! return end subroutine getmat ( a, ml, mu, neqn, nrow, entry, i, j ) ! !******************************************************************************* ! !! GETMAT returns the (I,J) entry of a matrix stored in LINPACK band storage. ! ! ! Modified: ! ! 21 January 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real A(NROW,NEQN), the matrix to be examined. ! ! Input, integer ML, MU, the lower and upper bandwidths of the matrix. ! ! Input, integer NEQN, the number of equations, and rows of the matrix. ! ! Input, integer NROW, the leading dimension of the matrix, which must ! be at least 2*ML+MU+1. ! ! Output, real ENTRY, the value of entry A(I,J). ! ! Input, integer I, J, the row and column of the matrix to be returned. ! implicit none ! integer neqn integer nrow ! real a(nrow,neqn) real entry integer i integer j integer ml integer mu ! entry = 0.0E+00 ! ! Check index values. ! if ( i <= 0 .or. i > neqn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GETMAT - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of index I = ', i stop end if if ( j <= 0 .or. j > neqn ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GETMAT - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of index J = ', j stop end if if ( ( j - i ) > mu ) then entry = 0.0E+00 return end if if ( ( i - j ) > ml ) then entry = 0.0E+00 return end if entry = a(i-j+ml+mu+1,j) return end subroutine gk_bandwidth ( maxelm, maxord, ml, mu, nelem, nodes, norder ) ! !******************************************************************************* ! !! GK_BANDWIDTH computes the bandwidth of the global stiffness matrix. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Output, integer ML, MU, the lower and upper bandwidths of the matrix. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxord ! integer i integer irow integer j integer jcol integer ml integer mu integer nel integer nelem integer nodes(maxord,maxelm) integer norder ! ml = 0 mu = 0 do nel = 1, nelem 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 ( a, ml, mu, neqn, nrow, irhs, rhs ) ! !******************************************************************************* ! !! GKF_WRITE prints out a matrix. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer neqn integer nrow ! real a(nrow,neqn) character ( len = 13 ) ctemp(10) character ( len = 13 ) ctemp2 real entry integer i integer icopy integer ihi integer ilo integer inc integer incp1 integer incx integer irhs integer j integer jcopy integer jhi integer jlo integer ml integer mu integer nonzer real rhs(neqn) ! if ( irhs == 0 ) then incx = 5 incp1 = incx else incx = 4 incp1 = incx+1 end if write ( *, '(a)' ) ' ' do jlo = 1, neqn, incx jhi = min(jlo+incx-1,neqn) inc = min(incx,neqn+1-jlo) write ( *, '(a)' ) ' ' write(*,1010)jlo,jhi write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ' ilo = max(jlo-mu,1) ihi = min(jhi+ml,neqn) do i = ilo,ihi icopy = i nonzer = 0 do j = 1,inc jcopy = jlo-1+j call getmat(a, ml, mu,neqn,nrow,entry,icopy,jcopy) if ( entry/=0.0)nonzer = 1 write(ctemp2,1000)entry ctemp(j) = ctemp2 if ( -mu<=i-jcopy.and.i-jcopy<=ml)go to 10 ctemp(j) = ' ' 10 continue end do if ( inc < incx ) then do j = inc+1,incp1 ctemp(j) = ' ' end do end if if ( irhs /= 0 ) then write(ctemp2,1000)rhs(i) ctemp(incp1) = ctemp2 end if if ( nonzer == 1 ) then write(*,1040)i,(ctemp(j),j = 1,incp1) end if end do end do write ( *, '(a)' ) ' ' return 1000 format(g13.5) 1010 format(' columns ',i5,' to ',i5) 1040 format(i4,1x,9a13) 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. ! 10 continue ir = mod ( ip, iq ) if ( ir /= 0 ) then ip = iq iq = ir go to 10 end if i_gcd = iq 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 ! ! The absolute value of the integer goes into S(ILO:IHI). ! ipos = ihi 10 continue ! ! Find the last digit of IVAL, strip it off, and stick it into ! the string. ! 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 if ( ival /= 0 ) then go to 10 end if ! ! Fill the empties with zeroes. ! do i = ilo, ipos s(i:i) = '0' end do return end subroutine indata ( op, var, value ) ! !******************************************************************************* ! !! INDATA stores and retrieves data. ! ! ! Discussion: ! ! The routine works like a sort of COMMON block. It stores or returns ! the values of certain variables. Thus, it allows two routines ! to "communicate" without having to have data passed up and ! down the calling tree in argument lists. ! ! Modified: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OP, the name of the operation to ! be carried out. ! ! 'GET' means set VALUE to the value of VAR. ! 'INC' means increment the value of VAR by VALUE. ! 'SET' means set the value of VAR to VALUE. ! ! Input, character ( len = * ) VAR, the name of the variable ! to be operated on. ! ! For this version of INDATA, VAR may be 'LPAGE' or 'NLINE'. ! ! Input/output, integer VALUE. ! ! If OP='GET', then VALUE is an output quantity, and ! is equal to the current value of the variable VAR. ! ! If OP='INC', then VALUE is an input quantity, and is ! the amount to be added to the value of VAR. ! ! If OP='SET', then VALUE is an input quantity, and is ! the value to be assigned to VAR. ! implicit none ! logical s_eqi integer, save :: lpage = 22 integer, save :: nline = 0 character ( len = * ) op integer value character ( len = * ) var ! if ( s_eqi ( var, 'NLINE' ) ) then if ( s_eqi ( op, 'SET' ) ) then nline = value else if ( s_eqi ( op, 'GET' ) ) then value = nline else if ( s_eqi ( op, 'INC' ) ) then nline = nline + value end if else if ( s_eqi ( var, 'LPAGE' ) ) then if ( s_eqi ( op, 'SET' ) ) then lpage = value else if ( s_eqi ( op, 'GET' ) ) then value = lpage else if ( s_eqi ( op, 'INC' ) ) then lpage = lpage + value end if end if return end subroutine inicom ( ifinis, ifree, indx1, indx2, ineg, infix, intsym, iopsym, & iprsym, ipval, irad, irpn, irtol, 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 irtol 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.0, 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) = 'RTOL' iopsym(nsym) = 0 iprsym(nsym) = 10 valsym(nsym) = epsilon ( 1.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 else if ( s_eqi ( symbol(i), 'RTOL' ) ) then irtol = i end if end do call random_seed value = 0.0E+00 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 factor 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 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 material_read ( infix, irpn, maxfix, maxmat, maxrpn, nmat ) ! !******************************************************************************* ! !! MATERIAL_READ reads coefficients K, B, and F for each material. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! implicit none ! integer maxfix integer maxrpn ! character com character ctemp integer i integer ierror integer ifrm character ( len = * ) infix integer irpn(maxrpn) integer j integer maxmat character ( len = 10 ) namvar integer nmat real value ! write ( *, '(a)' ) 'Define K, B and F in -del(k*del(u))+b*u=f+p' namvar = ' ' ! ! Get number of materials ! 10 continue write ( *, 1040 ) maxmat read ( *, * ) nmat if ( nmat < 1 .or. nmat > maxmat ) then write ( *, '(a)' ) 'Not a legal choice.' go to 10 end if do i = 1, nmat write ( *, 1060 ) i do j = 1, 3 if ( j == 1 ) then ctemp = 'K' else if ( j == 2) then ctemp = 'B' else if ( j == 3) then ctemp = 'F' end if 20 continue write ( *, 1000 ) ctemp read ( *, '(a)' ) infix com = 'F' ifrm = (i-1)*3+j call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) if ( ierror /= 0 ) then write ( *, '(a)' ) 'The formula does not compile.' ierror = 0 go to 20 end if end do end do return 1000 format('formula for ',a1) 1040 format('number of materials between 1 and ',i2) 1060 format('define material',i2) end subroutine nbc_apply ( dpsi, ef, ek, gf, gk, max_nbc, maxelm, maxnod, maxord, & ml, mu, num_nbc, nel_nbc, nnode, nodes, norder, nrow, side_nbc, psi, & val_nbc, x ) ! !******************************************************************************* ! !! NBC_APPLY enforces natural boundary conditions by modifying the system. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer max_nbc integer maxelm integer maxnod integer maxord integer nnode integer norder integer nrow ! real detj real detj1 real detj2 real dpsi(2,maxord) real ef(norder) real ek(norder,norder) real entry real fac real g1 real g2 real gf(maxnod) real gk(nrow,nnode) integer i integer i1 integer ieqn integer igauss integer inode integer j integer j1 integer jnode integer jvar integer ml integer mu integer num_nbc integer nc integer nel integer nel_nbc(max_nbc) integer nj integer njg integer nod(9) integer nodbc integer nodes(maxord,maxelm) integer npts integer ns integer side_nbc(max_nbc) real psi(maxord) real val_nbc(2,max_nbc) real wbcint(2) real x(2,maxnod) real xbc(2,3) real xbcint(2,2) ! ! Apply the natural boundary conditions. ! do i = 1, num_nbc ! ! Pick out the nodes that lie on the side of the element where ! the boundary condition applies. ! nel = nel_nbc(i) ns = side_nbc(i) if ( norder == 3 .or. norder == 6 ) then nc = 3 else if ( norder == 4 .or. norder == 8 .or. norder == 9 ) then nc = 4 end if if ( norder == 3 .or. norder == 4 ) then npts = 2 nod(1) = ns if ( ns == nc ) then nod(2) = 1 else nod(2) = ns + 1 end if else npts = 3 nod(1) = ns nod(2) = ns + nc if ( ns == nc ) then nod(3) = 1 else nod(3) = ns + 1 end if end if ! ! Pick out the integration points to use. ! if ( norder == 3 .or. norder == 6 ) then g1 = ( sqrt ( 3.0E+00 ) - 1.0E+00 ) / ( 2.0E+00 * sqrt ( 3.0E+00 ) ) g2 = ( sqrt ( 3.0E+00 ) + 1.0E+00 ) / ( 2.0E+00 * sqrt ( 3.0E+00 ) ) if ( nod(1) == 1 ) then wbcint(1) = 0.5 xbcint(1,1) = g1 xbcint(2,1) = 0.0E+00 wbcint(2) = 0.5 xbcint(1,2) = g2 xbcint(2,2) = 0.0E+00 else if ( nod(1) == 2 ) then wbcint(1) = sqrt ( 2.0E+00 ) / 2.0E+00 xbcint(1,1) = g1 xbcint(2,1) = 1.0E+00 - g1 wbcint(2) = sqrt ( 2.0E+00 ) / 2.0E+00 xbcint(1,2) = g2 xbcint(2,2) = 1.0E+00 - g2 else if ( nod(1) == 3 ) then wbcint(1) = 0.5 xbcint(1,1) = 0.0E+00 xbcint(2,1) = g2 wbcint(2) = 0.5 xbcint(1,2) = 0.0E+00 xbcint(2,2) = g1 end if else g1 = 1.0E+00 / sqrt ( 3.0E+00 ) if ( nod(1) == 1 ) then xbcint(1,1) = -g1 xbcint(2,1) = -1.0E+00 xbcint(1,2) = g1 xbcint(2,2) = -1.0E+00 else if ( nod(1) == 2 ) then xbcint(1,1) = 1.0E+00 xbcint(2,1) = -g1 xbcint(1,2) = 1.0E+00 xbcint(2,2) = g1 else if ( nod(1) == 3 ) then xbcint(1,1) = g1 xbcint(2,1) = 1.0E+00 xbcint(1,2) = -g1 xbcint(2,2) = 1.0E+00 else if ( nod(1) == 4 ) then xbcint(1,1) = -1.0E+00 xbcint(2,1) = g1 xbcint(1,2) = -1.0E+00 xbcint(2,2) = -g1 end if wbcint(1) = 1.0E+00 wbcint(2) = 1.0E+00 end if ! ! Pick out nodal coordinates. Compute global coordinates of the 2 or 3 ! nodes that lie along the boundary side. ! do j = 1,npts nj = nod(j) njg = nodes(nj,nel) xbc(1,j) = x(1,njg) xbc(2,j) = x(2,njg) end do ! ! Calculate boundary integrals EK and EF ! do i1 = 1, norder ef(i1) = 0.0E+00 do j1 = 1, norder ek(i1,j1) = 0.0E+00 end do end do ! ! Using Gauss quadrature order 2. ! do igauss = 1, 2 call shape ( dpsi, norder, psi, xbcint(1,igauss), xbcint(2,igauss) ) detj1 = 0.0E+00 detj2 = 0.0E+00 do i1 = 1,npts nodbc = nod(i1) detj1 = detj1 + xbc(1,i1) * dpsi(igauss,nodbc) detj2 = detj2 + xbc(2,i1) * dpsi(igauss,nodbc) end do detj = sqrt ( detj1**2 + detj2**2 ) fac = detj * wbcint(igauss) do i1 = 1, norder ef(i1) = ef(i1) + val_nbc(2,i) * psi(i1) * fac do j1 = 1, norder ek(i1,j1) = ek(i1,j1) + val_nbc(1,i) * psi(i1) * psi(j1) * fac end do end do end do ! ! Add EK to GK and EF to GF. ! do inode = 1, norder ieqn = nodes(inode,nel) gf(ieqn) = gf(ieqn) + ef(inode) do jnode = 1, norder jvar = nodes(jnode,nel) entry = ek(i,j) call addmat ( gk, ml, mu, nnode, nrow, entry, ieqn, jvar ) end do end do end do return end subroutine nbc_read ( infix, irpn, max_nbc, maxelm, maxfix, maxmat, maxnod, & maxord, maxrpn, num_nbc, nel_nbc, nelem, nodes, norder, side_nbc, nx, ny, & val_nbc, x ) ! !******************************************************************************* ! !! NBC_READ reads the natural boundary conditions. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer max_nbc integer maxelm integer maxfix integer maxnod integer maxord integer maxrpn ! character com integer ierror integer ifrm integer ihi integer ilo integer inc character ( len = * ) infix integer irpn(maxrpn) character isay integer maxmat character ( len = 10 ) namvar integer num_nbc integer nel_nbc(max_nbc) integer nelem integer nodes(maxord,maxelm) integer none integer norder integer ns integer side_nbc(max_nbc) integer ntwo integer nx integer ny logical s_eqi real value real val_nbc(2,max_nbc) real x(2,maxnod) ! num_nbc = 0 namvar = ' ' ! ! Read natural boundary condition data. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Specify the natural boundary conditions' write ( *, '(a)' ) 'of the form: -k*du/dn = v1*u - v2' 120 continue write ( *, '(a)' ) '0=no more, l=list singly, a=all have same' write ( *, '(a)' ) 'e, n, w, s = specify condition on given side.' read ( *, '(a)' ) isay if ( isay == '0' ) return 39 continue write ( *, '(a)' ) 'Enter formula for v1' read ( *, '(a)' ) infix com = 'f' ifrm = maxmat*3+4 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) if ( ierror/=0 ) then ierror = 0 write ( *, '(a)' ) ' The formula did not compile.' go to 39 end if 38 continue write ( *, '(a)' ) 'Enter formula for v2' read ( *, '(a)' ) infix com = 'f' ifrm = maxmat*3+5 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) if ( ierror/=0 ) then ierror = 0 write ( *, '(a)' ) 'formula did not compile.' go to 38 end if if ( s_eqi ( isay, 'l' ) )go to 130 if ( s_eqi ( isay, 'e' ) )go to 170 if ( s_eqi ( isay, 'n' ) )go to 170 if ( s_eqi ( isay, 's' ) )go to 170 if ( s_eqi ( isay, 'w' ) )go to 170 if ( s_eqi ( isay, 'a' ) )go to 170 go to 120 130 continue write ( *, '(a)' ) 'Enter element number' read ( *, * ) ilo ihi = ilo inc = 1 if ( ilo<=0 .or. ilo>nelem ) then write ( *, '(a)' ) 'illegal element' go to 120 end if ! ! Careful here. Let user specify N, S, E, or W ! and you should check for triangles ! write ( *, '(a)' ) 'Enter side n, w, e or s' read ( *, '(a)' ) isay if ( norder == 3 .or. norder == 6 ) then if ( s_eqi ( isay, 'n' ) )ns = 1 if ( s_eqi ( isay, 's' ) )ns = 1 if ( s_eqi ( isay, 'e' ) )ns = 3 if ( s_eqi ( isay, 'w' ) )ns = 3 none = 1 if ( ns == 1)ntwo = 2 if ( ns == 3)ntwo = 3 else if ( s_eqi ( isay, 's' ) ) ns = 1 if ( s_eqi ( isay, 'e' ) ) ns = 2 if ( s_eqi ( isay, 'n' ) )ns = 3 if ( s_eqi ( isay, 'w' ) )ns = 4 none = ns ntwo = none+1 if ( ntwo>4)ntwo = 1 end if call setnbc ( ihi, ilo, inc, infix, irpn, max_nbc, maxelm, maxfix, maxmat, & maxnod, maxord, maxrpn, num_nbc, nel_nbc, nodes, none, ns, side_nbc, & ntwo, val_nbc, x ) if ( num_nbc < max_nbc )go to 120 return 170 continue if ( ( .not. s_eqi ( isay, 'w' ) ) .and. ( .not. s_eqi ( isay, 'a' ) ) ) then go to 190 end if if ( norder == 3 ) then ilo = 2 ihi = 2*(ny-1) inc = 2 ns = 3 none = 1 ntwo = 3 else if ( norder == 4 ) then ilo = 1 ihi = ny-1 inc = 1 ns = 4 none = 4 ntwo = 1 else if ( norder == 6 ) then ilo = 2 ihi = ny-1 inc = 2 ns = 3 none = 1 ntwo = 3 else if ( norder == 8 .or. norder == 9 ) then ilo = 1 ihi = (ny-1)/2 inc = 1 ns = 4 none = 4 ntwo = 1 end if call setnbc ( ihi, ilo, inc, infix, irpn, max_nbc, & maxelm, maxfix, maxmat, maxnod, maxord, maxrpn, num_nbc, & nel_nbc, nodes, none, ns, side_nbc, ntwo, val_nbc, x ) 190 continue if ( ( .not. s_eqi ( isay, 'e' ) ) .and. ( .not. s_eqi ( isay, 'a' ) ) ) then go to 210 end if if ( norder == 3 ) then ilo = 2*(nx-2)*(ny-1)+1 ihi = 2*(nx-1)*(ny-1)-1 inc = 2 ns = 3 none = 1 ntwo = 3 else if ( norder == 4 ) then ilo = (nx-2)*(ny-1)+1 ihi = (nx-1)*(ny-1) inc = 1 ns = 2 none = 2 ntwo = 3 else if ( norder == 6 ) then ilo = 2*((nx-1)/2-1)*((ny-1)/2)+1 ihi = (nx-1)*((ny-1)/2)-1 inc = 2 ns = 3 none = 1 ntwo = 3 else if ( norder == 8 .or. norder==9 ) then ilo = ((nx-1)/2-1)*((ny-1)/2)+1 ihi = (nx-1)*(ny-1)/4 inc = 1 ns = 2 none = 2 ntwo = 3 end if call setnbc ( ihi, ilo, inc, infix, irpn, max_nbc, & maxelm, maxfix, maxmat, maxnod, maxord, maxrpn, num_nbc, & nel_nbc, nodes, none, ns, side_nbc, ntwo, val_nbc, x ) 210 continue if ( ( .not. s_eqi ( isay, 's' ) ) .and. ( .not. s_eqi ( isay, 'a' ) ) ) then go to 230 end if if ( norder == 3 ) then ilo = 1 ihi = 2*(nx-2)*(ny-1)+1 inc = 2*(ny-1) ns = 1 none = 1 ntwo = 2 else if ( norder == 4 ) then ilo = 1 ihi = (nx-2)*(ny-1)+1 inc = ny-1 ns = 1 none = 1 ntwo = 2 else if ( norder == 6 ) then ilo = 1 ihi = 2*((nx-1)/2-1)*((ny-1)/2)+1 inc = ny-1 ns = 1 none = 1 ntwo = 2 else if ( norder == 8 .or. norder==9 ) then ilo = 1 ihi = ((nx-1)/2-1)*((ny-1)/2)+1 inc = (ny-1)/2 ns = 1 none = 1 ntwo = 2 end if call setnbc ( ihi, ilo, inc, infix, irpn, max_nbc, & maxelm, maxfix, maxmat, maxnod, maxord, maxrpn, num_nbc, & nel_nbc, nodes, none, ns, side_nbc, ntwo, val_nbc, x ) 230 continue if ( ( .not. s_eqi ( isay, 'n' ) ) .and. ( .not. s_eqi ( isay, 'a' ) ) ) then go to 250 end if if ( norder == 3 ) then ilo = 2*(ny-1) ihi = 2*(nx-1)*(ny-1) inc = 2*(ny-1) ns = 1 none = 1 ntwo = 2 else if ( norder == 4 ) then ilo = ny-1 ihi = (nx-1)*(ny-1) inc = ny-1 ns = 3 none = 3 ntwo = 4 else if ( norder == 6 ) then ilo = ny-1 ihi = (nx-1)*(ny-1)/2 inc = ny-1 ns = 1 none = 1 ntwo = 2 else if ( norder == 8 .or. norder==9 ) then ilo = (ny-1)/2 ihi = (nx-1)*(ny-1)/4 inc = (ny-1)/2 ns = 3 none = 3 ntwo = 4 end if call setnbc ( ihi, ilo, inc, infix, irpn, max_nbc, maxelm, maxfix, maxmat, & maxnod, maxord, maxrpn, num_nbc, nel_nbc, nodes, none, ns, side_nbc, & ntwo, val_nbc, x ) 250 continue if ( s_eqi ( isay, 'a' ) ) then return end if go to 120 end subroutine nbc_write ( max_nbc, num_nbc, nel_nbc, norder, side_nbc, val_nbc ) ! !******************************************************************************* ! !! NBC_WRITE prints out the natural boundary conditions. ! ! ! Modified: ! ! 25 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NORDER, the element order. ! implicit none ! integer max_nbc ! character ( len = 5 ) dir(4) integer i integer ielem integer itemp integer jtemp integer num_nbc integer nel_nbc(max_nbc) integer norder integer side_nbc(max_nbc) real val_nbc(2,max_nbc) ! if ( num_nbc > 0 ) then dir(1) = 'south' dir(2) = 'east' dir(3) = 'north' dir(4) = 'west' write ( *, '(a)' )' ' write ( *, '(a)' ) 'Natural boundary conditions' write ( *, '(a)' ) ' ' write ( *, '(a)' ) '-k*(du/dn) = value1(x,y)*u - value2(x,y)' write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Element Side Value1 Value2' write ( *, '(a)' ) ' ' do i = 1, num_nbc itemp = side_nbc(i) if ( norder == 3 .or. norder == 6 ) then ielem = nel_nbc(i) if ( mod(ielem,2) == 0 ) then if ( itemp == 1)jtemp = 3 if ( itemp == 3)jtemp = 2 else if ( itemp == 1)jtemp = 1 if ( itemp == 3)jtemp = 4 end if itemp = jtemp end if write ( *, 1020 ) nel_nbc(i), dir(itemp), val_nbc(1,i), val_nbc(2,i) end do end if return 1020 format(i7,1x,a5,g14.6,g14.6) end subroutine node_read ( maxnod, nnode, norder, nx, ny, x ) ! !******************************************************************************* ! !! NODE_READ reads the node data. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Output, integer NNODE, the number of nodes. ! ! Input, integer NORDER, the element order. ! ! Output, integer NX, NY, the number of nodes in the X and Y directions. ! ! Output, real X(2,MAXNOD), the coordinates of the nodes. ! implicit none ! integer maxnod ! integer i integer inode integer j integer nnode integer norder integer nx integer nxm1h integer ny integer nym1h real wne real wnw real wse real wsw real x(2,maxnod) real xne real xnw real xse real xsw real yne real ynw real yse real ysw ! ! Get number of points in X, Y direction ! 10 continue write ( *, '(a)' ) 'Enter the number of nodes in the X and Y directions:' read ( *, * ) nx, ny if ( nx<2 .or. ny<2.or.(norder>4.and.(nx<3.or.ny<3))) then write ( *, '(a)' ) 'The number of nodes is too small.' go to 10 end if if ( norder>4.and.(mod(nx,2)/=1 .or. mod(ny,2)/=1) ) then write ( *, '(a)' ) 'This element requires odd number of nodes.' go to 10 end if if ( norder /= 8 ) then nnode = nx*ny else nxm1h = (nx-1)/2 nym1h = (ny-1)/2 nnode = nx*ny-nxm1h*nym1h end if if ( nnode > maxnod ) then write ( *, 1060 ) maxnod go to 10 end if write ( *, '(a)' ) 'Enter x, y for north-west corner' read ( *, * ) xnw, ynw write ( *, '(a)' ) 'Enter x, y for north-east corner' read ( *, * ) xne, yne write ( *, '(a)' ) 'Enter x, y for south-west corner' read ( *, * ) xsw, ysw write ( *, '(a)' ) 'Enter x, y for south-east corner' read ( *, * ) xse, yse inode = 0 do j = 1, nx do i = 1, ny if ( .not. ( mod(i,2) == 0 .and. mod(j,2) == 0 .and. norder == 8 ) ) then wsw = real((ny-i)*(nx-j)) / real((nx-1)*(ny-1)) wse = real((ny-i)*(j-1)) / real((nx-1)*(ny-1)) wnw = real((i-1)*(nx-j)) / real((nx-1)*(ny-1)) wne = real((i-1)*(j-1)) / real((nx-1)*(ny-1)) inode = inode+1 x(1,inode) = wnw * xnw + wne * xne + wsw * xsw + wse * xse x(2,inode) = wnw * ynw + wne * yne + wsw * ysw + wse * yse end if end do end do return 1060 format('sorry, i can only handle up to',i6,' nodes') end subroutine node_write ( maxnod, nnode, x ) ! !******************************************************************************* ! !! NODE_WRITE prints the coordinates of the nodes. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer NNODE, the number of nodes. ! ! Input, real X(2,MAXNOD), the coordinates of the nodes. ! implicit none ! integer maxnod ! integer n integer nnode real x(2,maxnod) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Node X-coordinate Y-coordinate' write ( *, '(a)' ) ' ' do n = 1, nnode write ( *, 1010 ) n, x(1,n), x(2,n) end do return 1010 format(i4,2x,2g16.4) end subroutine point_apply ( gf, maxnod, maxpt, npoint, npt, vpt ) ! !******************************************************************************* ! !! POINT_APPLY modifies the linear system by applying point loads. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXNOD, the maximum number of nodes. ! ! Input, integer MAXPT, the maximum number of point loads. ! implicit none ! integer maxnod integer maxpt ! real gf(maxnod) integer i integer npoint integer npt(maxpt) real vpt(maxpt) ! do i = 1, npoint gf(npt(i)) = gf(npt(i)) + vpt(i) end do return end subroutine point_read ( maxpt, nnode, npoint, npt, vpt ) ! !******************************************************************************* ! !! POINT_READ reads point loads. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXPT, the maximum number of point loads. ! implicit none ! integer maxpt ! integer nnode integer npoint integer npt(maxpt) integer nval real vpt(maxpt) real vval ! npoint = 0 10 continue write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Enter a node at which a point load is applied,' write ( *, '(a)' ) 'or "0" to quit.' read ( *, * ) nval if ( nval <= 0 ) then return end if if ( nval > nnode ) then write ( *, '(a)' ) ' ' write ( *, * ) 'Your node number was too large.' write ( *, * ) 'The maximum legal value is ', nnode go to 10 end if write ( *, * ) 'Enter the value of the load at this point:' read ( *, * ) vval npoint = npoint + 1 npt(npoint) = nval vpt(npoint) = vval if ( npoint < maxpt ) then go to 10 else write ( *, * ) ' ' write ( *, * ) 'The program can handle no more point loads.' end if return end subroutine point_write ( maxpt, npoint, npt, vpt ) ! !******************************************************************************* ! !! POINT_WRITE prints the point load data. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXPT, the maximum number of point loads. ! implicit none ! integer maxpt ! integer i integer npoint integer npt(maxpt) real vpt(maxpt) ! if ( npoint <= 0 ) then return end if write(*,*)' ' write(*,*)'Point loads' write(*,*)' ' write(*,*)'Node Load' do i = 1, npoint write ( *, '(i4,4x,g14.6)' ) npt(i), vpt(i) end do return end subroutine post ( dpsi, gu, infix, irpn, itrue, mat, maxelm, maxfix, & maxgauss, maxmat, maxnod, maxord, maxrpn, nelem, ngauss, nodes, norder, & nx, ny, psi, title, wgauss, x, xgauss ) ! !******************************************************************************* ! !! POST post-processes the solution. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxord integer maxrpn ! character com real dpsi(2,maxord) real error real gu(maxnod) integer i integer ielem integer ierror integer ifrm character ( len = * ) infix integer irpn(maxrpn) character isay integer itrue integer ix integer iy integer mat(maxelm) integer maxmat character ( len = 10 ) namvar integer nel integer nelem integer ngauss integer nodes(maxord,maxelm) integer norder integer nx integer ny real psi(maxord) logical s_eqi character ( len = 70 ) title real ucomp real utrue real value real wgauss(maxgauss) real x(2,maxnod) real xgauss(2,maxgauss) real xval real yval ! ifrm = 0 if ( itrue == 1 ) then call enorms ( dpsi, gu, infix, irpn, mat, maxelm, maxfix, maxgauss, & maxmat, maxnod, maxord, maxrpn, nelem, ngauss, nodes, norder, psi, & wgauss, x, xgauss ) end if ! ! Print out current solution ! write ( *, * ) ' ' write ( *, '(a)' ) title write ( *, * ) ' ' if ( itrue == 0 ) then write(*,1070) else write(*,1080) end if write ( *, * ) ' ' i = 0 do ix = 1, nx do iy = 1, ny if ( mod(ix,2) == 0 .and. mod(iy,2) == 0 .and. norder == 8 )go to 19 i = i+1 ucomp = gu(i) xval = x(1,i) yval = x(2,i) if ( itrue == 1 ) then com = 'v' namvar = 'x' value = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) com = 'v' namvar = 'y' value = yval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) com = 'e' ifrm = maxmat*3+1 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) utrue = value error = utrue - ucomp end if if ( itrue == 0 ) then write ( *, 1090 ) i, xval, yval, ucomp else write ( *, 1090 ) i, xval, yval, ucomp, utrue, error end if 19 continue end do write ( *, * ) ' ' end do write ( *, * ) ' ' write ( *, * ) 'Enter "Y" to display flux terms' read ( *, '(a)' ) isay if ( .not. s_eqi ( isay, 'Y' ) ) then return end if do nel = 1, nelem ielem = nel write ( *, * ) ' ' write ( *, '(a,i4)' ) ' Element ', nel write ( *, * ) ' ' write ( *, 1110 ) call flux ( dpsi, gu, ielem, infix, irpn, mat(nel), maxelm, & maxfix, maxgauss, maxnod, maxord, maxrpn, ngauss, nodes, & norder, psi, x, xgauss ) end do return 1070 format('node x',11x,'y',8x,'u computed') 1080 format('node x',11x,'y',8x,'u computed',6x,'u true',8x,'error') 1090 format(i4,2g12.4,3g14.6) 1110 format(5x,'x',19x,'y',13x,'-k*du/dx',11x,'-k*du/dy') end subroutine prep ( infix, irpn, mat, max_ebc, max_nbc, maxelm, maxfix, & maxgauss, maxmat, maxnod, maxord, maxpt, maxrpn, num_ebc, num_nbc, nel_nbc, & nelem, ngauss, nmat, nnode, node_ebc, nodes, norder, npoint, npt, side_nbc, & nx, ny, title, val_ebc, val_nbc, vpt, x, wgauss, xgauss ) ! !******************************************************************************* ! !! PREP preprocesses the problem. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAX_EBC, the maximum number of essential boundary ! conditions. ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NODE_EBC(MAX_EBC), a list of the nodes at which essential ! boundary conditions are applied. ! ! Output, integer NORDER, the element order. ! ! Input, integer NUM_EBC, the number of essential boundary conditions. ! ! Input, real VAL_EBC(MAX_EBC), the value of the essential boundary ! conditions. ! implicit none ! integer max_ebc integer max_nbc integer maxelm integer maxfix integer maxgauss integer maxnod integer maxord integer maxpt integer maxrpn ! character ( len = * ) infix integer irpn(maxrpn) integer mat(maxelm) integer maxmat integer num_ebc integer num_nbc integer nel_nbc(max_nbc) integer nelem integer ngauss integer nmat integer nnode integer node_ebc(max_ebc) integer nodes(maxord,maxelm) integer norder integer npoint integer npt(maxpt) integer side_nbc(max_nbc) integer nx integer ny character ( len = 70 ) title real val_ebc(max_ebc) real val_nbc(2,max_nbc) real vpt(maxpt) real wgauss(maxgauss) real x(2,maxnod) real xgauss(2,maxgauss) ! ! 1: Problem title ! write ( *, * ) 'Enter problem title.' read ( *, '(a)' ) title ! ! 2: Get the materials. ! write ( *, * ) 'Materials' call material_read ( infix, irpn, maxfix, maxmat, maxrpn, nmat ) ! ! 3: Get the element type. ! write ( *, * ) 'Elements' call element_read ( ngauss, norder ) call setint ( maxgauss, ngauss, norder, wgauss, xgauss ) ! ! 4: Get X, Y coordinates of nodes ! write ( *, * ) 'Node coordinates' call node_read ( maxnod, nnode, norder, nx, ny, x ) call selem ( mat, maxelm, nelem, nmat, nodes, norder, nx, ny ) ! ! 5: Define essential boundary conditions ! write ( *, * ) ' ' write ( *, * ) 'Essential boundary conditions' call ebc_read ( infix, irpn, max_ebc, maxfix, maxmat, maxnod, maxrpn, & num_ebc, nnode, node_ebc, norder, nx, ny, val_ebc, x ) ! ! 6: Define natural boundary conditions ! write ( *, * ) ' ' write ( *, * ) 'Natural boundary conditions' call nbc_read ( infix, irpn, max_nbc, maxelm, maxfix, maxmat, maxnod, & maxord, maxrpn, num_nbc, nel_nbc, nelem, nodes, norder, side_nbc, nx, & ny, val_nbc, x ) ! ! 7: Get point loads ! write ( *, * ) ' ' write ( *, * ) 'Point loads' npoint = 0 call point_read ( maxpt, nnode, npoint, npt, vpt ) return end subroutine proc ( dpsi, ef, ek, gf, gk, gu, ierror, infix, ipivot, irpn, mat, & max_ebc, max_nbc, maxelm, maxfix, maxgauss, maxnod, maxord, maxpt, maxban, & maxrpn, ml, mu, num_ebc, num_nbc, nel_nbc, nelem, ngauss, nnode, node_ebc, & nodes, norder, npoint, npt, nrow, side_nbc, psi, val_ebc, val_nbc, vpt, & wgauss, x, xgauss ) ! !******************************************************************************* ! !! PROC processes the problem, that is, it computes the solution. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxban integer max_ebc integer max_nbc integer maxelm integer maxfix integer maxgauss integer maxnod integer maxord integer maxpt integer maxrpn ! real dpsi(2,maxord) real ef(maxord) real ek(maxord,maxord) real gf(maxnod) real gk(maxban,maxnod) real gu(maxnod) integer i integer ierror integer ihave integer ineed character ( len = * ) infix integer info integer ipivot(maxnod) integer irhs integer irpn(maxrpn) character isay integer isolve integer j integer job integer mat(maxelm) integer ml integer mu integer num_ebc integer num_nbc integer nel_nbc(max_nbc) integer nelem integer ngauss integer nmax integer nnode integer node_ebc(max_ebc) integer nodes(maxord,maxelm) integer norder integer npoint integer npt(maxpt) integer nrow real psi(maxord) logical s_eqi integer side_nbc(max_nbc) real tol real val_ebc(max_ebc) real val_nbc(2,max_nbc) real vpt(maxpt) real wgauss(maxgauss) real x(2,maxnod) real xgauss(2,maxgauss) ! ! Compute the bandwidth of the stiffness matrix. ! call gk_bandwidth ( maxelm, maxord, ml, mu, nelem, nodes, norder ) write ( *, * ) ' ' write ( *, * ) 'The stiffness matrix bandwidths are:' write ( *, * ) ' Lower bandwidth = ', ml write ( *, * ) ' Upper bandwidth = ', mu write ( *, * ) ' ' write ( *, * ) 'Choose a solution method:' write ( *, * ) ' 0: direct solution, band storage' write ( *, * ) ' 1: Gauss-Seidel solution, band storage' read ( *, * ) isolve nmax = 0 tol = 0.0E+00 if ( isolve /= 0 ) then write ( *, * ) 'Enter maximum number of iterations?' read ( *, * ) nmax write ( *, * ) 'Enter convergence tolerance' read ( *, * ) tol end if write ( *, * ) ' ' write ( *, * ) 'Setup of system GK * U = GF.' ! ! Compare storage requirements ! write ( *, 3014 ) nnode * nnode write ( *, 3015 ) nnode * ( 2 * ml + mu + 1 ) ! ! Set NROW ! nrow = 2*ml+mu+1 ineed = nnode*nrow ihave = maxnod*maxban write ( *, 1040 ) ineed write ( *, 1050 ) ihave if ( ineed > ihave ) then write ( *, '(a)' ) 'Not enough matrix storage to solve problem' ierror = 1 return end if ! ! Zero out the matrix and right hand side. ! gf(1:maxnod) = 0.0E+00 gk(1:maxban,1:maxnod) = 0.0E+00 ! ! Set global matrix and right hand side ! call setgkf ( ef, ek, gf, gk, ierror, infix, irpn, mat, maxelm, maxfix, & maxgauss, maxnod, maxord, maxrpn, ml, mu, nelem, ngauss, nnode, nodes, & norder, nrow, wgauss, x, xgauss ) if ( ierror /= 0 ) then return end if ! ! Modify the right hand side for point loads. ! call point_apply ( gf, maxnod, maxpt, npoint, npt, vpt ) ! ! Modify the matrix and right hand side for essential boundary conditions. ! call ebc_apply ( gf, gk, max_ebc, maxnod, ml, mu, num_ebc, nnode, node_ebc, & nrow, val_ebc ) ! ! Modify the matrix and right hand side for natural boundary conditions. ! call nbc_apply ( dpsi, ef, ek, gf, gk, max_nbc, maxelm, maxnod, maxord, ml, & mu, num_nbc, nel_nbc, nnode, nodes, norder, nrow, side_nbc, psi, val_nbc, x ) ! ! Print matrix? ! write ( *, '(a)' ) 'Enter Y to print global matrix GK.' read ( *, '(a)' ) isay if ( s_eqi ( isay, 'y' ) ) then irhs = 1 call gkf_write ( gk, ml, mu, nnode, nrow, irhs, gf ) end if ! ! Solve linear system ! if ( isolve == 0 ) then call sgb_fa ( gk, nrow, nnode, ml, mu, ipivot, info ) if ( info /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PROC - Fatal error!' write ( *, '(a)' ) ' Matrix GK is singular' ierror = 1 return end if job = 0 gu(1:nnode) = gf(1:nnode) call sgb_sl ( gk, nrow, nnode, ml, mu, ipivot, gu, job ) else call seidel ( gf, gk, gu, ml, mu, nmax, nnode, nrow, tol ) end if return 1040 format('gk matrix storage requires ',i6,' entries.') 1050 format(i6,' entries are available.') 3014 format('full storage requires ',i6,' entries') 3015 format('band storage requires ',i6,' entries.') 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 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 ( *, * ) ' ' write ( *, * ) 'RPNSET - Error!' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'RPNSET - Error!' write ( *, * ) ' 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, irtol, 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 iarg3 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 irtol 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 iarg3 = 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, iarg3, ierror, ifree, iopsym, iprsym, & ipset, ipval, irad, irtol, itemp, maxsym, maxval, nsyms, & sym, symbol, valsym, value ) 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 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 ( string, sub ) ! !******************************************************************************* ! !! S_INDEXI is a case-insensitive INDEX function. ! ! ! Discussion: ! ! The runction returns the location in 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 STRING is of length 80, and SUB is of ! length 80, then if STRING = '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 STRING 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 = * ) STRING, 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 ! STRING. Otherwise STRING(S_INDEXI:S_INDEXI+LENS-1) = SUB, ! where LENS is the length of SUB, and is the first place ! this happens. However, note that S_INDEXI ignores case, ! unlike the standard FORTRAN INDEX function. ! implicit none ! integer i integer llen1 integer llen2 logical s_eqi integer s_indexi character ( len = * ) string character ( len = * ) sub ! s_indexi = 0 llen1 = len_trim ( string ) llen2 = len_trim ( sub ) ! ! In case STRING or SUB is blanks, use LEN. ! if ( llen1 == 0 ) then llen1 = len ( string ) 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 ( string(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 scopy ( n, sx, incx, sy, incy ) ! !******************************************************************************* ! !! SCOPY copies a vector, X, to a vector, Y. ! ! ! Parameters: ! ! Input, integer N, the number of entries to copy. ! ! Input, real SX(*), the vector to be copied. ! ! Input, integer INCX, the increment between successive ! entries of SX. ! ! Output, real SY(*), a copy of SX. ! ! Input, integer INCY, the increment between successive ! entries of SY. ! implicit none ! integer i integer incx integer incy integer ix integer iy integer n real sx(*) real sy(*) ! if ( n<=0 ) return ix = 1 iy = 1 if ( incx<0)ix = (-n+1)*incx + 1 if ( incy<0)iy = (-n+1)*incy + 1 do i = 1, n sy(iy) = sx(ix) ix = ix + incx iy = iy + incy end do return end subroutine seidel ( gf, gk, gu, ml, mu, nmax, nnode, nrow, tol ) ! !******************************************************************************* ! !! SEIDEL carries out the Gauss-Seidel iterative linear system solution. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer nnode integer nrow ! real entry real gf(nnode) real gk(nrow,nnode) real gu(nnode) integer i integer irow integer j integer jcol integer jhi integer jlo integer ml integer mu integer nmax integer nstep real res real resnew real reszer real tol ! ! Make a rough estimate of the solution. ! do i = 1, nnode call getmat ( gk, ml, mu, nnode, nrow, entry, i, i ) gu(i) = gf(i) / entry end do nstep = 0 10 continue ! ! Compute the Residual: ! ! RES(I) = GF(I) - GK(I,J) * GU(J) ! resnew = 0.0E+00 do i = 1, nnode jlo = max(1,i-ml) jhi = min(nnode,i+mu) res = gf(i) irow = i do j = jlo, jhi jcol = j call getmat ( gk, ml, mu, nnode, nrow, entry, irow, jcol ) res = res - entry * gu(j) end do resnew = max ( resnew, abs(res) ) end do if ( nstep == 0 ) then reszer = resnew end if write(*,1000) resnew if ( resnew < tol * (reszer+1.0E+00) ) then write ( *, '(a)' ) 'Convergence achieved.' go to 60 end if nstep = nstep+1 if ( nstep > nmax ) then write ( *, '(a)' ) 'Maximum number of steps reached' go to 60 end if do i = 1,nnode irow = i jlo = max(1,i-ml) jhi = min(nnode,i+mu) gu(i) = gf(i) do j = jlo,jhi if ( i/=j ) then jcol = j call getmat ( gk, ml, mu,nnode,nrow,entry,irow,jcol) gu(i) = gu(i) - entry * gu(j) end if end do jcol = i call getmat(gk, ml, mu,nnode,nrow,entry,irow,jcol) gu(i) = gu(i) / entry end do go to 10 60 continue return 1000 format('residual = ',g14.6) end subroutine selem ( mat, maxelm, nelem, nmat, nodes, norder, nx, ny ) ! !******************************************************************************* ! !! SELEM sets the node numbers and material types for each element. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm ! integer i integer ielem integer j integer mat(maxelm) integer mval integer nelem integer nmat integer nodes(9,maxelm) integer norder integer nx integer ny integer nymid ! ! 3-Node triangles ! if ( norder == 3 ) then nelem = 2*(nx-1)*(ny-1) ielem = 0 do i = 1,nx-1 do j = 1,ny-1 ielem = ielem+1 nodes(1,ielem) = (i-1)*ny+j nodes(2,ielem) = i*ny+j nodes(3,ielem) = (i-1)*ny+j+1 ielem = ielem+1 nodes(2,ielem) = (i-1)*ny+j+1 nodes(1,ielem) = i*ny+j+1 nodes(3,ielem) = i*ny+j end do end do ! ! 4-Node quadrilaterals ! else if ( norder == 4 ) then nelem = (nx-1)*(ny-1) ielem = 0 do i = 1,nx-1 do j = 1,ny-1 ielem = ielem+1 nodes(1,ielem) = (i-1)*ny+j nodes(2,ielem) = i*ny+j nodes(3,ielem) = i*ny+j+1 nodes(4,ielem) = (i-1)*ny+j+1 end do end do ! ! 6-Node triangles ! else if ( norder == 6 ) then nelem = (nx-1)*(ny-1)/2 ielem = 0 do i = 1,nx-2,2 do j = 1,ny-2,2 ielem = ielem+1 nodes(1,ielem) = (i-1)*ny+j nodes(2,ielem) = (i+1)*ny+j nodes(3,ielem) = (i-1)*ny+j+2 nodes(4,ielem) = i*ny+j nodes(5,ielem) = i*ny+j+1 nodes(6,ielem) = (i-1)*ny+j+1 ielem = ielem+1 nodes(1,ielem) = (i+1)*ny+j+2 nodes(2,ielem) = (i-1)*ny+j+2 nodes(3,ielem) = (i+1)*ny+j nodes(4,ielem) = i*ny+j+2 nodes(5,ielem) = i*ny+j+1 nodes(6,ielem) = (i+1)*ny+j+1 end do end do ! ! 8-Node quadrilaterals ! else if ( norder == 8 ) then nelem = (nx-1)*(ny-1)/4 ielem = 0 nymid = (ny+1)/2 do i = 1, nx-2, 2 do j = 1, ny-2, 2 ielem = ielem+1 nodes(1,ielem) = ((i-1)/2)*(nymid+ny)+j nodes(2,ielem) = nodes(1,ielem)+ny+nymid nodes(3,ielem) = nodes(2,ielem)+2 nodes(4,ielem) = nodes(1,ielem)+2 nodes(5,ielem) = nodes(1,ielem)+ny-((j-1)/2) nodes(6,ielem) = nodes(2,ielem)+1 nodes(7,ielem) = nodes(1,ielem)+ny+1-((j-1)/2) nodes(8,ielem) = nodes(1,ielem)+1 end do end do ! ! 9-Node quadrilaterals ! else if ( norder == 9 ) then nelem = (nx-1)*(ny-1)/4 ielem = 0 do i = 1, nx-2, 2 do j = 1, ny-2, 2 ielem = ielem+1 nodes(1,ielem) = (i-1)*ny+j nodes(2,ielem) = (i+1)*ny+j nodes(3,ielem) = (i+1)*ny+j+2 nodes(4,ielem) = (i-1)*ny+j+2 nodes(5,ielem) = i*ny+j nodes(6,ielem) = (i+1)*ny+j+1 nodes(7,ielem) = i*ny+j+2 nodes(8,ielem) = (i-1)*ny+j+1 nodes(9,ielem) = i*ny+j+1 end do end do else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SELEM - Fatal error!' write ( *, * ) ' Illegal element order = ', norder stop end if ! ! Material specification ! if ( nmat > 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Specify a material type for each element.' do i = 1, nelem write ( *, * ) 'Enter the material number for element ', i read ( *, * ) mval mat(i) = mval end do else do i = 1, nelem mat(i) = 1 end do end if return end subroutine setgkf ( ef, ek, gf, gk, ierror, infix, irpn, mat, maxelm, maxfix, & maxgauss, maxnod, maxord, maxrpn, ml, mu, nelem, ngauss, nnode, nodes, & norder, nrow, wgauss, x, xgauss ) ! !******************************************************************************* ! !! SETGKF sets up the global stiffness matrix and right hand side. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxelm integer maxfix integer maxgauss integer maxnod integer maxord integer maxrpn integer nnode integer norder integer nrow ! real ef(norder) real ek(norder,norder) real entry real gf(maxnod) real gk(nrow,nnode) integer ielem integer ieqn integer ierror integer imax integer in character ( len = * ) infix integer inode integer irpn(maxrpn) integer ishow integer j integer jnode integer jvar integer mat(maxelm) integer ml integer mu integer nel integer nelem integer ngauss integer nodes(maxord,maxelm) real wgauss(maxgauss) real x(2,maxnod) real xa(2,9) real xgauss(2,maxgauss) ! ! Want to see an element matrix? ! ierror = 0 ishow = 0 imax = nelem+1 do nel = 1, nelem ielem = nel if ( ishow == 0 .and. nel /= 1 )go to 44 if ( ishow /= imax ) then if ( ishow < nel ) then write ( *, '(a)' ) '0 to go on, or next element matrix to display.' read ( *, * ) ishow end if end if 44 continue do j = 1,norder in = nodes(j,nel) xa(1,j) = x(1,in) xa(2,j) = x(2,in) end do call ekf_set ( ef, ek, ierror, infix, irpn, mat(nel), & maxfix, maxgauss, maxrpn, norder, ngauss, wgauss, xa, xgauss ) if ( ierror/=0)return if ( ishow == nel .or. ishow==imax ) then call ekf_write ( ef, ek, maxelm, maxord, ielem, nodes, norder ) end if do inode = 1, norder ieqn = nodes(inode,nel) gf(ieqn) = gf(ieqn) + ef(inode) do jnode = 1, norder jvar = nodes(jnode,nel) entry = ek(inode,jnode) call addmat ( gk, ml, mu, nnode, nrow, entry, ieqn, jvar ) end do end do end do return end subroutine setint ( maxgauss, ngauss, norder, wgauss, xgauss ) ! !******************************************************************************* ! !! SETINT sets the quadrature points and weights. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NORDER, the element order. ! implicit none ! integer maxgauss ! integer ngauss integer norder real wgauss(maxgauss) real xgauss(2,maxgauss) if ( ngauss == 3 ) then wgauss(1) = 1.0E+00 / 3.0E+00 xgauss(1,1) = 0.0E+00 xgauss(2,1) = 0.5E+00 wgauss(2) = 1.0E+00 / 3.0E+00 xgauss(1,2) = 0.5E+00 xgauss(2,2) = 0.0E+00 wgauss(3) = 1.0E+00 / 3.0E+00 xgauss(1,3) = 0.5E+00 xgauss(2,3) = 0.5E+00 else if ( ngauss == 4.and. & (norder == 4 .or. norder==8.or.norder==9) ) then wgauss(1) = 0.25E+00 xgauss(1,1) = - 1.0E+00 / sqrt(3.0E+00) xgauss(2,1) = - 1.0E+00 / sqrt(3.0E+00) wgauss(2) = 0.25E+00 xgauss(1,2) = 1.0E+00 / sqrt(3.0E+00) xgauss(2,2) = - 1.0E+00 / sqrt(3.0E+00) wgauss(3) = 0.25E+00 xgauss(1,3) = - 1.0E+00 / sqrt(3.0E+00) xgauss(2,3) = 1.0E+00 / sqrt(3.0E+00) wgauss(4) = 0.25E+00 xgauss(1,4) = 1.0E+00 / sqrt(3.0E+00) xgauss(2,4) = 1.0E+00 / sqrt(3.0E+00) else if ( ngauss == 4.and.(norder==3 .or. norder==6) ) then wgauss(1) = -27.0E+00 / 96.0E+00 xgauss(1,1) = 1.0E+00 / 3.0E+00 xgauss(2,1) = 1.0E+00 / 3.0E+00 wgauss(2) = 25.0E+00 / 96.0E+00 xgauss(1,2) = 2.0E+00 / 15.0E+00 xgauss(2,2) = 11.0E+00 / 15.0E+00 wgauss(3) = 25.0E+00 / 96.0E+00 xgauss(1,3) = 2.0E+00 / 15.0E+00 xgauss(2,3) = 2.0E+00 / 15.0E+00 wgauss(4) = 25.0E+00 / 96.0E+00 xgauss(1,4) = 11.0E+00 / 15.0E+00 xgauss(2,4) = 2.0E+00 / 15.0E+00 else if ( ngauss == 7 ) then wgauss(1) = 3.0E+00/120.0E+00 xgauss(1,1) = 0.0E+00 xgauss(2,1) = 0.0E+00 wgauss(2) = 3.0E+00/120.0E+00 xgauss(1,2) = 1.0E+00 xgauss(2,2) = 0.0E+00 wgauss(3) = 3.0E+00/120.0E+00 xgauss(1,3) = 0.0E+00 xgauss(2,3) = 1.0E+00 wgauss(4) = 8.0E+00/120.0E+00 xgauss(1,4) = 0.5E+00 xgauss(2,4) = 0.0E+00 wgauss(5) = 8.0E+00/120.0E+00 xgauss(1,5) = 0.5E+00 xgauss(2,5) = 0.5E+00 wgauss(6) = 8.0E+00/120.0E+00 xgauss(1,6) = 0.0E+00 xgauss(2,6) = 0.5E+00 wgauss(7) = 27.0E+00/120.0E+00 xgauss(1,7) = 1.0E+00/3.0E+00 xgauss(2,7) = 1.0E+00/3.0E+00 else if ( ngauss == 9 ) then wgauss(1) = 25.0E+00/81.0E+00 xgauss(1,1) = -sqrt(0.6E+00) xgauss(2,1) = -sqrt(0.6E+00) wgauss(2) = 40.0E+00/81.0E+00 xgauss(1,2) = 0.0E+00 xgauss(2,2) = -sqrt(0.6E+00) wgauss(3) = 25.0E+00/81.0E+00 xgauss(1,3) = sqrt(0.6E+00) xgauss(2,3) = -sqrt(0.6E+00) wgauss(4) = 40.0E+00/81.0E+00 xgauss(1,4) = -sqrt(0.6E+00) xgauss(2,4) = 0.0E+00 wgauss(5) = 64.0E+00/81.0E+00 xgauss(1,5) = 0.0E+00 xgauss(2,5) = 0.0E+00 wgauss(6) = 40.0E+00/81.0E+00 xgauss(1,6) = sqrt(0.6E+00) xgauss(2,6) = 0.0E+00 wgauss(7) = 25.0E+00/81.0E+00 xgauss(1,7) = -sqrt(0.6E+00) xgauss(2,7) = sqrt(0.6E+00) wgauss(8) = 40.0E+00/81.0E+00 xgauss(1,8) = 0.0E+00 xgauss(2,8) = sqrt(0.6E+00) wgauss(9) = 25.0E+00/81.0E+00 xgauss(1,9) = sqrt(0.6E+00) xgauss(2,9) = sqrt(0.6E+00) end if return end subroutine setmat ( a, ml, mu, neqn, nrow, entry, i, j ) ! !******************************************************************************* ! !! SETMAT sets A(I,J) = ENTRY. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! implicit none ! integer neqn integer nrow ! real a(nrow,neqn) real entry integer i integer j integer ml integer mu ! ! Check range of i and j ! if ( i <= 0 .or. i > neqn ) go to 60 if(j<=0.or.j>neqn)go to 60 ! ! Banded LINPACK storage ! if((j-i)>mu)go to 50 if((i-j)>ml)go to 50 a(i-j+ml+mu+1,j) = entry return ! ! Indices i or j out of bandwidth ! 50 continue if(entry==0.0)return write ( *, * ) ' ' write ( *, * ) 'SETMAT - Fatal error!' write(*,1000)i,j,ml,mu,entry return ! ! indices i or j out of range ! 60 continue write ( *, * ) ' ' write ( *, * ) 'SETMAT - Fatal error!' write(*,1010)i,j,neqn,entry return 1000 format(' illegal i,j = ',2i6,' ml,mu=',2i6,' entry=',g14.6) 1010 format(' illegal i,j=',2i6,' neqn=',i6,' entry=',g14.6) end subroutine setnbc ( ihi, ilo, inc, infix, irpn, max_nbc, maxelm, maxfix, & maxmat, maxnod, maxord, maxrpn, num_nbc, nel_nbc, nodes, none, ns, & side_nbc, ntwo, val_nbc, x ) ! !******************************************************************************* ! !! SETNBC sets natural boundary conditions. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! implicit none ! integer max_nbc integer maxelm integer maxfix integer maxnod integer maxord integer maxrpn ! character com integer i integer ierror integer ifrm integer ihi integer ilo integer inc character ( len = * ) infix integer irpn(maxrpn) integer maxmat integer mone integer mtwo character ( len = 10 ) namvar integer num_nbc integer nel_nbc(max_nbc) integer nodes(maxord,maxelm) integer none integer ns integer side_nbc(max_nbc) integer ntwo real v1val real v2val real value real val_nbc(2,max_nbc) real x(2,maxnod) real xval real yval ! ifrm = 0 do i = ilo, ihi, inc mone = nodes(none,i) mtwo = nodes(ntwo,i) xval = 0.5 * (x(1,mone)+x(1,mtwo)) com = 'v' namvar = 'x' value = xval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) yval = 0.5 * (x(2,mone)+x(2,mtwo)) com = 'v' namvar = 'y' value = yval call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) com = 'e' ifrm = maxmat*3+4 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) v1val = value com = 'e' ifrm = maxmat*3+5 call comrpn ( com, ierror, ifrm, infix, irpn, maxrpn, & namvar, 1, 1, value ) v2val = value num_nbc = num_nbc+1 nel_nbc(num_nbc) = i side_nbc(num_nbc) = ns val_nbc(1,num_nbc) = v1val val_nbc(2,num_nbc) = v2val end do 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 ( *, * ) ' ' write ( *, * ) 'SGB_CHECK - Illegal LDA = ', lda end if if ( m < 1 ) then ierror = ierror + 2 write ( *, * ) ' ' write ( *, * ) 'SGB_CHECK - Illegal M = ', m end if if ( ml < 0 .or. ml > min ( m, n ) - 1 ) then ierror = ierror + 4 write ( *, * ) ' ' write ( *, * ) 'SGB_CHECK - Illegal ML = ', ml end if if ( mu < 0 .or. mu > min ( m, n ) - 1 ) then ierror = ierror + 8 write ( *, * ) ' ' write ( *, * ) 'SGB_CHECK - Illegal MU = ', mu end if if ( n < 1 ) then ierror = ierror + 16 write ( *, * ) ' ' write ( *, * ) '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 ! ! Author: ! ! John Burkardt ! ! 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 ( *, * ) ' ' write ( *, * ) 'SGB_FA - Fatal error!' write ( *, * ) ' 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 ! ! Author: ! ! John Burkardt ! ! 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, xsi, eta ) ! !******************************************************************************* ! !! SHAPE evaluates the shape functions on the reference element. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real DPSI(2,NORDER), DPSI(1,J) is the value of dPSI(J)/dX. ! DPSI(2,J) is the value of dPSI(J)/dY. ! ! Input, integer NORDER, the order of the element. ! ! Output, real PSI(NORDER), the value of the shape functions at (XSI,ETA). ! ! Input, real XSI, ETA, the reference element coordinates of the point. ! implicit none ! integer norder ! real dpsi(2,norder) real eta real psi(norder) real xsi ! ! 3 node linear triangles. ! if ( norder == 3 ) then psi(1) = ( 1.0E+00 - xsi - eta ) psi(2) = xsi psi(3) = eta dpsi(1,1) = -1.0E+00 dpsi(2,1) = -1.0E+00 dpsi(1,2) = 1.0E+00 dpsi(2,2) = 0.0E+00 dpsi(1,3) = 0.0E+00 dpsi(2,3) = 1.0E+00 ! ! 4 node linear quadrilaterals. ! else if ( norder == 4 ) then psi(1) = 0.25E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 - eta ) psi(2) = 0.25E+00 * ( 1.0E+00 + xsi ) * ( 1.0E+00 - eta ) psi(3) = 0.25E+00 * ( 1.0E+00 + xsi ) * ( 1.0E+00 + eta ) psi(4) = 0.25E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 + eta ) dpsi(1,1) = -0.25E+00 * ( 1.0E+00 - eta ) dpsi(2,1) = -0.25E+00 * ( 1.0E+00 - xsi ) dpsi(1,2) = 0.25E+00 * ( 1.0E+00 - eta ) dpsi(2,2) = -0.25E+00 * ( 1.0E+00 + xsi ) dpsi(1,3) = 0.25E+00 * ( 1.0E+00 + eta ) dpsi(2,3) = 0.25E+00 * ( 1.0E+00 + xsi ) dpsi(1,4) = -0.25E+00 * ( 1.0E+00 + eta ) dpsi(2,4) = 0.25E+00 * ( 1.0E+00 - xsi ) ! ! 6 node quadratic triangles. ! else if ( norder == 6 ) then psi(1) = 2.0E+00*(1.0-xsi-eta)*(0.5-xsi-eta) psi(2) = 2.0E+00*xsi*(xsi-0.5) psi(3) = 2.0E+00*eta*(eta-0.5) psi(4) = 4.0E+00*(1.0-xsi-eta)*xsi psi(5) = 4.0E+00*xsi*eta psi(6) = 4.0E+00*eta*(1.0-xsi-eta) dpsi(1,1) = -3.0+4.0*(xsi+eta) dpsi(2,1) = -3.0+4.0*(xsi+eta) dpsi(1,2) = 4.0*xsi-1.0E+00 dpsi(2,2) = 0.0E+00 dpsi(1,3) = 0.0E+00 dpsi(2,3) = 4.0*eta-1.0E+00 dpsi(1,4) = 4.0-8.0*xsi-4.0*eta dpsi(2,4) = -4.0*xsi dpsi(1,5) = 4.0*eta dpsi(2,5) = 4.0*xsi dpsi(1,6) = -4.0*eta dpsi(2,6) = 4.0-4.0*xsi-8.0*eta ! ! 8 node serendipity quadrilaterals. ! else if ( norder == 8 ) then psi(1) = -0.25*(1.0-xsi)*(1.0-eta)*(1.0+xsi+eta) psi(2) = -0.25*(1.0+xsi)*(1.0-eta)*(1.0-xsi+eta) psi(3) = -0.25*(1.0+xsi)*(1.0+eta)*(1.0-xsi-eta) psi(4) = -0.25*(1.0-xsi)*(1.0+eta)*(1.0+xsi-eta) psi(5) = 0.5*(1.0-xsi)*(1.0+xsi)*(1.0-eta) psi(6) = 0.5*(1.0+xsi)*(1.0+eta)*(1.0-eta) psi(7) = 0.5*(1.0-xsi)*(1.0+xsi)*(1.0+eta) psi(8) = 0.5*(1.0-xsi)*(1.0+eta)*(1.0-eta) dpsi(1,1) = 0.25*(1.0-eta)*(2.0*xsi+eta) dpsi(1,2) = 0.25*(1.0-eta)*(2.0*xsi-eta) dpsi(1,3) = 0.25*(1.0+eta)*(2.0*xsi+eta) dpsi(1,4) = 0.25*(1.0+eta)*(2.0*xsi-eta) dpsi(1,5) = -xsi*(1.0-eta) dpsi(1,6) = 0.5*(1.0+eta)*(1.0-eta) dpsi(1,7) = -xsi*(1.0+eta) dpsi(1,8) = -0.5*(1.0+eta)*(1.0-eta) dpsi(2,1) = 0.25*(1.0-xsi)*(2.0*eta+xsi) dpsi(2,2) = 0.25*(1.0+xsi)*(2.0*eta-xsi) dpsi(2,3) = 0.25*(1.0+xsi)*(2.0*eta+xsi) dpsi(2,4) = 0.25*(1.0-xsi)*(2.0*eta-xsi) dpsi(2,5) = -0.5*(1.0-xsi)*(1.0+xsi) dpsi(2,6) = -eta*(1.0+xsi) dpsi(2,7) = 0.5*(1.0-xsi)*(1.0+xsi) dpsi(2,8) = -eta*(1.0-xsi) ! ! 9 node quadratic quadrilaterals. ! else if ( norder == 9 ) then psi(1) = 0.25*xsi*(xsi-1.0)*eta*(eta-1.0) psi(2) = 0.25*xsi*(xsi+1.0)*eta*(eta-1.0) psi(3) = 0.25*xsi*(xsi+1.0)*eta*(eta+1.0) psi(4) = 0.25*xsi*(xsi-1.0)*eta*(eta+1.0) psi(5) = 0.5*(xsi+1.0)*(1.0-xsi)*xsi*(xsi-1.0) psi(6) = 0.5*xsi*(xsi+1.0)*(1.0+eta)*(1.0-eta) psi(7) = 0.5*(xsi+1.0)*(1.0-xsi)*eta*(eta+1.0) psi(8) = 0.5*xsi*(xsi-1.0)*(eta+1.0)*(1.0-eta) psi(9) = (xsi+1.0)*(1.0-xsi)*(eta+1.0)*(1.0-eta) dpsi(1,1) = (0.5*xsi-0.25)*eta*(eta-1.0) dpsi(2,1) = (0.5*eta-0.25)*xsi*(xsi-1.0) dpsi(1,2) = (0.5*xsi+0.25)*eta*(eta-1.0) dpsi(2,2) = (0.5*eta-0.25)*xsi*(xsi+1.0) dpsi(1,3) = (0.5*xsi+0.25)*eta*(eta+1.0) dpsi(2,3) = (0.5*eta+0.25)*xsi*(xsi+1.0) dpsi(1,4) = (0.5*xsi-0.25)*eta*(eta+1.0) dpsi(2,4) = (0.5*eta+0.25)*xsi*(xsi-1.0) dpsi(1,5) = -xsi*eta*(eta-1.0) dpsi(2,5) = -(eta-0.5)*(xsi-1.0)*(xsi+1.0) dpsi(1,6) = -(xsi+0.5)*(eta-1.0)*(eta+1.0) dpsi(2,6) = -eta*xsi*(xsi+1.0) dpsi(1,7) = -xsi*eta*(eta+1.0) dpsi(2,7) = -(eta+0.5)*(xsi-1.0)*(xsi+1.0) dpsi(1,8) = -(xsi-0.5)*(eta-1.0)*(eta+1.0) dpsi(2,8) = -eta*xsi*(xsi-1.0) dpsi(1,9) = 2.0*xsi*(eta-1.0)*(eta+1.0) dpsi(2,9) = 2.0*eta*(xsi-1.0)*(xsi+1.0) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SHAPE - Fatal error!' write ( *, '(a,i6)' ) ' Illegal element order = ', norder stop end if 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 ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'SYMADD - Error!' write ( *, * ) ' There is not enough memory for ', trim ( namvr ) write ( *, * ) ' 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 ( *, * ) ' 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: ! ! 17 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IFREE, the index of the next free memory ! location in VALSYM. ! ! Input, 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 IOUNIT(4). ! IOUNIT(1) is the FORTRAN input unit. ! IOUNIT(2) is the standard output unit, while IOUNIT(3) and ! IOUNIT(4), if nonzero, are auxilliary output units. ! ! 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. ! ! Input, integer MAXFRM, the maximum number of formulas ! that can be handled. ! ! 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. ! ! Input, character ( len = MAXLEN ) 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. ! ! Input, 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 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 = 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 maxrpn integer maxsym integer maxval ! integer i integer ierror integer ifree integer ifrm integer ihi integer ilo integer inc integer intsym(80) integer iopsym(maxsym) integer ipoint integer iprsym(maxsym) integer ipval(maxsym) integer irpn(maxrpn) integer itemp integer j integer jhi integer jinc 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 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) ! if ( s_eqi ( namvr, 'DEBUG' ) ) then nsyms = max ( nsym, nsyms ) write ( *, * ) ' Symbol Loc Prior Args Nrow Ncol' write ( *, * ) ' ' do i = 1, nsyms write ( *, '(a20,1x,5i6)' ) symbol(i), ipval(i), iprsym(i), & iopsym(i) end do write ( *, * ) 'Next free memory location at ', ifree write ( *, * ) ' ' write ( *, * ) 'Symbol Value' do i = 1, nsyms jlo = ipval(i) jhi = jlo jinc = 3 do klo = jlo, jhi, jinc khi = min(klo+jinc-1,jhi) write ( *, '(i4,3g20.8)' ) i, (valsym(k), k = klo, khi ) end do end do end if if ( s_eqi ( namvr, 'DEBUG' ) .or. s_eqi ( namvr, 'ALL' ) ) then write ( *, * ) ' ' write ( *, * ) 'Predefined symbols:' write ( *, * ) ' ' inc = 3 ihi = ( ( nsymp - 1 ) / inc ) + 1 do i = 1, ihi jhi = i + (inc-1)*ihi jhi = min ( jhi, nsymp ) write ( *, '(3(a20,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 ( *, * ) ' ' write ( *, * ) '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)' ) 'User RPN formulas:' write ( *, '(a)' ) ' ' do ifrm = 1, maxfrm+1 if ( ifrm == maxfrm+1 ) then if ( nints <= 0 ) then go to 100 end if write ( *, * ) ' ' write ( *, * ) 'Most recent user infix formula:' write ( *, * ) ' ' ilo = 1 ihi = nints else ilo = ( ifrm - 1 ) * 80 + 1 ihi = ilo + 79 end if jlo = 0 output = ' ' 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) lens = len_trim ( symbol(itemp) ) if ( (jlo+lens)>80 ) then write ( *, '(a)' ) output output = ' ' jlo = 0 end if jhi = jlo + lens jlo = jlo + 1 output(jlo:jhi) = sym(1:lens) jlo = jhi + 1 if ( sym == '$' ) then go to 90 end if 80 continue end do 90 continue write ( *, '(a)' ) output 100 continue end do write ( *, * ) ' ' 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 ( *, * ) 'Warning! There is no variable named ' // 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 ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' No variable named ' // 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 ( *, * ) ' ' write ( *, * ) 'SYMBOL_VALUE - Error!' write ( *, * ) ' You are trying to assign a value to' write ( *, * ) ' the name of a function or operator:' write ( *, * ) 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), 'RTOL' ) ) 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 ( *, * ) ' Removed the variable ' // trim ( namvar ) write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'TOKENS:' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'TOKENS - Error!' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) 'TOKENS - Error!' write ( *, * ) ' 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 ( *, * ) ' ' write ( *, * ) '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 ( *, * ) ' ' write ( *, * ) 'TOKENS:' write ( *, * ) ' The undeclared variable ', trim ( namvr ), & ' will be assumed to be a scalar.' go to 20 end subroutine xsitox ( norder, xnode, xsi, x, y ) ! !******************************************************************************* ! !! XSITOX converts coordinates in the master XSI system to the X, Y system. ! ! ! Modified: ! ! 22 May 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NORDER, the order of the element. ! ! Input, real XNODE(2,NORDER), the X and Y coordinates of the nodes ! that define the element. ! ! Input, real XSI(2), the XSI and ETA local coordinates of the node. ! ! Output, real X, Y, the X and Y global coordinates of the node. ! implicit none ! integer norder ! real wne real wnw real wse real wsw real x real xnode(2,norder) real xsi(2) real y ! if ( norder == 3 .or. norder == 6 ) then wsw = 1.0E+00 - xsi(1) - xsi(2) wnw = xsi(2) wse = xsi(1) x = wsw * xnode(1,1) + wse * xnode(1,2) + wnw * xnode(1,3) y = wsw * xnode(2,1) + wse * xnode(2,2) + wnw * xnode(2,3) else wsw = 0.25E+00 * ( 1.0E+00 - xsi(1) ) * ( 1.0E+00 - xsi(2) ) wnw = 0.25E+00 * ( 1.0E+00 - xsi(1) ) * ( 1.0E+00 + xsi(2) ) wse = 0.25E+00 * ( 1.0E+00 + xsi(1) ) * ( 1.0E+00 - xsi(2) ) wne = 0.25E+00 * ( 1.0E+00 + xsi(1) ) * ( 1.0E+00 + xsi(2) ) x = wsw * xnode(1,1) + wse * xnode(1,2) & + wne * xnode(1,3) + wnw * xnode(1,4) y = wsw * xnode(2,1) + wse * xnode(2,2) & + wne * xnode(2,3) + wnw * xnode(2,4) end if return end