subroutine fem2d_bvp_linear ( nx, ny, a, c, f, x, y, u ) c*********************************************************************72 c cc FEM2D_BVP_LINEAR solves boundary value problem on a rectangle. c c Discussion: c c The program uses the finite element method, with piecewise linear basis c functions to solve a 2D boundary value problem over a rectangle c c The following differential equation is imposed inside the region: c c - d/dx a(x,y) du/dx - d/dy a(x,y) du/dy + c(x,y) * u(x,y) = f(x,y) c c where a(x,y), c(x,y), and f(x,y) are given functions. c c On the boundary, the solution is constrained to have the value 0. c c The finite element method will use a regular grid of NX nodes in X, and c NY nodes in Y. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 20 June 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer NX, NY, the number of X and Y grid values. c c Input, function A ( X ), evaluates a(x); c c Input, function C ( X ), evaluates c(x); c c Input, function F ( X ), evaluates f(x); c c Input, double precision X(NX), Y(NY), the mesh points. c c Output, double precision U(NX,NY), the finite element coefficients, which c are also the value of the computed solution at the mesh points. c implicit none integer nx integer ny integer quad_num parameter ( quad_num = 3 ) double precision a external a double precision abscissa(quad_num) double precision amat(nx*ny,nx*ny) double precision aq double precision b(nx*ny) double precision c external c double precision cq integer e integer ex integer ey double precision f external f double precision fq integer i integer ierror integer ii integer j integer jj integer k integer mn integer n integer ne integer nw integer qx integer qy integer s integer se integer sw integer w double precision weight(quad_num) double precision wq double precision u(nx,ny) double precision vne double precision vnex double precision vney double precision vnw double precision vnwx double precision vnwy double precision vse double precision vsex double precision vsey double precision vsw double precision vswx double precision vswy double precision x(nx) double precision xe double precision xq double precision xw double precision y(ny) double precision yn double precision yq double precision ys save abscissa save weight data abscissa / & -0.774596669241483377035853079956D+00, & 0.000000000000000000000000000000D+00, & 0.774596669241483377035853079956D+00 / data weight / & 0.555555555555555555555555555556D+00, & 0.888888888888888888888888888889D+00, & 0.555555555555555555555555555556D+00 / c c Zero out the matrix and right hand side. c mn = nx * ny do j = 1, mn do i = 1, mn amat(i,j) = 0.0D+00 end do end do do i = 1, mn b(i) = 0.0D+00 end do do ex = 1, nx - 1 w = ex e = ex + 1 xw = x(w) xe = x(e) do ey = 1, ny - 1 s = ey n = ey + 1 ys = y(s) yn = y(n) sw = ( ey - 1 ) * nx + ex se = ( ey - 1 ) * nx + ex + 1 nw = ey * nx + ex ne = ey * nx + ex + 1 do qx = 1, quad_num xq = ( ( 1.0D+00 - abscissa(qx) ) * xw & + ( 1.0D+00 + abscissa(qx) ) * xe ) & / 2.0D+00 do qy = 1, quad_num yq = ( ( 1.0D+00 - abscissa(qy) ) * ys & + ( 1.0D+00 + abscissa(qy) ) * yn ) & / 2.0D+00 wq = weight(qx) * ( xe - xw ) / 2.0D+00 & * weight(qy) * ( yn - ys ) / 2.0D+00 aq = a ( xq, yq ) cq = c ( xq, yq ) fq = f ( xq, yq ) vsw = ( xe - xq ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vswx = (-1.0D+00 ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vswy = ( xe - xq ) / ( xe - xw ) & * (-1.0D+00 ) / ( yn - ys ) vse = ( xq - xw ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vsex = ( 1.0D+00 ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vsey = ( xq - xw ) / ( xe - xw ) & * (-1.0D+00 ) / ( yn - ys ) vnw = ( xe - xq ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vnwx = (-1.0D+00 ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vnwy = ( xe - xq ) / ( xe - xw ) & * ( 1.0D+00 ) / ( yn - ys ) vne = ( xq - xw ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vnex = ( 1.0D+00 ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vney = ( xq - xw ) / ( xe - xw ) & * ( 1.0D+00 ) / ( yn - ys ) amat(sw,sw) = amat(sw,sw) + wq * ( vswx * aq * vswx & + vswy * aq * vswy & + vsw * cq * vsw ) amat(sw,se) = amat(sw,se) + wq * ( vswx * aq * vsex & + vswy * aq * vsey & + vsw * cq * vse ) amat(sw,nw) = amat(sw,nw) + wq * ( vswx * aq * vnwx & + vswy * aq * vnwy & + vsw * cq * vnw ) amat(sw,ne) = amat(sw,ne) + wq * ( vswx * aq * vnex & + vswy * aq * vney & + vsw * cq * vne ) b(sw) = b(sw) + wq * ( vsw * fq ) amat(se,sw) = amat(se,sw) + wq * ( vsex * aq * vswx & + vsey * aq * vswy & + vse * cq * vsw ) amat(se,se) = amat(se,se) + wq * ( vsex * aq * vsex & + vsey * aq * vsey & + vse * cq * vse ) amat(se,nw) = amat(se,nw) + wq * ( vsex * aq * vnwx & + vsey * aq * vnwy & + vse * cq * vnw ) amat(se,ne) = amat(se,ne) + wq * ( vsex * aq * vnex & + vsey * aq * vney & + vse * cq * vne ) b(se) = b(se) + wq * ( vse * fq ) amat(nw,sw) = amat(nw,sw) + wq * ( vnwx * aq * vswx & + vnwy * aq * vswy & + vnw * cq * vsw ) amat(nw,se) = amat(nw,se) + wq * ( vnwx * aq * vsex & + vnwy * aq * vsey & + vnw * cq * vse ) amat(nw,nw) = amat(nw,nw) + wq * ( vnwx * aq * vnwx & + vnwy * aq * vnwy & + vnw * cq * vnw ) amat(nw,ne) = amat(nw,ne) + wq * ( vnwx * aq * vnex & + vnwy * aq * vney & + vnw * cq * vne ) b(nw) = b(nw) + wq * ( vnw * fq ) amat(ne,sw) = amat(ne,sw) + wq * ( vnex * aq * vswx & + vney * aq * vswy & + vne * cq * vsw ) amat(ne,se) = amat(ne,se) + wq * ( vnex * aq * vsex & + vney * aq * vsey & + vne * cq * vse ) amat(ne,nw) = amat(ne,nw) + wq * ( vnex * aq * vnwx & + vney * aq * vnwy & + vne * cq * vnw ) amat(ne,ne) = amat(ne,ne) + wq * ( vnex * aq * vnex & + vney * aq * vney & + vne * cq * vne ) b(ne) = b(ne) + wq * ( vne * fq ) end do end do end do end do c c Where a node is on the boundary, c replace the finite element equation by a boundary condition. c k = 0 do j = 1, ny do i = 1, nx k = k + 1 if ( i .eq. 1 .or. & i .eq. nx .or. & j .eq. 1 .or. & j .eq. ny ) then do jj = 1, mn amat(k,jj) = 0.0D+00 end do do ii = 1, mn amat(ii,k) = 0.0D+00 end do amat(k,k) = 1.0D+00 b(k) = 0.0D+00 end if end do end do c c Solve the linear system. c call r8mat_solve2 ( mn, amat, b, u, ierror ) return end subroutine fem2d_h1s_error_linear ( nx, ny, x, y, u, exact_ux, & exact_uy, h1s ) c*********************************************************************72 c cc FEM2D_H1S_ERROR_LINEAR: seminorm error of a finite element solution. c c Discussion: c c The finite element method has been used, over a rectangle, c involving a grid of NX*NY nodes, with piecewise linear elements used c for the basis. c c The finite element solution U(x,y) has been computed, and formulas for the c exact derivatives Vx(x,y) and Vy(x,y) are known. c c This function estimates the H1 seminorm of the error: c c H1S = sqrt ( integral ( x, y ) ( Ux(x,y) - Vx(x,y) )^2 c + ( Uy(x,y) - Vy(x,y) )^2 dx dy ) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 20 June 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer NX, NY, the number of nodes in the X c and Y directions. c c Input, double precision X(NX), Y(NY), the grid coordinates. c c Input, double precision U(NX,NY), the finite element coefficients. c c Input, function EXACT_UX(X,Y), EXACT_UY(X,Y), return the exact c derivatives with respect to X and Y, at the point (X,Y). c c Output, double precision H1S, the estimated seminorm of c the error. c implicit none integer nx integer ny integer quad_num parameter ( quad_num = 3 ) double precision abscissa(quad_num) integer e integer ex integer ey double precision exact_ux external exact_ux double precision exact_uy external exact_uy double precision exq double precision eyq double precision h1s integer n integer ne integer nw integer qx integer qy integer s integer se integer sw double precision u(nx,ny) double precision uxq double precision uyq double precision vnex double precision vney double precision vnwx double precision vnwy double precision vsex double precision vsey double precision vswx double precision vswy integer w double precision weight(quad_num) double precision wq double precision x(nx) double precision xe double precision xq double precision xw double precision y(ny) double precision yn double precision yq double precision ys save abscissa save weight data abscissa / & -0.774596669241483377035853079956D+00, & 0.000000000000000000000000000000D+00, & 0.774596669241483377035853079956D+00 / data weight / & 0.555555555555555555555555555556D+00, & 0.888888888888888888888888888889D+00, & 0.555555555555555555555555555556D+00 / h1s = 0.0D+00 do ex = 1, nx - 1 w = ex e = ex + 1 xw = x(w) xe = x(e) do ey = 1, ny - 1 s = ey n = ey + 1 ys = y(s) yn = y(n) sw = ( ey - 1 ) * nx + ex se = ( ey - 1 ) * nx + ex + 1 nw = ey * nx + ex ne = ey * nx + ex + 1 do qx = 1, quad_num xq = ( ( 1.0D+00 - abscissa(qx) ) * xw & + ( 1.0D+00 + abscissa(qx) ) * xe ) & / 2.0D+00 do qy = 1, quad_num yq = ( ( 1.0D+00 - abscissa(qy) ) * ys & + ( 1.0D+00 + abscissa(qy) ) * yn ) & / 2.0D+00 wq = weight(qx) * ( xe - xw ) / 2.0D+00 & * weight(qy) * ( yn - ys ) / 2.0D+00 vswx = (-1.0D+00 ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vswy = ( xe - xq ) / ( xe - xw ) & * (-1.0D+00 ) / ( yn - ys ) vsex = ( 1.0D+00 ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vsey = ( xq - xw ) / ( xe - xw ) & * (-1.0D+00 ) / ( yn - ys ) vnwx = (-1.0D+00 ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vnwy = ( xe - xq ) / ( xe - xw ) & * ( 1.0D+00 ) / ( yn - ys ) vnex = ( 1.0D+00 ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vney = ( xq - xw ) / ( xe - xw ) & * ( 1.0D+00 ) / ( yn - ys ) c c Note that the south-west component of U is stored in U(W,S), not U(S,W)! c uxq = u(w,s) * vswx + u(e,s) * vsex & + u(w,n) * vnwx + u(e,n) * vnex uyq = u(w,s) * vswy + u(e,s) * vsey & + u(w,n) * vnwy + u(e,n) * vney exq = exact_ux ( xq, yq ) eyq = exact_uy ( xq, yq ) h1s = h1s + wq * ( ( uxq - exq ) ** 2 & + ( uyq - eyq ) ** 2 ) end do end do end do end do h1s = sqrt ( h1s ) return end subroutine fem2d_l1_error ( nx, ny, x, y, u, exact, e1 ) c*********************************************************************72 c cc FEM2D_L1_ERROR estimates the l1 error norm of a finite element solution. c c Discussion: c c The finite element method has been used, over a rectangle, c involving a grid of NX*NY nodes, with piecewise linear elements used c for the basis. c c The finite element coefficients have been computed, and a formula for the c exact solution is known. c c This function estimates the little l1 norm of the error: c E1 = sum ( 1 <= I <= NX, 1 <= J <= NY ) c abs ( U(i,j) - EXACT(X(i),Y(j)) ) c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 20 June 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer NX, NY, the number of nodes in the X c and Y directions. c c Input, double precision X(NX), Y(NY), the grid coordinates. c c Input, double precision U(NX,NY), the finite element coefficients. c c Input, function EQ = EXACT ( X, Y ), returns the value of the exact c solution at the point (X,Y). c c Output, double precision E1, the estimated L2 norm of the error. c implicit none integer nx integer ny double precision e1 double precision exact external exact integer i integer j double precision u(nx,ny) double precision x(nx) double precision y(ny) e1 = 0.0D+0 do j = 1, ny do i = 1, nx e1 = e1 + abs ( u(i,j) - exact ( x(i), y(j) ) ) end do end do e1 = e1 / dble ( nx ) / dble ( ny ) return end subroutine fem2d_l2_error_linear ( nx, ny, x, y, u, exact, e2 ) c*********************************************************************72 c cc FEM2D_L2_ERROR_LINEAR: L2 error norm of a finite element solution. c c Discussion: c c The finite element method has been used, over a rectangle, c involving a grid of NX*NY nodes, with piecewise linear elements used c for the basis. c c The finite element coefficients have been computed, and a formula for the c exact solution is known. c c This function estimates E2, the L2 norm of the error: c c E2 = Integral ( X, Y ) ( U(X,Y) - EXACT(X,Y) )^2 dX dY c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 20 June 2014 c c Author: c c John Burkardt c c Parameters: c c Input, integer NX, NY, the number of nodes in the X c and Y directions. c c Input, double precision X(NX), Y(NY), the grid coordinates. c c Input, double precision U(NX,NY), the finite element coefficients. c c Input, double precision, external EXACT(X,Y), returns the value of c the exact solution at the point (X,Y). c c Output, double precision E2, the estimated L2 norm of the error. c implicit none integer nx integer ny integer quad_num parameter ( quad_num = 3 ) double precision abscissa(quad_num) integer e double precision e2 double precision eq integer ex double precision exact external exact integer ey integer n integer ne integer nw integer qx integer qy integer s integer se integer sw double precision u(nx,ny) double precision uq double precision vne double precision vnw double precision vse double precision vsw integer w double precision weight(quad_num) double precision wq double precision x(nx) double precision xe double precision xq double precision xw double precision y(ny) double precision yn double precision yq double precision ys save abscissa save weight data abscissa / & -0.774596669241483377035853079956D+00, & 0.000000000000000000000000000000D+00, & 0.774596669241483377035853079956D+00 / data weight / & 0.555555555555555555555555555556D+00, & 0.888888888888888888888888888889D+00, & 0.555555555555555555555555555556D+00 / e2 = 0.0D+00 c c Integrate over each interval. c do ex = 1, nx - 1 w = ex e = ex + 1 xw = x(w) xe = x(e) do ey = 1, ny - 1 s = ey n = ey + 1 ys = y(s) yn = y(n) sw = ( ey - 1 ) * nx + ex se = ( ey - 1 ) * nx + ex + 1 nw = ey * nx + ex ne = ey * nx + ex + 1 do qx = 1, quad_num xq = ( ( 1.0D+00 - abscissa(qx) ) * xw & + ( 1.0D+00 + abscissa(qx) ) * xe ) & / 2.0D+00 do qy = 1, quad_num yq = ( ( 1.0D+00 - abscissa(qy) ) * ys & + ( 1.0D+00 + abscissa(qy) ) * yn ) & / 2.0D+00 wq = weight(qx) * ( xe - xw ) / 2.0D+00 & * weight(qy) * ( yn - ys ) / 2.0D+00 vsw = ( xe - xq ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vse = ( xq - xw ) / ( xe - xw ) & * ( yn - yq ) / ( yn - ys ) vnw = ( xe - xq ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) vne = ( xq - xw ) / ( xe - xw ) & * ( yq - ys ) / ( yn - ys ) c c Note that the south-west component of U is stored in U(W,S), not U(S,W)! c uq = u(w,s) * vsw + u(e,s) * vse & + u(w,n) * vnw + u(e,n) * vne eq = exact ( xq, yq ) e2 = e2 + wq * ( uq - eq ) ** 2 end do end do end do end do e2 = sqrt ( e2 ) return end subroutine r8mat_solve2 ( n, a, b, x, ierror ) c*********************************************************************72 c cc R8MAT_SOLVE2 computes the solution of an N by N linear system. c c Discussion: c c An R8MAT is an array of R8 values. c c The linear system may be represented as c c A*X = B c c If the linear system is singular, but consistent, then the routine will c still produce a solution. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 19 August 2010 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of equations. c c Input/output, double precision A(N,N). c On input, A is the coefficient matrix to be inverted. c On output, A has been overwritten. c c Input/output, double precision B(N). c On input, B is the right hand side of the system. c On output, B has been overwritten. c c Output, double precision X(N), the solution of the linear system. c c Output, integer IERROR. c 0, no error detected. c 1, consistent singularity. c 2, inconsistent singularity. c implicit none integer n double precision a(n,n) double precision amax double precision b(n) integer i integer ierror integer imax integer ipiv(n) integer j integer k double precision x(n) ierror = 0 do i = 1, n ipiv(i) = 0 end do do i = 1, n x(i) = 0.0D+00 end do c c Process the matrix. c do k = 1, n c c In column K: c Seek the row IMAX with the properties that: c IMAX has not already been used as a pivot; c A(IMAX,K) is larger in magnitude than any other candidate. c amax = 0.0D+00 imax = 0 do i = 1, n if ( ipiv(i) .eq. 0 ) then if ( amax .lt. abs ( a(i,k) ) ) then imax = i amax = abs ( a(i,k) ) end if end if end do c c If you found a pivot row IMAX, then, c eliminate the K-th entry in all rows that have not been used for pivoting. c if ( imax .ne. 0 ) then ipiv(imax) = k do j = k + 1, n a(imax,j) = a(imax,j) / a(imax,k) end do b(imax) = b(imax) / a(imax,k) a(imax,k) = 1.0D+00 do i = 1, n if ( ipiv(i) .eq. 0 ) then do j = k + 1, n a(i,j) = a(i,j) - a(i,k) * a(imax,j) end do b(i) = b(i) - a(i,k) * b(imax) a(i,k) = 0.0D+00 end if end do end if end do c c Now, every row with nonzero IPIV begins with a 1, and c all other rows are all zero. Begin solution. c do j = n, 1, -1 imax = 0 do k = 1, n if ( ipiv(k) .eq. j ) then imax = k end if end do if ( imax .eq. 0 ) then x(j) = 0.0D+00 if ( b(j) .eq. 0.0D+00 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_SOLVE2 - Warning:' write ( *, '(a,i8)' ) & ' Consistent singularity, equation = ', j else ierror = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_SOLVE2 - Error:' write ( *, '(a,i8)' ) & ' Inconsistent singularity, equation = ', j end if else x(j) = b(imax) do i = 1, n if ( i .ne. imax ) then b(i) = b(i) - a(i,j) * x(j) end if end do end if end do return end subroutine r8vec_even ( n, alo, ahi, a ) c*********************************************************************72 c cc R8VEC_EVEN returns an R8VEC of evenly spaced values. c c Discussion: c c An R8VEC is a vector of R8 values. c c If N is 1, then the midpoint is returned. c c Otherwise, the two endpoints are returned, and N-2 evenly c spaced points between them. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 09 December 2004 c c Author: c c John Burkardt c c Parameters: c c Input, integer N, the number of values. c c Input, double precision ALO, AHI, the low and high values. c c Output, double precision A(N), N evenly spaced values. c Normally, A(1) = ALO and A(N) = AHI. c However, if N = 1, then A(1) = 0.5*(ALO+AHI). c implicit none integer n double precision a(n) double precision ahi double precision alo integer i if ( n .eq. 1 ) then a(1) = 0.5D+00 * ( alo + ahi ) else do i = 1, n a(i) = ( dble ( n - i ) * alo & + dble ( i - 1 ) * ahi ) & / dble ( n - 1 ) end do end if return end subroutine timestamp ( ) c*********************************************************************72 c cc TIMESTAMP prints out the current YMDHMS date as a timestamp. c c Licensing: c c This code is distributed under the MIT license. c c Modified: c c 12 January 2007 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