program main c*********************************************************************72 c cc fem2d_poisson_rectangle() solves the Poisson equation on a rectangle. c c Discussion: c c The code solves c c -Laplacian U(X,Y) = F(X,Y) c c in a rectangular region in the plane. Homogeneous Dirichlet c boundary conditions are imposed. The code uses continuous c piecewise quadratic basis functions on triangles determined by c a uniform grid of NX by NY points. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Local parameters: c c Local, double precision A(3*IB+1,NUNK), the coefficient matrix. c c Local, double precision AREA(NEL), the area of each element. c c Local, double precision EH1, the H1 seminorm error. c c Local, double precision EL2, the L2 error. c c Local, double precision F(NUNK), the right hand side. c c Local, integer IB, the half-bandwidth of the matrix. c c Local, integer INDX(NP), gives the index of the unknown quantity c associated with the given node. c c Local, integer LDA, the "leading dimension of A". This is the c value used in the first dimension of the declaration statement for A. c The first dimension of A must be at least 3 * IB + 1, where IB is c the half bandwidth of A, but when we start the program, we don't know c what IB is. The number used for LDA is an intelligent guess. Because c LDA is generally different from the actual first dimension of the data, c we need to pass information about both quantities to any subroutines c that want to store into or retrieve from the array. See ASSEM and c DGB_FA and DGB_SL. c c Local, integer NEL, the number of elements. c c Local, integer NNODES, the number of nodes used to form one element. c c Local, integer NODE(NEL,NNODES); NODE(I,J) is the global node index of c the local node J in element I. c c Local, integer NQ, the number of quadrature points used for assembly. c c Local, integer NQE, the number of quadrature points used for error c estimation. c c Local, integer NUNK, the number of unknowns. c c Local, integer NX, the number of points in the X direction. c c Local, integer NY, the number of points in the Y direction. c c Local, double precision WQ(NQ), quadrature weights. c c Local, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c c Local, double precision XL, XR, YB, YT, the X coordinates of c the left and right sides of the rectangle, and the Y coordinates c of the bottom and top of the rectangle. c c Local, double precision XQ(NEL,NQ), YQ(NEL,NQ), the X and Y c coordinates of the quadrature points in each element. c implicit none integer nnodes integer nq integer nqe integer nx integer ny parameter ( nnodes = 6 ) parameter ( nq = 3 ) parameter ( nqe = 13 ) parameter ( nx = 7 ) parameter ( ny = 7 ) integer nel integer np parameter ( nel = ( nx - 1 ) * ( ny - 1 ) * 2 ) parameter ( np = ( 2 * nx - 1 ) * ( 2 * ny - 1 ) ) integer lda parameter ( lda = 6 * ( nx + ny ) + 1 ) double precision a(lda,np) double precision area(nel) double precision eh1 double precision el2 double precision f(np) integer ib integer ierr integer indx(np) integer job integer node(nel,nnodes) integer nunk integer pivot(np) double precision wq(nq) double precision xc(np) double precision xl double precision xq(nel,nq) double precision xr double precision yb double precision yc(np) double precision yq(nel,nq) double precision yt xl = 0.0D+00 xr = 1.0D+00 yb = 0.0D+00 yt = 1.0D+00 call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Fortran77 version' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Solution of the Poisson equation on a' write ( *, '(a)' ) ' unit box in 2 dimensions.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' - Uxx - Uyy = F(x,y) in the box' write ( *, '(a)' ) ' U(x,y) = 0 on the boundary.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The finite element method is used, with' write ( *, '(a)' ) ' piecewise quadratic basis functions on' write ( *, '(a)' ) ' 6 node triangular elements.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The corner nodes of the triangles are' write ( *, '(a)' ) ' generated by an underlying grid whose ' write ( *, '(a)' ) ' dimensions are' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' NX = ', nx write ( *, '(a,i8)' ) ' NY = ', ny c c Set up the geometry information. c call geom ( xc, yc, xq, yq, wq, area, node, indx, nel, & nunk, np, ib, nnodes, xl, xr, yb, yt, nx, ny ) c c Make a picture of the mesh. c if ( nx .le. 10 .and. ny .le. 10 ) then call mesh_eps ( 'fem2d_poisson_rectangle.eps', np, xc, yc, & nnodes, nel, node ) end if c c Allocate space for the coefficient matrix A and right hand side F. c if ( lda < 3 * ib + 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_RECTANGLE - Fatal error!' write ( *, '(a)' ) ' The dimension LDA is too small!' write ( *, '(a,i8)' ) ' The current value is LDA = ', lda write ( *, '(a,i8)' ) ' We need ', 3 * ib + 1 stop end if c c Assemble the coefficient matrix A and the right-hand side F of the c finite element equations. c call assem ( lda, a, f, ib, xc, yc, xq, yq, wq, area, node, & indx, nnodes, nunk, nq, nel, np ) call dgb_print_some ( lda, nunk, nunk, ib, ib, a, 1, 1, 5, 5, & ' Initial 5 x 5 block of matrix A:' ) call r8vec_print_some ( nunk, f, 10, & ' Part of right hand side vector:' ) c c Solve the linear system using a banded solver. c call dgb_fa ( lda, nunk, ib, ib, a, pivot, ierr ) if ( ierr .ne. 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_POISSON_RECTANGLE - Fatal error!' write ( *, '(a)' ) ' DGB_FA returned an error condition.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The linear system was not factored,' write ( *, '(a)' ) ' and the algorithm cannot proceed.' stop end if job = 0 call dgb_sl ( lda, nunk, ib, ib, a, pivot, f, job ) call r8vec_print_some ( nunk, f, 10, & ' Part of the solution vector:' ) c c Calculate error using 13 point quadrature rule. c call eror ( area, node, indx, xc, yc, f, nel, nnodes, & nqe, nunk, np, el2, eh1 ) c c Write an ASCII file that can be read into MATLAB. c call solution_write ( f, indx, np, nunk, & 'fem2d_poisson_rectangle_solution.txt', xc, yc ) c c Terminate. c write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine assem ( lda, a, f, ib, xc, yc, xq, yq, wq, area, & node, indx, nnodes, nunk, nq, nel, np ) c*********************************************************************72 c cc ASSEM assembles the matrix and right-hand side using piecewise quadratics. c c Discussion: c c The matrix is known to be banded. A special matrix storage format c is used to reduce the space required. Details of this format are c discussed in the routine DGB_FA. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer LDA, the leading dimension of A. c c Output, double precision A(LDA,NUNK), the NUNK by NUNK coefficient c matrix, stored in a compressed format. c c Output, double precision F(NUNK), the right hand side. c c Input, integer IB, the half-bandwidth of the matrix. c c Input, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c c Input, double precision XQ(NEL,NQ), YQ(NEL,NQ), the X and Y c coordinates of the quadrature points in each element. c c Input, double precision WQ(NQ), quadrature weights. c c Input, double precision AREA(NEL), the area of each element. c c Input, integer NODE(NEL,NNODES); NODE(I,J) is the global node index of c the local node J in element I. c c Input, integer INDX(NP), gives the index of the unknown quantity c associated with the given node. c c Input, integer NNODES, the number of nodes used to form one element. c c Input, integer NUNK, the number of unknowns. c c Input, integer NQ, the number of quadrature points used in assembly. c c Input, integer NEL, the number of elements. c c Input, integer NP, the number of nodes. c c Local parameters: c c Local, double precision BB, BX, BY, the value of some basis function c and its first derivatives at a quadrature point. c c Local, double precision BBB, BBX, BBY, the value of another basis c function and its first derivatives at a quadrature point. c c Local, integer NP, the number of global nodes. c implicit none integer ib integer lda integer nel integer nnodes integer np integer nq integer nunk double precision ar double precision area(nel) double precision a(lda,nunk) double precision bb double precision bbb double precision bbx double precision bby double precision bx double precision by double precision f(nunk) integer i integer indx(np) integer ip integer ipp integer iq integer it integer itest integer itrial integer j integer node(nel,nnodes) double precision rhs double precision wq(nq) double precision x double precision xc(np) double precision xq(nel,nq) double precision y double precision yc(np) double precision yq(nel,nq) c c Initialize the arrays to zero. c do i = 1, nunk f(i) = 0.0D+00 end do do j = 1, nunk do i = 1, 3*ib + 1 a(i,j) = 0.0D+00 end do end do c c The actual values of A and F are determined by summing up c contributions from all the elements. c do it = 1, nel do iq = 1, nq x = xq(it,iq) y = yq(it,iq) ar = area(it) * wq(iq) do itest = 1, nnodes ip = node(it,itest) i = indx(ip) if ( 0 .lt. i ) then call qbf ( x, y, it, itest, bb, bx, by, xc, yc, node, & nel, nnodes, np ) f(i) = f(i) + ar * rhs ( x, y ) * bb c c We are about to compute a contribution associated with the c I-th test function and the J-th basis function, and add this c to the entry A(I,J). c c Because of the compressed storage of the matrix, the element c will actually be stored in A(I-J+2*IB+1,J). c do itrial = 1, nnodes ipp = node(it,itrial) j = indx(ipp) if ( 0 .lt. j ) then call qbf ( x, y, it, itrial, bbb, bbx, bby, & xc, yc, node, nel, nnodes, np ) a(i-j+2*ib+1,j) = a(i-j+2*ib+1,j) & + ar * ( bx * bbx + by * bby ) end if end do end if end do end do end do return end subroutine dgb_fa ( lda, n, ml, mu, a, pivot, info ) c*********************************************************************72 c cc DGB_FA performs a LINPACK-style PLU factorization of an DGB matrix. c c Discussion: c c The DGB storage format is for an M by N banded matrix, with lower c bandwidth ML and upper bandwidth MU. Storage includes room for ML c extra superdiagonals, which may be required to store nonzero entries c generated during Gaussian elimination. c c The original M by N matrix is "collapsed" downward, so that diagonals c become rows of the storage array, while columns are preserved. The c collapsed array is logically 2*ML+MU+1 by N. c c The following program segment will set up the input. c c m = ml + mu + 1 c do j = 1, n c i1 = max ( 1, j-mu ) c i2 = min ( n, j+ml ) c do i = i1, i2 c k = i - j + m c a(k,j) = afull(i,j) c end do c end do c c This uses rows ML+1 through 2*ML+MU+1 of the array A. c In addition, the first ML rows in the array are used for c elements generated during the triangularization. c c The ML+MU by ML+MU upper left triangle and the c ML by ML lower right triangle are not referenced. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c Fortran77 original version by Dongarra, Bunch, Moler, Stewart. c This version by John Burkardt c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979 c c Parameters: c c Input, integer LDA, the leading dimension of A. c c Input, integer N, the order of the matrix. c N must be positive. c c Input, integer ML, MU, the lower and upper bandwidths. c ML and MU must be nonnegative, and no greater than N-1. c c Input/output, double precision A(2*ML+MU+1,N), on input, the matrix c in band storage, on output, information about the LU factorization. c c Output, integer PIVOT(N), the pivot vector. c c Output, integer INFO, singularity flag. c 0, no singularity detected. c nonzero, the factorization failed on the INFO-th step. c implicit none integer lda integer ml integer mu integer n double precision a(lda,n) integer i integer i0 integer info integer pivot(n) integer j integer j0 integer j1 integer ju integer jz integer k integer l integer lm integer m integer mm m = ml + mu + 1 info = 0 c c Zero out the initial fill-in columns. c j0 = mu + 2 j1 = min ( n, m ) - 1 do jz = j0, j1 i0 = m + 1 - jz do i = i0, ml a(i,jz) = 0.0D+00 end do end do jz = j1 ju = 0 do k = 1, n-1 c c Zero out the next fill-in column. c jz = jz + 1 if ( jz .le. n ) then do i = 1, ml a(i,jz) = 0.0D+00 end do end if c c Find L = pivot index. c lm = min ( ml, n-k ) l = m do j = m+1, m+lm if ( abs ( a(l,k) ) .lt. abs ( a(j,k) ) ) then l = j end if end do pivot(k) = l + k - m c c Zero pivot implies this column already triangularized. c if ( a(l,k) .eq. 0.0D+00 ) then info = k write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DGB_FA - Fatal error!' write ( *, '(a,i6)' ) ' Zero pivot on step ', info return end if c c Interchange if necessary. c call r8_swap ( a(l,k), a(m,k) ) c c Compute multipliers. c do i = m+1, m+lm a(i,k) = - a(i,k) / a(m,k) end do c c Row elimination with column indexing. c ju = max ( ju, mu+pivot(k) ) ju = min ( ju, n ) mm = m do j = k+1, ju l = l - 1 mm = mm - 1 if ( l .ne. mm ) then call r8_swap ( a(l,j), a(mm,j) ) end if do i = 1, lm a(mm+i,j) = a(mm+i,j) + a(mm,j) * a(m+i,k) end do end do end do pivot(n) = n if ( a(m,n) .eq. 0.0D+00 ) then info = n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DGB_FA - Fatal error!' write ( *, '(a,i6)' ) ' Zero pivot on step ', info end if return end subroutine dgb_print_some ( lda, m, n, ml, mu, a, ilo, jlo, & ihi, jhi, title ) c*********************************************************************72 c cc DGB_PRINT_SOME prints some of a DGB matrix. c c Discussion: c c The DGB storage format is for an M by N banded matrix, with lower c bandwidth ML and upper bandwidth MU. Storage includes room for ML c extra superdiagonals, which may be required to store nonzero entries c generated during Gaussian elimination. c c The original M by N matrix is "collapsed" downward, so that diagonals c become rows of the storage array, while columns are preserved. The c collapsed array is logically 2*ML+MU+1 by N. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer LDA, the leading dimension of A. c c Input, integer M, the number of rows of the matrix. c M must be positive. c c Input, integer N, the number of columns of the matrix. c N must be positive. c c Input, integer ML, MU, the lower and upper bandwidths. c ML and MU must be nonnegative, and no greater than min(M,N)-1.. c c Input, double precision A(2*ML+MU+1,N), the SGB matrix. c c Input, integer ILO, JLO, IHI, JHI, the first row and c column, and the last row and column to be printed. c c Input, character * ( * ) TITLE, the title. c implicit none integer incx parameter ( incx = 5 ) integer lda integer ml integer mu integer n double precision a(lda,n) character ( len = 14 ) ctemp(incx) integer i integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2 integer j2hi integer j2lo integer jhi integer jlo integer m character * ( * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) title c c Print the columns of the matrix, in strips of 5. c do j2lo = jlo, jhi, incx j2hi = j2lo + incx - 1 j2hi = min ( j2hi, n ) j2hi = min ( j2hi, jhi ) inc = j2hi + 1 - j2lo write ( *, '(a)' ) ' ' do j = j2lo, j2hi j2 = j + 1 - j2lo write ( ctemp(j2), '(i7,7x)' ) j end do write ( *, '(a,5a14)' ) ' Col: ', ( ctemp(j2), j2 = 1, inc ) write ( *, '(a)' ) ' Row' write ( *, '(a)' ) ' ---' c c Determine the range of the rows in this strip. c i2lo = max ( ilo, 1 ) i2lo = max ( i2lo, j2lo - mu ) i2hi = min ( ihi, m ) i2hi = min ( i2hi, j2hi + ml ) do i = i2lo, i2hi c c Print out (up to) 5 entries in row I, that lie in the current strip. c do j2 = 1, inc j = j2lo - 1 + j2 if ( ml .lt. i-j .or. mu .lt. j-i ) then ctemp(j2) = ' ' else write ( ctemp(j2), '(g14.6)' ) a(i-j+ml+mu+1,j) end if end do write ( *, '(i5,1x,5a14)' ) i, ( ctemp(j2), j2 = 1, inc ) end do end do return end subroutine dgb_sl ( lda, n, ml, mu, a, pivot, b, job ) c*********************************************************************72 c cc DGB_SL solves a system factored by DGB_FA. c c Discussion: c c The DGB storage format is for an M by N banded matrix, with lower c bandwidth ML and upper bandwidth MU. Storage includes room for ML c extra superdiagonals, which may be required to store nonzero entries c generated during Gaussian elimination. c c The original M by N matrix is "collapsed" downward, so that diagonals c become rows of the storage array, while columns are preserved. The c collapsed array is logically 2*ML+MU+1 by N. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c Fortran77 original version by Dongarra, Bunch, Moler, Stewart. c This version by John Burkardt c c Reference: c c Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, c LINPACK User's Guide, c SIAM, 1979 c c Parameters: c c Input, integer LDA, the leading dimension of A. c c Input, integer N, the order of the matrix. c N must be positive. c c Input, integer ML, MU, the lower and upper bandwidths. c ML and MU must be nonnegative, and no greater than N-1. c c Input, double precision A(2*ML+MU+1,N), the LU factors from DGB_FA. c c Input, integer PIVOT(N), the pivot vector from DGB_FA. c c Input/output, double precisionl B(N). c On input, the right hand side vector. c On output, the solution. c c Input, integer JOB. c 0, solve A * x = b. c nonzero, solve A' * x = b. c implicit none integer lda integer ml integer mu integer n double precision a(lda,n) double precision b(n) integer i integer job integer k integer l integer la integer lb integer lm integer m integer pivot(n) m = mu + ml + 1 c c Solve A * x = b. c if ( job .eq. 0 ) then c c Solve L * Y = B. c if ( 1 .le. ml ) then do k = 1, n-1 lm = min ( ml, n-k ) l = pivot(k) if ( l .ne. k ) then call r8_swap ( b(l), b(k) ) end if do i = 1, lm b(k+i) = b(k+i) + b(k) * a(m+i,k) end do end do end if c c Solve U * X = Y. c do k = n, 1, -1 b(k) = b(k) / a(m,k) lm = min ( k, m ) - 1 la = m - lm lb = k - lm do i = 0, lm-1 b(lb+i) = b(lb+i) - b(k) * a(la+i,k) end do end do c c Solve A' * X = B. c else c c Solve U' * Y = B. c do k = 1, n lm = min ( k, m ) - 1 la = m - lm lb = k - lm do i = 0, lm-1 b(k) = b(k) - a(la+i,k) * b(lb+i) end do b(k) = b(k) / a(m,k) end do c c Solve L' * X = Y. c if ( 1 .le. ml ) then do k = n-1, 1, -1 lm = min ( ml, n-k ) do i = 1, lm b(k) = b(k) + a(m+i,k) * b(k+i) end do l = pivot(k) if ( l .ne. k ) then call r8_swap ( b(l), b(k) ) end if end do end if end if return end subroutine eror ( area, node, indx, xc, yc, f, nel, & nnodes, nqe, nunk, np, el2, eh1 ) c*********************************************************************72 c cc EROR calculates the L2 and H1-seminorm errors. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, double precision AREA(NEL), the area of each element. c c Input, integer NODE(NEL,NNODES); NODE(I,J) is the global node index of c the local node J in element I. c c Input, integer INDX(NP), gives the index of the unknown quantity c associated with the given node. c c Input, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c c Input, double precision F(NUNK), the coefficients of the solution. c c Input, integer NEL, the number of elements. c c Input, integer NNODES, the number of nodes used to form one element. c c Input, integer NQE, the number of points in the quadrature rule. c This is actually fixed at 13. c c Input, integer NUNK, the number of unknowns. c c Input, integer NP, the number of nodes. c c Output, double precision EL2, the L2 error. c c Output, double precision EH1, the H1 seminorm error. c c Local Parameters: c c Local, double precision AR, the weight for a given quadrature point c in a given element. c c Local, double precision BB, BX, BY, a basis function and its first c derivatives evaluated at a particular quadrature point. c c Local, double precision EH1, the H1 seminorm error. c c Local, double precision EL2, the L2 error. c c Local, double precision UEX, UEXX, UEXY, the exact solution and its first c derivatives evaluated at a particular quadrature point. c c Local, double precision UH, UHX, UHY, the computed solution and its first c derivatives evaluated at a particular quadrature point. c c Local, double precision WQE(NQE), stores the quadrature weights. c c Local, double precision X, Y, the coordinates of a particular c quadrature point. c c Local, double precision XQE(NQE), YQE(NQE), stores the location c of quadrature points in a given element. c implicit none integer nel integer nnodes integer np integer nqe integer nunk double precision ar double precision area(nel) double precision bb double precision bx double precision by double precision eh1 double precision el2 double precision f(nunk) integer i integer in1 integer indx(np) integer ip integer iq integer it integer node(nel,nnodes) double precision uex double precision uexx double precision uexy double precision uh double precision uhx double precision uhy double precision wqe(nqe) double precision x double precision x1 double precision xc(np) double precision xqe(nqe) double precision y double precision y1 double precision yc(np) double precision yqe(nqe) el2 = 0.0D+00 eh1 = 0.0D+00 do it = 1, nel call quad13 ( xc, yc, node, it, wqe, xqe, yqe, nel, nqe, & nnodes, np ) do iq = 1, nqe ar = area(it) * wqe(iq) x = xqe(iq) y = yqe(iq) uh = 0.0D+00 uhx = 0.0D+00 uhy = 0.0D+00 do in1 = 1, nnodes ip = node(it,in1) call qbf ( x, y, it, in1, bb, bx, by, xc, yc, node, nel, & nnodes, np ) x1 = xc(ip) y1 = yc(ip) i = indx(ip) if ( 0 .lt. i ) then uh = uh + bb * f(i) uhx = uhx + bx * f(i) uhy = uhy + by * f(i) end if end do call exact ( x, y, uex, uexx, uexy ) el2 = el2 + ( uh - uex )**2 * ar eh1 = eh1 + ( ( uhx - uexx )**2 + ( uhy - uexy )**2 ) * ar end do end do el2 = sqrt ( el2 ) eh1 = sqrt ( eh1 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) '*********************************************' write ( *, '(a)' ) '* *' write ( *, '(a)' ) '* EROR: *' write ( *, '(a,g14.6,a)' ) & '* L2 error = ', el2, ' *' write ( *, '(a,g14.6,a)' ) & '* H1-seminorm error = ', eh1, ' *' write ( *, '(a)' ) '* *' write ( *, '(a)' ) '*********************************************' return end subroutine exact ( x, y, uex, uexx, uexy ) c*********************************************************************72 c cc EXACT calculates the exact solution and its first derivatives. c c Discussion: c c The function specified here depends on the problem being c solved. This is one of the routines that a user will c normally want to change. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, double precision X, Y, the coordinates of a point c in the region, at which the right hand side of the c differential equation is to be evaluated. c c Output, double precision UEX, UEXX, UEXY, the value of c the exact solution U and its derivatives dUdX c and dUdY at the point (X,Y). c implicit none double precision pi parameter ( pi = 3.1415926535897932384D+00 ) double precision uex double precision uexx double precision uexy double precision x double precision y uex = sin ( pi * x ) * sin ( pi * y ) uexx = pi * cos ( pi * x ) * sin ( pi * y ) uexy = pi * sin ( pi * x ) * cos ( pi * y ) return end subroutine geom ( xc, yc, xq, yq, wq, area, node, indx, & nel, nunk, np, ib, nnodes, xl, xr, yb, yt, nx, ny ) c*********************************************************************72 c cc GEOM sets up geometric information for a rectangular domain. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Output, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c c Output, double precision XQ(NEL,3), YQ(NEL,3), the X and Y c coordinates of the quadrature points in each element. c c Output, double precision WQ(3), quadrature weights. c c Output, double precision AREA(NEL), the area of each element. c c Output, integer NODE(NEL,NNODES); NODE(I,J) is the global node index of c the local node J in element I. c c Output, integer INDX(NP), gives the index of the unknown quantity c associated with the given node. c c Input, integer NEL, the number of elements. c c Output, integer NUNK, the number of unknowns. c c Input, integer NP, the number of nodes. c c Output, integer IB, the half-bandwidth of the matrix. c c Input, integer NNODES, the number of nodes used to form one element. c c Input, double precision XL, XR, YB, YT, the X coordinates of c the left and right sides of the rectangle, and the Y coordinates c of the bottom and top of the rectangle. c c Input, integer NX, NY, the number of points in the X and Y directions. c c Local parameters: c c Local, integer IP, counts the nodes. c c Local, integer IT1, counts the elements. c c Local, integer IU, counts the unknowns. c implicit none integer nel integer nnodes integer np double precision area(nel) double precision hx double precision hy integer i integer ib integer ic integer ij integer in1 integer indx(np) integer inn integer ip integer ip1 integer ip2 integer ip3 integer ipp integer it integer it1 integer it2 integer iu integer ix integer iy integer j integer jc integer node(nel,nnodes) integer nunk integer nx integer nx2 integer nxm integer ny integer ny2 integer nym double precision wq(3) double precision x1 double precision x2 double precision x3 double precision xc(np) double precision xl double precision xq(nel,3) double precision xr double precision xx double precision y1 double precision y2 double precision y3 double precision yb double precision yc(np) double precision yq(nel,3) double precision yt double precision yy c c Set up the grid information. c nxm = nx - 1 hx = ( xr - xl ) / dble ( nxm ) nym = ny-1 hy = ( yt - yb ) / dble ( nym ) yy = -hy / 2.0D+00 nx2 = nx + nxm ny2 = ny + nym iu = 0 ip = 0 it1 = -1 do jc = 1, ny2 yy = yy + hy / 2.0D+00 xx = -hx / 2.0D+00 iy = mod ( jc, 2 ) do ic = 1, nx2 xx = xx + hx / 2.0D+00 ix = mod ( ic, 2 ) ip = ip + 1 xc(ip) = xx yc(ip) = yy c c If JC is even, then we are at an even-numbered node in the Y direction. c These are midside nodes, and we don't need to generate a new set of elements. c c If IC is even, then we are at an even-numbered node in the X direction. c These are midside nodes, and we don't need to generate a new set of elements. c c If IC and JC are both odd, and we are not in the last row or column, then c we need to generate a new pair of elements. c if ( ix .eq. 1 .and. iy .eq. 1 ) then if ( ic .lt. nx2 .and. jc < ny2 ) then it1 = it1 + 2 node(it1,1) = ip node(it1,2) = ip + nx2 + nx2 node(it1,3) = ip + nx2 + nx2 + 2 node(it1,4) = ip + nx2 node(it1,5) = ip + nx2 + nx2 + 1 node(it1,6) = ip + nx2 + 1 it2 = it1 + 1 node(it2,1) = ip node(it2,2) = ip + 2 node(it2,3) = ip + nx2 + nx2 + 2 node(it2,4) = ip + 1 node(it2,5) = ip + nx2 + 2 node(it2,6) = ip + nx2 + 1 end if end if c c Assign an index for the unknown associated with this node. c We assume homogeneous Dirichlet boundary data, so nodes on c the boundary do NOT have an associated unknown. c indx(ip) = 0 if ( 1 .lt. ic .and. ic < nx2 ) then if ( 1 .lt. jc .and. jc < ny2 ) then iu = iu + 1 indx(ip) = iu end if end if end do end do nunk = iu write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GEOM:' write ( *, '(a,i6)' ) ' Number of unknowns = ', nunk write ( *, '(a,i6)' ) ' Number of nodes = ', np write ( *, '(a,i6)' ) ' Number of elements = ', nel c c Set up quadrature information needed for assembly. c wq(1) = 1.0D+00 / 3.0D+00 wq(2) = wq(1) wq(3) = wq(1) do it = 1, nel ip1 = node(it,1) ip2 = node(it,2) ip3 = node(it,3) x1 = xc(ip1) x2 = xc(ip2) x3 = xc(ip3) y1 = yc(ip1) y2 = yc(ip2) y3 = yc(ip3) xq(it,1) = 0.5D+00 * ( x1 + x2 ) xq(it,2) = 0.5D+00 * ( x2 + x3 ) xq(it,3) = 0.5D+00 * ( x1 + x3 ) yq(it,1) = 0.5D+00 * ( y1 + y2 ) yq(it,2) = 0.5D+00 * ( y2 + y3 ) yq(it,3) = 0.5D+00 * ( y1 + y3 ) area(it) = 0.5D+00 * abs ( & y1 * ( x2 - x3 ) & + y2 * ( x3 - x1 ) & + y3 * ( x1 - x2 ) ) end do c c Determine the bandwidth of the coefficient matrix. c ib = 0 do it = 1, nel do in1 = 1, nnodes ip = node(it,in1) i = indx(ip) if ( 0 .lt. i ) then do inn = 1, nnodes ipp = node(it,inn) j = indx(ipp) ij = j - i if ( ib .lt. ij ) then ib = ij end if end do end if end do end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GEOM:' write ( *, '(a,i8)' ) ' Total bandwidth is ', 3 * ib + 1 return end subroutine get_unit ( iunit ) c*********************************************************************72 c cc GET_UNIT returns a free Fortran unit number. c c Discussion: c c A "free" Fortran unit number is a value between 1 and 99 which c is not currently associated with an I/O device. A free Fortran unit c number is needed in order to open a file with the OPEN command. c c If IUNIT = 0, then no free Fortran unit could be found, although c all 99 units were checked (except for units 5, 6 and 9, which c are commonly reserved for console I/O). c c Otherwise, IUNIT is a value between 1 and 99, representing a c free Fortran unit. Note that GET_UNIT assumes that units 5 and 6 c are special, and will never return those values. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 02 September 2013 c c Author: c c John Burkardt c c Parameters: c c Output, integer IUNIT, the free unit number. c implicit none integer i integer iunit logical value iunit = 0 do i = 1, 99 if ( i .ne. 5 .and. i .ne. 6 .and. i .ne. 9 ) then inquire ( unit = i, opened = value, err = 10 ) if ( .not. value ) then iunit = i return end if end if 10 continue end do return end function i4_huge ( ) c*********************************************************************72 c cc I4_HUGE returns a "huge" I4. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Output, integer I4_HUGE, a huge number. c implicit none integer i4_huge i4_huge = 2147483647 return end subroutine i4vec_print_some ( n, a, max_print, title ) c*********************************************************************72 c cc I4VEC_PRINT_SOME prints "some" of an integer vector. c c Discussion: c c The user specifies MAX_PRINT, the maximum number of lines to print. c c If N, the size of the vector, is no more than MAX_PRINT, then c the entire vector is printed, one entry per line. c c Otherwise, if possible, the first MAX_PRINT-2 entries are printed, c followed by a line of periods suggesting an omission, c and the last entry. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries of the vector. c c Input, integer A(N), the vector to be printed. c c Input, integer MAX_PRINT, the maximum number of lines to print. c c Input, character ( len = * ) TITLE, an optional title. c implicit none integer n integer a(n) integer i integer max_print integer s_len_trim character ( len = * ) title if ( max_print .le. 0 ) then return end if if ( n .le. 0 ) then return end if if ( 0 .lt. s_len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) title write ( *, '(a)' ) ' ' end if if ( n .le. max_print ) then do i = 1, n write ( *, '(i6,2x,i10)' ) i, a(i) end do else if ( 3 .le. max_print ) then do i = 1, max_print-2 write ( *, '(i6,2x,i10)' ) i, a(i) end do write ( *, '(a)' ) '...... ..............' i = n write ( *, '(i6,2x,i10)' ) i, a(i) else do i = 1, max_print - 1 write ( *, '(i6,2x,i10)' ) i, a(i) end do i = max_print write ( *, '(i6,2x,i10,2x,a)' ) i, a(i), '...more entries...' end if return end subroutine mesh_eps ( file_name, np, xc, yc, nnodes, & nel, node ) c*********************************************************************72 c cc MESH_EPS creates an EPS file containing an image of the mesh. c c Discussion: c c This routine has been specialized to deal correctly ONLY with c a mesh of 6 node elements, with the property that starting c from local node 1 and traversing the edges of the element will c result in encountering local nodes 1, 4, 2, 5, 3, 6 in that order. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character ( len = * ) FILE_NAME, the name of the file to create. c c Input, integer NP, the number of nodes. c c Input, double precision XC(NP), YC(NP), the XY coordinates of the nodes. c c Input, integer NNODES, the number of local nodes per element. c c Input, integer NEL, the number of elements. c c Input, integer NODE(NEL,NNODES), the element->node data. c implicit none integer nel integer nnodes integer np integer eps_unit character ( len = * ) file_name integer i integer ip1 integer j integer k integer node(nel,nnodes) integer order(6) character ( len = 20 ) string double precision xc(np) double precision x_ave integer x_eps double precision yc(np) double precision y_ave integer y_eps order(1) = 1 order(2) = 4 order(3) = 2 order(4) = 5 order(5) = 3 order(6) = 6 call get_unit ( eps_unit ) open ( unit = eps_unit, file = file_name, & status = 'unknown' ) write ( eps_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( eps_unit, '(a)' ) '%%Creator: mesh_eps' write ( eps_unit, '(a,a)' ) '%%Title: ', file_name write ( eps_unit, '(a)' ) '%%Pages: 1' write ( eps_unit, '(a)' ) '%%BoundingBox: 36 36 576 756' write ( eps_unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( eps_unit, '(a)' ) '%%LanguageLevel: 1' write ( eps_unit, '(a)' ) '%%EndComments' write ( eps_unit, '(a)' ) '%%BeginProlog' write ( eps_unit, '(a)' ) '/inch {72 mul} def' write ( eps_unit, '(a)' ) '%%EndProlog' write ( eps_unit, '(a)' ) '%%Page: 1 1' write ( eps_unit, '(a)' ) 'save' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Set RGB line color.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.9000 0.9000 0.9000 setrgbcolor' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw a gray border around the page.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) 'newpath' write ( eps_unit, '(a)' ) ' 36 126 moveto' write ( eps_unit, '(a)' ) ' 576 126 lineto' write ( eps_unit, '(a)' ) ' 576 666 lineto' write ( eps_unit, '(a)' ) ' 36 666 lineto' write ( eps_unit, '(a)' ) ' 36 126 lineto' write ( eps_unit, '(a)' ) 'stroke' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Set RGB line color.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Define a clipping polygon' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 36 126 moveto' write ( eps_unit, '(a)' ) ' 576 126 lineto' write ( eps_unit, '(a)' ) ' 576 666 lineto' write ( eps_unit, '(a)' ) ' 36 666 lineto' write ( eps_unit, '(a)' ) ' 36 126 lineto' write ( eps_unit, '(a)' ) 'clip newpath' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw filled dots at each node:' write ( eps_unit, '(a)' ) '%' do i = 1, np x_eps = int ( ( 1.0 - xc(i) ) * 61.0 + xc(i) * 551.0 ) y_eps = int ( ( 1.0 - yc(i) ) * 151.0 + yc(i) * 641.0 ) write ( eps_unit, '(a,i4,2x,i4,a)' ) & 'newpath ', x_eps, y_eps, ' 5 0 360 arc closepath fill' end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the nodes:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) & '/Times-Roman findfont 0.20 inch scalefont setfont' do i = 1, np x_eps = int ( ( 1.0 - xc(i) ) * 61.0 + xc(i) * 551.0 ) y_eps = int ( ( 1.0 - yc(i) ) * 151.0 + yc(i) * 641.0 ) write ( string, '(i4)' ) i write ( eps_unit, '(i4,2x,i4,a,a,a)' ) x_eps, y_eps+5, & ' moveto (', string, ') show' end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw the element sides:' write ( eps_unit, '(a)' ) '%' do j = 1, nel k = node(j,order(1)) x_eps = int ( ( 1.0 - xc(k) ) * 61.0 + xc(k) * 551.0 ) y_eps = int ( ( 1.0 - yc(k) ) * 151.0 + yc(k) * 641.0 ) write ( eps_unit, '(a,i4,2x,i4,a)' ) & 'newpath ', x_eps, y_eps, ' moveto' do i = 1, nnodes ip1 = mod ( i, nnodes ) + 1; k = node(j,order(ip1)) x_eps = int ( ( 1.0 - xc(k) ) * 61.0 + xc(k) * 551.0 ) y_eps = int ( ( 1.0 - yc(k) ) * 151.0 + yc(k) * 641.0 ) write ( eps_unit, '(i4,2x,i4,a)' ) x_eps, y_eps, ' lineto' end do write ( eps_unit, '(a)' ) 'stroke' end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the elements:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) & '/Times-Roman findfont 0.30 inch scalefont setfont' do j = 1, nel x_ave = 0.0 y_ave = 0.0 do i = 1, nnodes k = node(j,i) x_ave = x_ave + xc(k) y_ave = y_ave + yc(k) end do x_ave = x_ave / real ( nnodes ) y_ave = y_ave / real ( nnodes ) x_eps = int ( ( 1.0 - x_ave ) * 61.0 + x_ave * 551.0 ) y_eps = int ( ( 1.0 - y_ave ) * 151.0 + y_ave * 641.0 ) write ( string, '(i4)' ) j write ( eps_unit, '(i4,2x,i4,a,a,a)' ) x_eps, y_eps, & ' moveto (', string, ') show' end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) 'restore showpage' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% End of page' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '%%Trailer' write ( eps_unit, '(a)' ) '%%EOF' close ( unit = eps_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MESH_EPS:' write ( *, '(a)' ) ' An encapsulated PostScript file was created' write ( *, '(a)' ) ' with an image of the nodes and elements.' write ( *, '(a)' ) ' The file is named "', file_name, '".' return end subroutine qbf ( x, y, it, inode, bb, bx, by, xc, yc, & node, nel, nnodes, np ) c*********************************************************************72 c cc QBF evaluates the quadratic basis functions. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, double precision X, Y, the (global) coordinates of the point c at which the basis function is to be evaluated. c c Input, integer IT, the index of the element which contains the point. c c Input, integer INODE, the local index (between 1 and 6) that c specifies which basis function is to be evaluated. c c Output, double precision BB, BX, BY, the value of the basis function c and its X and Y derivatives at (X,Y). c c Input, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c c Input, integer NODE(NEL,NNODES); NODE(I,J) is the global node index of c the local node J in element I. c c Input, integer NEL, the number of elements. c c Input, integer NNODES, the number of nodes used to form one element. c c Input, integer NP, the number of nodes. c implicit none integer nel integer nnodes integer np double precision bb double precision bx double precision by double precision c double precision d integer i1 integer i2 integer i3 integer in1 integer in2 integer in3 integer inn integer inode integer it integer j1 integer j2 integer j3 integer node(nel,nnodes) double precision s double precision t double precision x double precision xc(np) double precision y double precision yc(np) c c Basis functions associated with "corner" nodes. c if ( inode .le. 3 ) then in1 = inode in2 = mod ( inode, 3 ) + 1 in3 = mod ( inode+1, 3 ) + 1 i1 = node(it,in1) i2 = node(it,in2) i3 = node(it,in3) d = ( xc(i2) - xc(i1) ) * ( yc(i3) - yc(i1) ) & - ( xc(i3) - xc(i1) ) * ( yc(i2) - yc(i1) ) t = 1.0D+00 & + ( ( yc(i2) - yc(i3) ) * ( x - xc(i1) ) & + ( xc(i3) - xc(i2) ) * ( y - yc(i1) ) ) / d bb = t * ( 2.0D+00 * t - 1.0D+00 ) bx = ( yc(i2) - yc(i3) ) * ( 4.0D+00 * t - 1.0D+00 ) / d by = ( xc(i3) - xc(i2) ) * ( 4.0D+00 * t - 1.0D+00 ) / d c c Basis functions associated with midside nodes. c else if ( inode .le. 6 ) then inn = inode - 3 in1 = inn in2 = mod ( inn, 3 ) + 1 in3 = mod ( inn+1,3 ) + 1 i1 = node(it,in1) i2 = node(it,in2) i3 = node(it,in3) j1 = i2 j2 = i3 j3 = i1 d = ( xc(i2) - xc(i1) ) * ( yc(i3) - yc(i1) ) & - ( xc(i3) - xc(i1) ) * ( yc(i2) - yc(i1) ) c = ( xc(j2) - xc(j1) ) * ( yc(j3) - yc(j1) ) & - ( xc(j3) - xc(j1) ) * ( yc(j2) - yc(j1) ) t = 1.0D+00 & + ( ( yc(i2) - yc(i3) ) * ( x - xc(i1) ) & + ( xc(i3) - xc(i2) ) * ( y - yc(i1) ) ) / d s = 1.0D+00 & + ( ( yc(j2) - yc(j3) ) * ( x - xc(j1) ) & + ( xc(j3) - xc(j2) ) * ( y - yc(j1) ) ) / c bb = 4.0D+00 * s * t bx = 4.0D+00 * ( t * ( yc(j2) - yc(j3) ) / c & + s * ( yc(i2) - yc(i3) ) / d ) by = 4.0D+00 * ( t * ( xc(j3) - xc(j2) ) / c & + s * ( xc(i3) - xc(i2) ) / d ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'QBF - Fatal error!' write ( *, '(a,i6)' ) & ' Request for local basis function INODE = ', inode stop end if return end subroutine quad13 ( xc, yc, node, it, wqe, xqe, yqe, nel, & nqe, nnodes, np ) c*********************************************************************72 c cc QUAD13 sets up quadrature information for a 13-point rule in a given element. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c c Input, integer NODE(NEL,NNODES); NODE(I,J) is the global node index of c the local node J in element I. c c Input, integer IT, the index of the element for which the quadrature c points are to be computed. c c Output, double precision WQE(13), the quadrature weights. c c Output, double precision XQE(NQE), YQE(NQE), the X and Y coordinates c of the quadrature points. c c Input, integer NEL, the number of elements. c c Input, integer NQE, the number of points in the quadrature rule. c This is actually fixed at 13. c c Input, integer NNODES, the number of nodes used to form one element. c c Input, integer NP, the number of nodes. c implicit none integer nel integer nnodes integer np integer nqe integer i integer ii integer iii integer ip1 integer ip2 integer ip3 integer it integer node(nel,nnodes) double precision wqe(nqe) double precision xc(np) double precision xqe(nqe) double precision yc(np) double precision yqe(nqe) double precision x1 double precision x2 double precision x3 double precision y1 double precision y2 double precision y3 double precision z1 double precision z2 double precision z3 double precision z4 double precision z5 double precision z6 double precision z7 do i = 1, 3 wqe(i) = 0.175615257433204D+00 ii = i + 3 wqe(ii) = 0.053347235608839D+00 ii = i + 6 iii = ii + 3 wqe(ii) = 0.077113760890257D+00 wqe(iii) = wqe(ii) end do wqe(13) = -0.14957004446767D+00 z1 = 0.479308067841923D+00 z2 = 0.260345966079038D+00 z3 = 0.869739794195568D+00 z4 = 0.065130102902216D+00 z5 = 0.638444188569809D+00 z6 = 0.312865496004875D+00 z7 = 0.048690315425316D+00 ip1 = node(it,1) ip2 = node(it,2) ip3 = node(it,3) x1 = xc(ip1) x2 = xc(ip2) x3 = xc(ip3) y1 = yc(ip1) y2 = yc(ip2) y3 = yc(ip3) xqe( 1) = z1 * x1 + z2 * x2 + z2 * x3 yqe( 1) = z1 * y1 + z2 * y2 + z2 * y3 xqe( 2) = z2 * x1 + z1 * x2 + z2 * x3 yqe( 2) = z2 * y1 + z1 * y2 + z2 * y3 xqe( 3) = z2 * x1 + z2 * x2 + z1 * x3 yqe( 3) = z2 * y1 + z2 * y2 + z1 * y3 xqe( 4) = z3 * x1 + z4 * x2 + z4 * x3 yqe( 4) = z3 * y1 + z4 * y2 + z4 * y3 xqe( 5) = z4 * x1 + z3 * x2 + z4 * x3 yqe( 5) = z4 * y1 + z3 * y2 + z4 * y3 xqe( 6) = z4 * x1 + z4 * x2 + z3 * x3 yqe( 6) = z4 * y1 + z4 * y2 + z3 * y3 xqe( 7) = z5 * x1 + z6 * x2 + z7 * x3 yqe( 7) = z5 * y1 + z6 * y2 + z7 * y3 xqe( 8) = z5 * x1 + z7 * x2 + z6 * x3 yqe( 8) = z5 * y1 + z7 * y2 + z6 * y3 xqe( 9) = z6 * x1 + z5 * x2 + z7 * x3 yqe( 9) = z6 * y1 + z5 * y2 + z7 * y3 xqe(10) = z6 * x1 + z7 * x2 + z5 * x3 yqe(10) = z6 * y1 + z7 * y2 + z5 * y3 xqe(11) = z7 * x1 + z5 * x2 + z6 * x3 yqe(11) = z7 * y1 + z5 * y2 + z6 * y3 xqe(12) = z7 * x1 + z6 * x2 + z5 * x3 yqe(12) = z7 * y1 + z6 * y2 + z5 * y3 xqe(13) = ( x1 + x2 + x3 ) / 3.0D+00 yqe(13) = ( y1 + y2 + y3 ) / 3.0D+00 return end function r8_is_int ( r ) c*********************************************************************72 c cc R8_IS_INT determines if an R8 represents an integer value. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Author: c c John Burkardt c c Parameters: c c Input, double precision R, the number to be checked. c c Output, logical R8_IS_INT, is TRUE if R is an integer value. c implicit none integer i4_huge double precision r logical r8_is_int if ( dble ( i4_huge ( ) ) .lt. r ) then r8_is_int = .false. else if ( r .lt. - dble ( i4_huge ( ) ) ) then r8_is_int = .false. else if ( r .eq. dble ( int ( r ) ) ) then r8_is_int = .true. else r8_is_int = .false. end if return end subroutine r8_swap ( x, y ) c*********************************************************************72 c cc R8_SWAP switches two R8's. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input/output, double precision X, Y. On output, the values of X and c Y have been interchanged. c implicit none double precision x double precision y double precision z z = x x = y y = z return end subroutine r8vec_print_some ( n, a, max_print, title ) c*********************************************************************72 c cc R8VEC_PRINT_SOME prints "some" of a double precision vector. c c Discussion: c c The user specifies MAX_PRINT, the maximum number of lines to print. c c If N, the size of the vector, is no more than MAX_PRINT, then c the entire vector is printed, one entry per line. c c Otherwise, if possible, the first MAX_PRINT-2 entries are printed, c followed by a line of periods suggesting an omission, c and the last entry. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of entries of the vector. c c Input, double precision A(N), the vector to be printed. c c Input, integer MAX_PRINT, the maximum number of lines to print. c c Input, character ( len = * ) TITLE, an optional title. c implicit none integer n double precision a(n) integer i integer max_print integer s_len_trim character ( len = * ) title if ( max_print .le. 0 ) then return end if if ( n .le. 0 ) then return end if if ( 0 .lt. s_len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) title write ( *, '(a)' ) ' ' end if if ( n .le. max_print ) then do i = 1, n write ( *, '(i6,2x,g14.6)' ) i, a(i) end do else if ( 3 .le. max_print ) then do i = 1, max_print-2 write ( *, '(i6,2x,g14.6)' ) i, a(i) end do write ( *, '(a)' ) '...... ..............' i = n write ( *, '(i6,2x,g14.6)' ) i, a(i) else do i = 1, max_print-1 write ( *, '(i6,2x,g14.6)' ) i, a(i) end do i = max_print write ( *, '(i6,2x,g14.6,a)' ) i, a(i), '...more entries...' end if return end function rhs ( x, y ) c*********************************************************************72 c cc RHS gives the right-hand side of the differential equation. c c Discussion: c c The function specified here depends on the problem being c solved. This is one of the routines that a user will c normally want to change. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, double precision X, Y, the coordinates of a point c in the region, at which the right hand side of the c differential equation is to be evaluated. c c Output, double precision RHS, the value of the right c hand side of the differential equation at (X,Y). c implicit none double precision pi parameter ( pi = 3.1415926535897932384 ) double precision rhs double precision x double precision y rhs = 2.0D+00 * pi * pi * sin ( pi * x ) * sin ( pi * y ) return end function s_len_trim ( s ) c*********************************************************************72 c cc S_LEN_TRIM returns the length of a string to the last nonblank. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, character ( len = * ) S, a string. c c Output, integer S_LEN_TRIM, the length of the string to the last nonblank. c implicit none integer i character ( len = * ) s integer s_len_trim do i = len ( s ), 1, -1 if ( s(i:i) .ne. ' ' ) then s_len_trim = i return end if end do s_len_trim = 0 return end subroutine solution_write ( f, indx, np, nunk, & output_filename, xc, yc ) c*********************************************************************72 c cc SOLUTION_WRITE writes the solution to a file. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c Input, double precision F(NUNK), the coefficients of the solution. c c Input, integer INDX(NP), gives the index of the unknown quantity c associated with the given node. c c Input, integer NP, the number of nodes. c c Input, integer NUNK, the number of unknowns. c c Input, character ( len = * ) OUTPUT_FILENAME, the name of the file c in which the data should be stored. c c Input, double precision XC(NP), YC(NP), the X and Y coordinates of nodes. c implicit none integer np integer nunk double precision f(nunk) integer i integer indx(np) character ( len = * ) :: output_filename integer output_unit double precision u double precision ux double precision uy double precision x double precision xc(np) double precision y double precision yc(np) call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'unknown' ) do i = 1, np if ( 0 .lt. indx(i) ) then u = f(indx(i)) else x = xc(i) y = yc(i) call exact ( x, y, u, ux, uy ) end if write ( output_unit, '(f8.4,2x,f8.4,2x,g14.6)' ) xc(i), yc(i), u end do close ( unit = output_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SOLUTION_WRITE:' write ( *, '(a)' ) ' Wrote an ASCII file' write ( *, '(a,a,a)' ) ' "', output_filename, '"' write ( *, '(a)' ) ' of the form' write ( *, '(a)' ) ' X(I), Y(I), U(I)' write ( *, '(a)' ) ' which can be contour-plotted using the' write ( *, '(a)' ) & ' MATLAB routine FEM2D_POISSON_RECTANGLE_CONTOUR.' return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Discussion: c c This Fortran77 version is made available for cases where the c Fortran90 version cannot be used. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 24 September 2008 c c Author: c c John Burkardt c c Parameters: c c None c implicit none character * ( 8 ) ampm integer d character * ( 8 ) date integer h integer m integer mm character * ( 9 ) month(12) integer n integer s character * ( 10 ) time integer y save month data month / & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' / call date_and_time ( date, time ) read ( date, '(i4,i2,i2)' ) y, m, d read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm if ( h .lt. 12 ) then ampm = 'AM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h .lt. 12 ) then ampm = 'PM' else if ( h .eq. 12 ) then if ( n .eq. 0 .and. s .eq. 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, & '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, month(m), y, h, ':', n, ':', s, '.', mm, ampm return end