program main !*****************************************************************************80 ! !! fem2d_poisson_rectangle() solves the Poisson equation on a rectangle. ! ! Discussion: ! ! This program solves ! ! -Laplacian U(X,Y) = F(X,Y) ! ! in a rectangular region in the plane. ! ! Along the boundary of the region, Dirichlet conditions ! are imposed: ! ! U(X,Y) = G(X,Y) ! ! The code uses continuous piecewise quadratic basis functions on ! triangles determined by a uniform grid of NX by NY points. ! ! Local: ! ! real ( kind = rk ) A(3*IB+1,NUNK), the coefficient matrix. ! ! real ( kind = rk ) EH1, the H1 seminorm error. ! ! real ( kind = rk ) EL2, the L2 error. ! ! real ( kind = rk ) ELEMENT_AREA(ELEMENT_NUM), the area of elements. ! ! integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! integer ELEMENT_NUM, the number of elements. ! ! real ( kind = rk ) F(NUNK), the right hand side. ! ! integer IB, the half-bandwidth of the matrix. ! ! integer INDX(NODE_NUM), gives the index of the ! unknown quantity associated with the given node. ! ! integer NNODES, the number of nodes used to form ! one element. ! ! integer NQ, the number of quadrature points used ! for assembly. ! ! integer NQE, the number of quadrature points used ! for error estimation. ! ! integer NUNK, the number of unknowns. ! ! integer NX, the number of points in the X direction. ! ! integer NY, the number of points in the Y direction. ! ! real ( kind = rk ) WQ(NQ), quadrature weights. ! ! real ( kind = rk ) NODE_XY(2,NODE_NUM), the ! coordinates of nodes. ! ! real ( kind = rk ) XL, XR, YB, YT, the X coordinates of ! the left and right sides of the rectangle, and the Y coordinates ! of the bottom and top of the rectangle. ! ! real ( kind = rk ) XQ(NQ,ELEMENT_NUM), YQ(NQ,ELEMENT_NUM), the ! coordinates of the quadrature points in each element. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: nnodes = 6 integer, parameter :: nq = 3 integer, parameter :: nqe = 13 integer, parameter :: nx = 7 integer, parameter :: ny = 7 integer, parameter :: element_num = ( nx - 1 ) * ( ny - 1 ) * 2 integer, parameter :: node_num = ( 2 * nx - 1 ) * ( 2 * ny - 1 ) real ( kind = rk ), allocatable, dimension (:,:) :: a logical, parameter :: debug = .true. real ( kind = rk ) eh1 real ( kind = rk ) el2 real ( kind = rk ), dimension(element_num) :: element_area integer, dimension(nnodes,element_num) :: element_node real ( kind = rk ), allocatable, dimension (:) :: f integer ib integer ierr integer, dimension(node_num) :: indx integer job character ( len = 255 ) node_eps_file_name character ( len = 255 ) node_txt_file_name logical node_label integer node_show real ( kind = rk ), dimension(2,node_num) :: node_xy integer nunk integer, allocatable, dimension (:) :: pivot character ( len = 255 ) solution_txt_file_name integer triangle_show character ( len = 255 ) triangulation_eps_file_name character ( len = 255 ) triangulation_txt_file_name real ( kind = rk ), dimension(nq) :: wq real ( kind = rk ), parameter :: xl = 0.0D+00 real ( kind = rk ), dimension(nq,element_num) :: xq real ( kind = rk ), parameter :: xr = 1.0D+00 real ( kind = rk ), parameter :: yb = 0.0D+00 real ( kind = rk ), dimension(nq,element_num) :: yq real ( kind = rk ), parameter :: yt = 1.0D+00 call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle()' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Solution of the Poisson equation on a unit rectangle' write ( *, '(a)' ) ' in 2 dimensions.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' - Uxx - Uyy = F(x,y) in the box' write ( *, '(a)' ) ' U(x,y) = G(x,y) on the boundary.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The finite element method is used, with piecewise' write ( *, '(a)' ) ' quadratic basis functions on 6 node triangular' write ( *, '(a)' ) ' elements.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The corner nodes of the triangles are generated by an' write ( *, '(a)' ) ' underlying grid whose dimensions are' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' NX = ', nx write ( *, '(a,i8)' ) ' NY = ', ny write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a,i8)' ) ' Number of elements = ', element_num ! ! Set the coordinates of the nodes. ! call xy_set ( nx, ny, node_num, xl, xr, yb, yt, node_xy ) ! ! Organize the nodes into a grid of 6-node triangles. ! call grid_t6 ( nx, ny, nnodes, element_num, element_node ) ! ! Set the quadrature rule for assembly. ! call quad_a ( node_xy, element_node, element_num, node_num, & nnodes, wq, xq, yq ) ! ! Determine the areas of the elements. ! call area_set ( node_num, node_xy, nnodes, element_num, & element_node, element_area ) ! ! Determine which nodes have an associated finite element unknown. ! call indx_set ( nx, ny, node_num, indx, nunk ) write ( *, '(a,i8)' ) ' Number of unknowns = ', nunk ! ! Determine the bandwidth of the coefficient matrix. ! call bandwidth ( nnodes, element_num, element_node, node_num, indx, ib ) write ( *, '(a,i8)' ) ' The matrix bandwidth is ', 3 * ib + 1 ! ! Make an EPS picture of the nodes. ! if ( nx <= 10 .and. ny <= 10 ) then node_eps_file_name = 'rectangle_nodes.eps' node_label = .true. call nodes_plot ( node_eps_file_name, node_num, node_xy, node_label ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Wrote an EPS file' write ( *, '(a)' ) ' "' // trim ( node_eps_file_name ) // '"' write ( *, '(a)' ) ' containing a picture of the nodes.' end if ! ! Write the nodes to an ASCII file that can be read into MATLAB. ! node_txt_file_name = 'rectangle_nodes.txt' call nodes_write ( node_num, node_xy, node_txt_file_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Wrote an ASCII node file' write ( *, '(a)' ) ' "' // trim ( node_txt_file_name ) // '"' write ( *, '(a)' ) ' of the form' write ( *, '(a)' ) ' X(I), Y(I)' write ( *, '(a)' ) ' which can be used for plotting.' ! ! Make an EPS picture of the elements. ! if ( nx <= 10 .and. ny <= 10 ) then triangulation_eps_file_name = 'rectangle_elements.eps' node_show = 2 triangle_show = 2 call triangulation_order6_plot ( triangulation_eps_file_name, node_num, & node_xy, element_num, element_node, node_show, triangle_show ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Wrote an EPS file' write ( *, '(a)' ) ' "' // trim ( triangulation_eps_file_name ) // '"' write ( *, '(a)' ) ' containing a picture of the elements.' end if ! ! Write the elements to a file that can be read into MATLAB. ! triangulation_txt_file_name = 'rectangle_elements.txt' call element_write ( nnodes, element_num, element_node, & triangulation_txt_file_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Wrote an ASCII element file' write ( *, '(a)' ) ' "' // trim ( triangulation_txt_file_name ) // '"' write ( *, '(a)' ) ' of the form' write ( *, '(a)' ) ' Node(1) Node(2) Node(3) Node(4) Node(5) Node(6)' write ( *, '(a)' ) ' which can be used for plotting.' ! ! Allocate space for the coefficient matrix A and right hand side F. ! allocate ( a(3*ib+1,nunk) ) allocate ( f(nunk) ) allocate ( pivot(nunk) ) ! ! Assemble the coefficient matrix A and the right-hand side F of the ! finite element equations. ! call assemble ( node_num, node_xy, nnodes, element_num, & element_node, nq, wq, xq, yq, element_area, indx, ib, nunk, a, f ) if ( debug ) then call dgb_print_some ( nunk, nunk, ib, ib, a, 1, 1, 5, 5, & ' Initial 5 x 5 block of matrix A:' ) call r8vec_print_some ( nunk, f, 1, 10, ' Part of right hand side F:' ) end if ! ! Modify the coefficient matrix and right hand side to account for ! boundary conditions. ! call boundary ( nx, ny, node_num, node_xy, indx, ib, nunk, a, f ) if ( debug ) then call dgb_print_some ( nunk, nunk, ib, ib, a, 1, 1, 5, 5, & ' A after BC adjustment:' ) call r8vec_print_some ( nunk, f, 1, 10, ' F after BC adjustment:' ) end if ! ! Solve the linear system using a banded solver. ! call dgb_fa ( nunk, ib, ib, a, pivot, ierr ) if ( ierr /= 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, and the' write ( *, '(a)' ) ' algorithm cannot proceed.' stop end if job = 0 call dgb_sl ( nunk, ib, ib, a, pivot, f, job ) call r8vec_print_some ( nunk, f, 1, 10, ' Part of the solution vector:' ) ! ! Calculate error using 13 point quadrature rule. ! call errors ( element_area, element_node, indx, node_xy, f, & element_num, nnodes, nqe, nunk, node_num, el2, eh1 ) call compare ( node_num, node_xy, indx, nunk, f ) ! ! Write an ASCII file that can be read into MATLAB. ! solution_txt_file_name = 'rectangle_solution.txt' call solution_write ( f, indx, node_num, nunk, & solution_txt_file_name, node_xy ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Wrote an ASCII solution file' write ( *, '(a)' ) ' "' // trim ( solution_txt_file_name ) // '"' write ( *, '(a)' ) ' of the form' write ( *, '(a)' ) ' U( X(I), Y(I) )' write ( *, '(a)' ) ' which can be used for plotting.' ! ! Free memory. ! deallocate ( a ) deallocate ( f ) deallocate ( pivot ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'fem2d_poisson_rectangle():' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop ( 0 ) end subroutine area_set ( node_num, node_xy, nnodes, element_num, & element_node, element_area ) !*****************************************************************************80 ! !! AREA_SET sets the area of each element. ! ! Discussion: ! ! The areas of the elements are needed in order to adjust ! the integral estimates produced by the quadrature formulas. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the ! coordinates of the nodes. ! ! Input, integer NNODES, the number of local nodes per element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Output, real ( kind = rk ) ELEMENT_AREA(ELEMENT_NUM), the area of elements. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer node_num integer element real ( kind = rk ) element_area(element_num) integer element_node(nnodes,element_num) integer i1 integer i2 integer i3 real ( kind = rk ) x1 real ( kind = rk ) x2 real ( kind = rk ) x3 real ( kind = rk ) node_xy(2,node_num) real ( kind = rk ) y1 real ( kind = rk ) y2 real ( kind = rk ) y3 do element = 1, element_num i1 = element_node(1,element) x1 = node_xy(1,i1) y1 = node_xy(2,i1) i2 = element_node(2,element) x2 = node_xy(1,i2) y2 = node_xy(2,i2) i3 = element_node(3,element) x3 = node_xy(1,i3) y3 = node_xy(2,i3) element_area(element) = 0.5D+00 * abs & ( y1 * ( x2 - x3 ) & + y2 * ( x3 - x1 ) & + y3 * ( x1 - x2 ) ) end do return end subroutine assemble ( node_num, node_xy, nnodes, element_num, & element_node, nq, wq, xq, yq, element_area, indx, ib, nunk, a, f ) !*****************************************************************************80 ! !! ASSEMBLE assembles the coefficient matrix A and right hand side F. ! ! Discussion: ! ! The matrix is known to be banded. A special matrix storage format ! is used to reduce the space required. Details of this format are ! discussed in the routine DGB_FA. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of nodes. ! ! Input, integer NNODES, the number of nodes used to form ! one element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer NQ, the number of quadrature points used ! in assembly. ! ! Input, real ( kind = rk ) WQ(NQ), quadrature weights. ! ! Input, real ( kind = rk ) XQ(NQ,ELEMENT_NUM), YQ(NQ,ELEMENT_NUM), the ! coordinates of the quadrature points in each element. ! ! Input, real ( kind = rk ) ELEMENT_AREA(ELEMENT_NUM), the area of elements. ! ! Input, integer INDX(NODE_NUM), gives the index of the ! unknown quantity associated with the given node. ! ! Input, integer IB, the half-bandwidth of the matrix. ! ! Input, integer NUNK, the number of unknowns. ! ! Output, real ( kind = rk ) A(3*IB+1,NUNK), the NUNK by NUNK ! coefficient matrix, stored in a compressed format. ! ! Output, real ( kind = rk ) F(NUNK), the right hand side. ! ! Local parameters: ! ! Local, real ( kind = rk ) BB, BX, BY, the value of some basis function ! and its first derivatives at a quadrature point. ! ! Local, real ( kind = rk ) BBB, BBX, BBY, the value of another basis ! function and its first derivatives at a quadrature point. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ib integer element_num integer nnodes integer node_num integer nq integer nunk real ( kind = rk ), dimension(3*ib+1,nunk) :: a real ( kind = rk ) aij integer basis real ( kind = rk ) bi real ( kind = rk ) bj real ( kind = rk ) dbidx real ( kind = rk ) dbidy real ( kind = rk ) dbjdx real ( kind = rk ) dbjdy integer element real ( kind = rk ), dimension(element_num) :: element_area integer, dimension(nnodes,element_num) :: element_node real ( kind = rk ), dimension(nunk) :: f integer i integer, dimension(node_num) :: indx integer ip integer ipp integer j real ( kind = rk ), dimension(2,node_num) :: node_xy integer quad real ( kind = rk ) :: rhs integer test real ( kind = rk ) w real ( kind = rk ), dimension(nq) :: wq real ( kind = rk ) x real ( kind = rk ), dimension(nq,element_num) :: xq real ( kind = rk ) y real ( kind = rk ), dimension(nq,element_num) :: yq ! ! Initialize the arrays to zero. ! f(1:nunk) = 0.0D+00 a(1:3*ib+1,1:nunk) = 0.0D+00 ! ! The actual values of A and F are determined by summing up ! contributions from all the elements. ! do element = 1, element_num ! ! Consider a quadrature point QUAD, with coordinates (X,Y). ! do quad = 1, nq x = xq(quad,element) y = yq(quad,element) w = element_area(element) * wq(quad) ! ! Consider one of the basis functions, which will play the ! role of test function in the integral. ! ! We will generate an integral for the test function; that is, a basis ! function associated with a degree of freedom. ! ! If the degree of freedom associated with the node is ! constrained by a boundary condition, then the finite element ! integral we set up here will be modified or replaced later ! (see subroutine BOUNDARY). ! do test = 1, nnodes ip = element_node(test,element) i = indx(ip) call qbf ( x, y, element, test, node_xy, element_node, & element_num, nnodes, node_num, bi, dbidx, dbidy ) f(i) = f(i) + w * rhs ( x, y ) * bi ! ! Consider a basis function, which is used to form the ! value of the solution function. ! ! If this basis function is associated with a boundary condition, ! then subtract the term from the right hand side. ! ! Otherwise, this term is included in the system matrix. ! ! Logically, this term goes in entry A(I,J). Because of the ! band matrix storage, entry (I,J) is actually stored in ! A(I-J+2*NHBA+1,J). ! do basis = 1, nnodes ipp = element_node(basis,element) j = indx(ipp) call qbf ( x, y, element, basis, node_xy, element_node, & element_num, nnodes, node_num, bj, dbjdx, dbjdy ) aij = dbidx * dbjdx + dbidy * dbjdy a(i-j+2*ib+1,j) = a(i-j+2*ib+1,j) + w * aij end do end do end do end do return end subroutine bandwidth ( nnodes, element_num, element_node, node_num, indx, nhba ) !*****************************************************************************80 ! !! BANDWIDTH determines the bandwidth of the coefficient matrix. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NNODES, the number of local nodes per element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer INDX(NODE_NUM), the index of the unknown in the ! finite element linear system. ! ! Output, integer NHBA, the half bandwidth of the matrix. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer node_num integer element integer element_node(nnodes,element_num) integer i integer iln integer in integer indx(node_num) integer j integer jln integer jn integer nhba nhba = 0 do element = 1, element_num do iln = 1, nnodes in = element_node(iln,element) i = indx(in) if ( 0 < i ) then do jln = 1, nnodes jn = element_node(jln,element) j = indx(jn) nhba = max ( nhba, j - i ) end do end if end do end do return end subroutine boundary ( nx, ny, node_num, node_xy, indx, ib, nunk, a, f ) !*****************************************************************************80 ! !! BOUNDARY modifies the linear system for boundary conditions. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NX, NY, controls the number of elements along ! the X and Y directions. The number of elements will be ! 2 * ( NX - 1 ) * ( NY - 1 ). ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of nodes. ! ! Input, integer INDX(NODE_NUM), gives the index of the ! unknown quantity associated with the given node. ! ! Input, integer IB, the half-bandwidth of the matrix. ! ! Input, integer NUNK, the number of unknowns. ! ! Input/output, real ( kind = rk ) A(3*IB+1,NUNK), the NUNK by NUNK ! coefficient matrix, stored in a compressed format. On output, ! A has been adjusted for boundary conditions. ! ! Input/output, real ( kind = rk ) F(NUNK), the right hand side. ! On output, F has been adjusted for boundary conditions. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ib integer node_num integer nunk real ( kind = rk ), dimension(3*ib+1,nunk) :: a integer col real ( kind = rk ) dudx real ( kind = rk ) dudy real ( kind = rk ), dimension(nunk) :: f integer i integer, dimension(node_num) :: indx integer j integer jhi integer jlo integer node real ( kind = rk ), dimension(2,node_num) :: node_xy integer nx integer ny integer row real ( kind = rk ) u real ( kind = rk ) x real ( kind = rk ) y ! ! Consider each node. ! node = 0 do row = 1, 2 * ny - 1 do col = 1, 2 * nx - 1 node = node + 1 if ( row == 1 .or. & row == 2 * ny - 1 .or. & col == 1 .or. & col == 2 * nx - 1 ) then i = indx(node) x = node_xy(1,node) y = node_xy(2,node) call exact ( x, y, u, dudx, dudy ) jlo = max ( i - ib, 1 ) jhi = min ( i + ib, nunk ) do j = jlo, jhi a(i-j+2*ib+1,j) = 0.0D+00 end do a(i-i+2*ib+1,i) = 1.0D+00 f(i) = u end if end do end do return end subroutine compare ( node_num, node_xy, indx, nunk, f ) !*****************************************************************************80 ! !! COMPARE compares the exact and computed solution at the nodes. ! ! Discussion: ! ! This is a rough comparison, done only at the nodes. Such a pointwise ! comparison is easy, because the value of the finite element ! solution is exactly the value of the finite element coefficient ! associated with that node. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the ! coordinates of the nodes. ! ! Input, integer INDX(NODE_NUM), the index of the unknown in the ! finite element linear system. ! ! Input, integer NUNK, the number of unknowns in the finite ! element system. ! ! Input, real ( kind = rk ) F(NUNK), the solution vector of the finite ! element system. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer nunk real ( kind = rk ) dudx real ( kind = rk ) dudy real ( kind = rk ) f(nunk) integer i integer indx(node_num) integer node real ( kind = rk ) node_xy(2,node_num) real ( kind = rk ) u real ( kind = rk ) uh real ( kind = rk ) x real ( kind = rk ) y write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'COMPARE:' write ( *, '(a)' ) ' Compare computed and exact solutions at the nodes.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' X Y U U' write ( *, '(a)' ) ' computed exact' write ( *, '(a)' ) ' ' do node = 1, node_num x = node_xy(1,node) y = node_xy(2,node) call exact ( x, y, u, dudx, dudy ) i = indx(node) uh = f(i) write ( *, '(2x,4f12.4)' ) x, y, uh, u end do return end subroutine dgb_fa ( n, ml, mu, a, pivot, info ) !*****************************************************************************80 ! !! DGB_FA performs a LINPACK-style PLU factorization of an DGB matrix. ! ! Discussion: ! ! The DGB storage format is for an M by N banded matrix, with lower ! bandwidth ML and upper bandwidth MU. Storage includes room for ML ! extra superdiagonals, which may be required to store nonzero entries ! generated during Gaussian elimination. ! ! The original M by N matrix is "collapsed" downward, so that diagonals ! become rows of the storage array, while columns are preserved. The ! collapsed array is logically 2*ML+MU+1 by N. ! ! 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 ML+MU by ML+MU upper left triangle and the ! ML by ML lower right triangle are not referenced. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979 ! ! Parameters: ! ! 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/output, real ( kind = rk ) A(2*ML+MU+1,N), on input, the matrix ! in band storage, on output, information about the LU factorization. ! ! Output, integer PIVOT(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, parameter :: rk = kind ( 1.0D+00 ) integer ml integer mu integer n real ( kind = rk ) a(2*ml+mu+1,n) 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 real ( kind = rk ) temp 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 a(i0:ml,jz) = 0.0D+00 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 a(1:ml,jz) = 0.0D+00 end if ! ! Find L = pivot index. ! lm = min ( ml, n-k ) l = m do j = m+1, m+lm if ( abs ( a(l,k) ) < abs ( a(j,k) ) ) then l = j end if end do pivot(k) = l + k - m ! ! Zero pivot implies this column already triangularized. ! if ( a(l,k) == 0.0D+00 ) then info = k write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DGB_FA - Fatal error!' write ( *, '(a,i8)' ) ' Zero pivot on step ', info return end if ! ! Interchange if necessary. ! if ( m /= l ) then temp = a(l,k) a(l,k) = a(m,k) a(m,k) = temp end if ! ! Compute multipliers. ! a(m+1:m+lm,k) = - a(m+1:m+lm,k) / a(m,k) ! ! Row elimination with column indexing. ! 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 /= mm ) then temp = a(l,j) a(l,j) = a(mm,j) a(mm,j) = temp end if a(mm+1:mm+lm,j) = a(mm+1:mm+lm,j) + a(mm,j) * a(m+1:m+lm,k) end do end do pivot(n) = n if ( a(m,n) == 0.0D+00 ) then info = n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DGB_FA - Fatal error!' write ( *, '(a,i8)' ) ' Zero pivot on step ', info end if return end subroutine dgb_print_some ( m, n, ml, mu, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! DGB_PRINT_SOME prints some of a DGB matrix. ! ! Discussion: ! ! The DGB storage format is for an M by N banded matrix, with lower ! bandwidth ML and upper bandwidth MU. Storage includes room for ML ! extra superdiagonals, which may be required to store nonzero entries ! generated during Gaussian elimination. ! ! The original M by N matrix is "collapsed" downward, so that diagonals ! become rows of the storage array, while columns are preserved. The ! collapsed array is logically 2*ML+MU+1 by N. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! 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.. ! ! Input, real ( kind = rk ) A(2*ML+MU+1,N), the DGB matrix. ! ! Input, integer ILO, JLO, IHI, JHI, the first row and ! column, and the last row and column to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer ml integer mu integer n real ( kind = rk ) a(2*ml+mu+1,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 ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if ! ! Print the columns of the matrix, in strips of 5. ! 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)' ) ' ---' ! ! Determine the range of the rows in this strip. ! i2lo = max ( ilo, 1 ) i2lo = max ( i2lo, j2lo - mu - ml ) i2hi = min ( ihi, m ) i2hi = min ( i2hi, j2hi + ml ) do i = i2lo, i2hi ! ! Print out (up to) 5 entries in row I, that lie in the current strip. ! do j2 = 1, inc j = j2lo - 1 + j2 if ( i < j - ml - mu .or. j + ml < 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 ( n, ml, mu, a, pivot, b, job ) !*****************************************************************************80 ! !! DGB_SL solves a system factored by DGB_FA. ! ! Discussion: ! ! The DGB storage format is for an M by N banded matrix, with lower ! bandwidth ML and upper bandwidth MU. Storage includes room for ML ! extra superdiagonals, which may be required to store nonzero entries ! generated during Gaussian elimination. ! ! The original M by N matrix is "collapsed" downward, so that diagonals ! become rows of the storage array, while columns are preserved. The ! collapsed array is logically 2*ML+MU+1 by N. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979 ! ! Parameters: ! ! 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, real ( kind = rk ) A(2*ML+MU+1,N), the LU factors from DGB_FA. ! ! Input, integer PIVOT(N), the pivot vector from DGB_FA. ! ! Input/output, real ( kind = rk )l B(N). ! On input, the right hand side vector. ! On output, the solution. ! ! Input, integer JOB. ! 0, solve A * x = b. ! nonzero, solve A' * x = b. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ml integer mu integer n real ( kind = rk ) a(2*ml+mu+1,n) real ( kind = rk ) b(n) integer pivot(n) integer job integer k integer l integer la integer lb integer lm integer m real ( kind = rk ) temp m = mu + ml + 1 ! ! Solve A * x = b. ! if ( job == 0 ) then ! ! Solve L * Y = B. ! if ( 1 <= ml ) then do k = 1, n-1 lm = min ( ml, n-k ) l = pivot(k) if ( l /= k ) then temp = b(l) b(l) = b(k) b(k) = temp end if b(k+1:k+lm) = b(k+1:k+lm) + b(k) * a(m+1:m+lm,k) 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 b(lb:lb+lm-1) = b(lb:lb+lm-1) - b(k) * a(la:la+lm-1,k) end do ! ! Solve A' * X = B. ! else ! ! Solve U' * Y = B. ! do k = 1, n lm = min ( k, m ) - 1 la = m - lm lb = k - lm b(k) = ( b(k) - sum ( a(la:la+lm-1,k) * b(lb:lb+lm-1) ) ) & / a(m,k) end do ! ! Solve L' * X = Y. ! if ( 1 <= ml ) then do k = n-1, 1, -1 lm = min ( ml, n-k ) b(k) = b(k) + sum ( a(m+1:m+lm,k) * b(k+1:k+lm) ) l = pivot(k) if ( l /= k ) then temp = b(l) b(l) = b(k) b(k) = temp end if end do end if end if return end subroutine element_write ( nnodes, element_num, element_node, file_name ) !*****************************************************************************80 ! !! ELEMENT_WRITE writes the element information to a text file. ! ! Discussion: ! ! The element information and the node/solution information can be read ! into another program to make plots. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NNODES, the number of local nodes per element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer element integer element_node(nnodes,element_num) character ( len = * ) file_name integer output_status integer output_unit call get_unit ( output_unit ) open ( unit = output_unit, file = file_name, status = 'replace', & iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENTS_WRITE - Warning!' write ( *, '(a)' ) ' Could not write the element file.' return end if do element = 1, element_num write ( output_unit, '(2x,i4,2x,i4,2x,i4,2x,i4,2x,i4,2x,i4)' ) & element_node(1:nnodes,element) end do close ( unit = output_unit ) return end subroutine errors ( element_area, element_node, indx, node_xy, f, & element_num, nnodes, nqe, nunk, node_num, el2, eh1 ) !*****************************************************************************80 ! !! ERRORS calculates the error in the L2 norm and H1 seminorm. ! ! Discussion: ! ! This routine uses a 13 point quadrature rule in each element, ! in order to estimate the values of ! ! EL2 = Sqrt ( Integral ( U(x,y) - Uh(x,y) )**2 dx dy ) ! ! EH1 = Sqrt ( Integral ( Ux(x,y) - Uhx(x,y) )**2 + ! ( Uy(x,y) - Uhy(x,y) )**2 dx dy ) ! ! Here U is the exact solution, and Ux and Uy its spatial derivatives, ! as evaluated by a user-supplied routine. ! ! Uh, Uhx and Uhy are the computed solution and its spatial derivatives, ! as specified by the computed finite element solution. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) ELEMENT_AREA(ELEMENT_NUM), the area of elements. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer INDX(NODE_NUM), the index of the unknown ! quantity associated with the given node. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the X and Y ! coordinates of nodes. ! ! Input, real ( kind = rk ) F(NUNK), the coefficients of the solution. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer NNODES, the number of nodes used to form ! one element. ! ! Input, integer NQE, the number of points in the quadrature ! rule. This is actually fixed at 13. ! ! Input, integer NUNK, the number of unknowns. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Output, real ( kind = rk ) EL2, the L2 error. ! ! Output, real ( kind = rk ) EH1, the H1 seminorm error. ! ! Local Parameters: ! ! Local, real ( kind = rk ) AR, the weight for a given quadrature point ! in a given element. ! ! Local, real ( kind = rk ) BB, BX, BY, a basis function and its first ! derivatives evaluated at a particular quadrature point. ! ! Local, real ( kind = rk ) EH1, the H1 seminorm error. ! ! Local, real ( kind = rk ) EL2, the L2 error. ! ! Local, real ( kind = rk ) UEX, UEXX, UEXY, the exact solution and its first ! derivatives evaluated at a particular quadrature point. ! ! Local, real ( kind = rk ) UH, UHX, UHY, the computed solution and its first ! derivatives evaluated at a particular quadrature point. ! ! Local, real ( kind = rk ) WQE(NQE), stores the quadrature weights. ! ! Local, real ( kind = rk ) X, Y, the coordinates of a particular ! quadrature point. ! ! Local, real ( kind = rk ) XQE(NQE), YQE(NQE), stores the location ! of quadrature points in a given element. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer node_num integer nqe integer nunk real ( kind = rk ) :: ar real ( kind = rk ) bi real ( kind = rk ) dbidx real ( kind = rk ) dbidy real ( kind = rk ) :: dudx real ( kind = rk ) :: dudxh real ( kind = rk ) :: dudy real ( kind = rk ) :: dudyh real ( kind = rk ) :: eh1 real ( kind = rk ) :: el2 integer element real ( kind = rk ), dimension(element_num) :: element_area integer, dimension(nnodes,element_num) :: element_node real ( kind = rk ), dimension(nunk) :: f integer i integer in1 integer, dimension(node_num) :: indx integer ip real ( kind = rk ), dimension(2,node_num) :: node_xy integer quad real ( kind = rk ) :: u real ( kind = rk ) :: uh real ( kind = rk ), dimension(nqe) :: wqe real ( kind = rk ) :: x real ( kind = rk ), dimension(nqe) :: xqe real ( kind = rk ) :: y real ( kind = rk ), dimension(nqe) :: yqe el2 = 0.0D+00 eh1 = 0.0D+00 ! ! For each element, retrieve the nodes, area, quadrature weights, ! and quadrature points. ! do element = 1, element_num call quad_e ( node_xy, element_node, element, & element_num, nnodes, node_num, nqe, wqe, xqe, yqe ) ! ! For each quadrature point, evaluate the computed solution and its X and ! Y derivatives. ! do quad = 1, nqe ar = element_area(element) * wqe(quad) x = xqe(quad) y = yqe(quad) uh = 0.0D+00 dudxh = 0.0D+00 dudyh = 0.0D+00 do in1 = 1, nnodes ip = element_node(in1,element) call qbf ( x, y, element, in1, node_xy, element_node, & element_num, nnodes, node_num, bi, dbidx, dbidy ) i = indx(ip) uh = uh + bi * f(i) dudxh = dudxh + dbidx * f(i) dudyh = dudyh + dbidy * f(i) end do ! ! Evaluate the exact solution and its X and Y derivatives. ! call exact ( x, y, u, dudx, dudy ) ! ! Add the weighted value at this quadrature point to the quadrature sum. ! el2 = el2 + ( uh - u )**2 * ar eh1 = eh1 + ( ( dudxh - dudx )**2 + ( dudyh - dudy )**2 ) * ar end do end do el2 = sqrt ( el2 ) eh1 = sqrt ( eh1 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' *********************************************' write ( *, '(a)' ) ' * *' write ( *, '(a)' ) ' * ERRORS: *' 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, u, dudx, dudy ) !*****************************************************************************80 ! !! EXACT calculates the exact solution and its first derivatives. ! ! Discussion: ! ! The function specified here depends on the problem being ! solved. The user must be sure to change both EXACT and RHS ! or the program will have inconsistent data. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) X, Y, the coordinates of a point ! in the region, at which the exact solution is to be evaluated. ! ! Output, real ( kind = rk ) U, DUDX, DUDY, the value of ! the exact solution U and its derivatives dUdX ! and dUdY at the point (X,Y). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) dudx real ( kind = rk ) dudy real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 real ( kind = rk ) u real ( kind = rk ) x real ( kind = rk ) y u = sin ( pi * x ) * sin ( pi * y ) + x dudx = pi * cos ( pi * x ) * sin ( pi * y ) + 1.0D+00 dudy = pi * sin ( pi * x ) * cos ( pi * y ) return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is a value between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is a value between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT, the free unit number. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine grid_t6 ( nx, ny, nnodes, element_num, element_node ) !*****************************************************************************80 ! !! GRID_T6 produces a grid of pairs of 6 node triangles. ! ! Example: ! ! Input: ! ! NX = 4, NY = 3 ! ! Output: ! ! ELEMENT_NODE = ! 1, 3, 15, 2, 9, 8; ! 17, 15, 3, 16, 9, 10; ! 3, 5, 17, 4, 11, 10; ! 19, 17, 5, 18, 11, 12; ! 5, 7, 19, 6, 13, 12; ! 21, 19, 7, 20, 13, 14; ! 15, 17, 29, 16, 23, 22; ! 31, 29, 17, 30, 23, 24; ! 17, 19, 31, 18, 25, 24; ! 33, 31, 19, 32, 25, 26; ! 19, 21, 33, 20, 27, 26; ! 35, 33, 21, 34, 27, 28. ! ! Diagram: ! ! 29-30-31-32-33-34-35 ! |\ 8 |\10 |\12 | ! | \ | \ | \ | ! 22 23 24 25 26 27 28 ! | \ | \ | \ | ! | 7 \| 9 \| 11 \| ! 15-16-17-18-19-20-21 ! |\ 2 |\ 4 |\ 6 | ! | \ | \ | \ | ! 8 9 10 11 12 13 14 ! | \ | \ | \ | ! | 1 \| 3 \| 5 \| ! 1--2--3--4--5--6--7 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NX, NY, controls the number of elements ! along the X and Y directions. The number of elements will be ! 2 * ( NX - 1 ) * ( NY - 1 ). ! ! Input, integer NNODES, the number of local nodes per element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Output, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the index of the I-th node of the J-th element. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer c integer e integer element integer element_node(nnodes,element_num) integer i integer j integer n integer ne integer nw integer nx integer ny integer s integer se integer sw integer w element = 0 do j = 1, ny - 1 do i = 1, nx - 1 sw = ( j - 1 ) * 2 * ( 2 * nx - 1 ) + 2 * i - 1 w = sw + 1 nw = sw + 2 s = sw + 2 * nx - 1 c = s + 1 n = s + 2 se = s + 2 * nx - 1 e = se + 1 ne = se + 2 element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = nw element_node(4,element) = s element_node(5,element) = c element_node(6,element) = w element = element + 1 element_node(1,element) = ne element_node(2,element) = nw element_node(3,element) = se element_node(4,element) = n element_node(5,element) = c element_node(6,element) = e end do end do return end subroutine indx_set ( nx, ny, node_num, indx, nunk ) !*****************************************************************************80 ! !! INDX_SET assigns an unknown value index at each node. ! ! Discussion: ! ! Every finite element node will is assigned an index which ! indicates the finite element basis function and its coefficient ! which are associated with that node. ! ! Example: ! ! On a simple 5 by 5 grid, where the nodes are numbered starting ! at the lower left, and increasing in X first, we would have the ! following values of INDX: ! ! 21 22 23 24 25 ! 16 17 18 19 20 ! 11 12 13 14 15 ! 6 7 8 9 10 ! 1 2 3 4 5 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NX, NY, the number of elements in the X and ! Y directions. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Output, integer INDX(NODE_NUM), the index of the unknown ! in the finite element linear system. ! ! Output, integer NUNK, the number of unknowns. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer i integer in integer indx(node_num) integer j integer nunk integer nx integer ny nunk = 0 in = 0 do j = 1, 2 * ny - 1 do i = 1, 2 * nx - 1 in = in + 1 nunk = nunk + 1 indx(in) = nunk end do end do return end subroutine nodes_plot ( file_name, node_num, node_xy, node_label ) !*****************************************************************************80 ! !! NODES_PLOT plots a pointset. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer NODE_NUM, the number of points. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the nodes. ! ! Input, logical NODE_LABEL, is TRUE if the nodes should be labeled. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer :: circle_size integer delta character ( len = * ) file_name integer file_unit integer ios integer node logical node_label real ( kind = rk ) node_xy(2,node_num) character ( len = 40 ) string real ( kind = rk ) x_max real ( kind = rk ) x_min integer x_ps integer :: x_ps_max = 576 integer :: x_ps_max_clip = 594 integer :: x_ps_min = 36 integer :: x_ps_min_clip = 18 real ( kind = rk ) x_scale real ( kind = rk ) y_max real ( kind = rk ) y_min integer y_ps integer :: y_ps_max = 666 integer :: y_ps_max_clip = 684 integer :: y_ps_min = 126 integer :: y_ps_min_clip = 108 real ( kind = rk ) y_scale ! ! We need to do some figuring here, so that we can determine ! the range of the data, and hence the height and width ! of the piece of paper. ! x_max = maxval ( node_xy(1,1:node_num) ) x_min = minval ( node_xy(1,1:node_num) ) x_scale = x_max - x_min x_max = x_max + 0.05D+00 * x_scale x_min = x_min - 0.05D+00 * x_scale x_scale = x_max - x_min y_max = maxval ( node_xy(2,1:node_num) ) y_min = minval ( node_xy(2,1:node_num) ) y_scale = y_max - y_min y_max = y_max + 0.05D+00 * y_scale y_min = y_min - 0.05D+00 * y_scale y_scale = y_max - y_min if ( x_scale < y_scale ) then delta = nint ( real ( x_ps_max - x_ps_min, kind = rk ) & * ( y_scale - x_scale ) / ( 2.0D+00 * y_scale ) ) x_ps_max = x_ps_max - delta x_ps_min = x_ps_min + delta x_ps_max_clip = x_ps_max_clip - delta x_ps_min_clip = x_ps_min_clip + delta x_scale = y_scale else if ( y_scale < x_scale ) then delta = nint ( real ( y_ps_max - y_ps_min, kind = rk ) & * ( x_scale - y_scale ) / ( 2.0D+00 * x_scale ) ) y_ps_max = y_ps_max - delta y_ps_min = y_ps_min + delta y_ps_max_clip = y_ps_max_clip - delta y_ps_min_clip = y_ps_min_clip + delta y_scale = x_scale end if call get_unit ( file_unit ) open ( unit = file_unit, file = file_name, status = 'replace', & iostat = ios ) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NODES_PLOT - Fatal error!' write ( *, '(a)' ) ' Can not open output file.' return end if write ( file_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( file_unit, '(a)' ) '%%Creator: nodes_plot.f90' write ( file_unit, '(a)' ) '%%Title: ' // trim ( file_name ) write ( file_unit, '(a)' ) '%%Pages: 1' write ( file_unit, '(a,i3,2x,i3,2x,i3,2x,i3)' ) '%%BoundingBox: ', & x_ps_min, y_ps_min, x_ps_max, y_ps_max write ( file_unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( file_unit, '(a)' ) '%%LanguageLevel: 1' write ( file_unit, '(a)' ) '%%EndComments' write ( file_unit, '(a)' ) '%%BeginProlog' write ( file_unit, '(a)' ) '/inch {72 mul} def' write ( file_unit, '(a)' ) '%%EndProlog' write ( file_unit, '(a)' ) '%%Page: 1 1' write ( file_unit, '(a)' ) 'save' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB line color to very light gray.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.900 0.900 0.900 setrgbcolor' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Draw a gray border around the page.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) 'newpath' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_min, ' moveto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_max, y_ps_min, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_max, y_ps_max, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_max, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_min, ' lineto' write ( file_unit, '(a)' ) 'stroke' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB line color to black.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.000 0.000 0.000 setrgbcolor' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the font and its size.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '/Times-Roman findfont' write ( file_unit, '(a)' ) '0.50 inch scalefont' write ( file_unit, '(a)' ) 'setfont' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Print a title.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% 210 702 moveto' write ( file_unit, '(a)' ) '% (Pointset) show' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Define a clipping polygon.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) 'newpath' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_min_clip, ' moveto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_max_clip, y_ps_min_clip, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_max_clip, y_ps_max_clip, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_max_clip, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_min_clip, ' lineto' write ( file_unit, '(a)' ) 'clip newpath' ! ! Draw the nodes. ! if ( node_num <= 200 ) then circle_size = 5 else if ( node_num <= 500 ) then circle_size = 4 else if ( node_num <= 1000 ) then circle_size = 3 else if ( node_num <= 5000 ) then circle_size = 2 else circle_size = 1 end if write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Draw filled dots at each node.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB color to blue.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.000 0.150 0.750 setrgbcolor' write ( file_unit, '(a)' ) '%' do node = 1, node_num x_ps = int ( & ( ( x_max - node_xy(1,node) ) * real ( x_ps_min, kind = rk ) & + ( node_xy(1,node) - x_min ) * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( file_unit, '(a,i4,2x,i4,2x,i4,2x,a)' ) 'newpath ', x_ps, y_ps, & circle_size, '0 360 arc closepath fill' end do ! ! Label the nodes. ! if ( node_label ) then write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Label the nodes:' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB color to darker blue.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.000 0.250 0.850 setrgbcolor' write ( file_unit, '(a)' ) '/Times-Roman findfont' write ( file_unit, '(a)' ) '0.20 inch scalefont' write ( file_unit, '(a)' ) 'setfont' do node = 1, node_num x_ps = int ( & ( ( x_max - node_xy(1,node) ) * real ( x_ps_min, kind = rk ) & + ( + node_xy(1,node) - x_min ) * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( string, '(i4)' ) node string = adjustl ( string ) write ( file_unit, '(i4,2x,i4,a)' ) x_ps, y_ps+5, & ' moveto (' // trim ( string ) // ') show' end do end if write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) 'restore showpage' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% End of page.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '%%Trailer' write ( file_unit, '(a)' ) '%%EOF' close ( unit = file_unit ) return end subroutine nodes_write ( node_num, node_xy, output_filename ) !*****************************************************************************80 ! !! NODES_WRITE writes the nodes to a file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the nodes. ! ! Input, character ( len = * ) OUTPUT_FILENAME, the name of the file ! in which the data should be stored. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer node real ( kind = rk ), dimension(2,node_num) :: node_xy character ( len = * ) :: output_filename integer output_status integer output_unit call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, status = 'replace', & iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NODES_WRITE - Warning!' write ( *, '(a)' ) ' Could not write the node file.' return end if do node = 1, node_num write ( output_unit, '(f8.4,2x,f8.4,2x,g14.6)' ) node_xy(1:2,node) end do close ( unit = output_unit ) return end subroutine qbf ( x, y, element, inode, node_xy, element_node, & element_num, nnodes, node_num, b, dbdx, dbdy ) !*****************************************************************************80 ! !! QBF evaluates the quadratic basis functions. ! ! Discussion: ! ! This routine assumes that the "midpoint" nodes are, in fact, ! exactly the average of the two extreme nodes. This is NOT true ! for a general quadratic triangular element. ! ! Assuming this property of the midpoint nodes makes it easy to ! determine the values of (R,S) in the reference element that ! correspond to (X,Y) in the physical element. ! ! Once we know the (R,S) coordinates, it's easy to evaluate the ! basis functions and derivatives. ! ! The physical element T6: ! ! In this picture, we don't mean to suggest that the bottom of ! the physical triangle is horizontal. However, we do assume that ! each of the sides is a straight line, and that the intermediate ! points are exactly halfway on each side. ! ! | ! | ! | 3 ! | / \ ! | / \ ! Y 6 5 ! | / \ ! | / \ ! | 1-----4-----2 ! | ! +--------X--------> ! ! Reference element T6: ! ! In this picture of the reference element, we really do assume ! that one side is vertical, one horizontal, of length 1. ! ! | ! | ! 1 3 ! | |\ ! | | \ ! S 6 5 ! | | \ ! | | \ ! 0 1--4--2 ! | ! +--0--R--1--------> ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) X, Y, the (global) coordinates of the point ! at which the basis function is to be evaluated. ! ! Input, integer ELEMENT, the index of the element which ! contains the point. ! ! Input, integer INODE, the local index (between 1 and 6) that ! specifies which basis function is to be evaluated. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the X and Y ! coordinates of nodes. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer NNODES, the number of nodes used to form ! one element. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Output, real ( kind = rk ) B, DBDX, DBDY, the value of the basis function ! and its X and Y derivatives at (X,Y). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer node_num real ( kind = rk ) b real ( kind = rk ) dbdr real ( kind = rk ) dbds real ( kind = rk ) dbdx real ( kind = rk ) dbdy real ( kind = rk ) det real ( kind = rk ) drdx real ( kind = rk ) drdy real ( kind = rk ) dsdx real ( kind = rk ) dsdy integer element integer, dimension(nnodes,element_num) :: element_node integer inode real ( kind = rk ), dimension(2,node_num) :: node_xy real ( kind = rk ) r real ( kind = rk ) s real ( kind = rk ) x real ( kind = rk ) xn(6) real ( kind = rk ) y real ( kind = rk ) yn(6) xn(1:6) = node_xy(1,element_node(1:6,element)) yn(1:6) = node_xy(2,element_node(1:6,element)) ! ! Determine the (R,S) coordinates corresponding to (X,Y). ! ! What is happening here is that we are solving the linear system: ! ! ( X2-X1 X3-X1 ) * ( R ) = ( X - X1 ) ! ( Y2-Y1 Y3-Y1 ) ( S ) ( Y - Y1 ) ! ! by computing the inverse of the coefficient matrix and multiplying ! it by the right hand side to get R and S. ! ! The values of dRdX, dRdY, dSdX and dSdY are easily from the formulas ! for R and S. ! det = ( xn(2) - xn(1) ) * ( yn(3) - yn(1) ) & - ( xn(3) - xn(1) ) * ( yn(2) - yn(1) ) r = ( ( yn(3) - yn(1) ) * ( x - xn(1) ) & + ( xn(1) - xn(3) ) * ( y - yn(1) ) ) / det drdx = ( yn(3) - yn(1) ) / det drdy = ( xn(1) - xn(3) ) / det s = ( ( yn(1) - yn(2) ) * ( x - xn(1) ) & + ( xn(2) - xn(1) ) * ( y - yn(1) ) ) / det dsdx = ( yn(1) - yn(2) ) / det dsdy = ( xn(2) - xn(1) ) / det ! ! The basis functions can now be evaluated in terms of the ! reference coordinates R and S. It's also easy to determine ! the values of the derivatives with respect to R and S. ! if ( inode == 1 ) then b = 2.0D+00 * ( 1.0D+00 - r - s ) * ( 0.5D+00 - r - s ) dbdr = - 3.0D+00 + 4.0D+00 * r + 4.0D+00 * s dbds = - 3.0D+00 + 4.0D+00 * r + 4.0D+00 * s else if ( inode == 2 ) then b = 2.0D+00 * r * ( r - 0.5D+00 ) dbdr = - 1.0D+00 + 4.0D+00 * r dbds = 0.0D+00 else if ( inode == 3 ) then b = 2.0D+00 * s * ( s - 0.5D+00 ) dbdr = 0.0D+00 dbds = - 1.0D+00 + 4.0D+00 * s else if ( inode == 4 ) then b = 4.0D+00 * r * ( 1.0D+00 - r - s ) dbdr = 4.0D+00 - 8.0D+00 * r - 4.0D+00 * s dbds = - 4.0D+00 * r else if ( inode == 5 ) then b = 4.0D+00 * r * s dbdr = 4.0D+00 * s dbds = 4.0D+00 * r else if ( inode == 6 ) then b = 4.0D+00 * s * ( 1.0D+00 - r - s ) dbdr = - 4.0D+00 * s dbds = 4.0D+00 - 4.0D+00 * r - 8.0D+00 * s else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'QBF - Fatal error!' write ( *, '(a,i8)' ) ' Request for local basis function INODE = ', inode stop end if ! ! We need to convert the derivative information from (R(X,Y),S(X,Y)) ! to (X,Y) using the chain rule. ! dbdx = dbdr * drdx + dbds * dsdx dbdy = dbdr * drdy + dbds * dsdy return end subroutine quad_a ( node_xy, element_node, element_num, node_num, & nnodes, wq, xq, yq ) !*****************************************************************************80 ! !! QUAD_A sets the quadrature rule for assembly. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the X and ! Y coordinates of nodes. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer NNODES, the number of nodes used to form ! one element. ! ! Output, real ( kind = rk ) WQ(3), quadrature weights. ! ! Output, real ( kind = rk ) XQ(3,ELEMENT_NUM), YQ(3,ELEMENT_NUM), the ! coordinates of the quadrature points in each element. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer node_num integer element integer, dimension(nnodes,element_num) :: element_node integer ip1 integer ip2 integer ip3 real ( kind = rk ), dimension(2,node_num) :: node_xy real ( kind = rk ), dimension(3) :: wq real ( kind = rk ) x1 real ( kind = rk ) x2 real ( kind = rk ) x3 real ( kind = rk ), dimension(3,element_num) :: xq real ( kind = rk ) y1 real ( kind = rk ) y2 real ( kind = rk ) y3 real ( kind = rk ), dimension(3,element_num) :: yq wq(1) = 1.0D+00 / 3.0D+00 wq(2) = wq(1) wq(3) = wq(1) do element = 1, element_num ip1 = element_node(1,element) ip2 = element_node(2,element) ip3 = element_node(3,element) x1 = node_xy(1,ip1) x2 = node_xy(1,ip2) x3 = node_xy(1,ip3) y1 = node_xy(2,ip1) y2 = node_xy(2,ip2) y3 = node_xy(2,ip3) xq(1,element) = 0.5D+00 * ( x1 + x2 ) xq(2,element) = 0.5D+00 * ( x2 + x3 ) xq(3,element) = 0.5D+00 * ( x1 + x3 ) yq(1,element) = 0.5D+00 * ( y1 + y2 ) yq(2,element) = 0.5D+00 * ( y2 + y3 ) yq(3,element) = 0.5D+00 * ( y1 + y3 ) end do return end subroutine quad_e ( node_xy, element_node, element, element_num, nnodes, & node_num, nqe, wqe, xqe, yqe ) !*****************************************************************************80 ! !! QUAD_E sets up the quadrature rule for error integration. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the X and Y ! coordinates of nodes. ! ! Input, integer ELEMENT_NODE(NNODES,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer ELEMENT, the index of the element for which ! the quadrature points are to be computed. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer NNODES, the number of nodes used to form ! one element. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer NQE, the number of points in the quadrature ! rule. This is actually fixed at 13. ! ! Output, real ( kind = rk ) WQE(13), the quadrature weights. ! ! Output, real ( kind = rk ) XQE(NQE), YQE(NQE), the X and Y coordinates ! of the quadrature points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer nnodes integer node_num integer nqe integer element integer, dimension(nnodes,element_num) :: element_node integer i integer ii integer iii integer ip1 integer ip2 integer ip3 real ( kind = rk ), dimension(2,node_num) :: node_xy real ( kind = rk ), dimension(nqe) :: wqe real ( kind = rk ), dimension(nqe) :: xqe real ( kind = rk ), dimension(nqe) :: yqe real ( kind = rk ) x1 real ( kind = rk ) x2 real ( kind = rk ) x3 real ( kind = rk ) y1 real ( kind = rk ) y2 real ( kind = rk ) y3 real ( kind = rk ) z1 real ( kind = rk ) z2 real ( kind = rk ) z3 real ( kind = rk ) z4 real ( kind = rk ) z5 real ( kind = rk ) z6 real ( kind = rk ) 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 = element_node(1,element) ip2 = element_node(2,element) ip3 = element_node(3,element) x1 = node_xy(1,ip1) x2 = node_xy(1,ip2) x3 = node_xy(1,ip3) y1 = node_xy(2,ip1) y2 = node_xy(2,ip2) y3 = node_xy(2,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 subroutine r8vec_print_some ( n, a, i_lo, i_hi, title ) !*****************************************************************************80 ! !! R8VEC_PRINT_SOME prints "some" of an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of entries of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, integer I_LO, I_HI, the first and last indices ! to print. The routine expects 1 <= I_LO <= I_HI <= N. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) integer i integer i_hi integer i_lo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = max ( i_lo, 1 ), min ( i_hi, n ) write ( *, '(2x,i8,2x,g14.8)' ) i, a(i) end do return end function rhs ( x, y ) !*****************************************************************************80 ! !! RHS gives the right-hand side of the differential equation. ! ! Discussion: ! ! The function specified here depends on the problem being ! solved. The user must be sure that RHS and EXACT are consistent. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) X, Y, the coordinates of a point ! in the region, at which the right hand side of the ! differential equation is to be evaluated. ! ! Output, real ( kind = rk ) RHS, the value of the right ! hand side of the differential equation at (X,Y). ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 real ( kind = rk ) rhs real ( kind = rk ) x real ( kind = rk ) y rhs = 2.0D+00 * pi * pi * sin ( pi * x ) * sin ( pi * y ) return end subroutine solution_write ( f, indx, node_num, nunk, output_filename, & node_xy ) !*****************************************************************************80 ! !! SOLUTION_WRITE writes the solution to a file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) F(NUNK), the coefficients of the solution. ! ! Input, integer INDX(NODE_NUM), gives the index of the unknown ! quantity associated with the given node. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer NUNK, the number of unknowns. ! ! Input, character ( len = * ) OUTPUT_FILENAME, the name of the file ! in which the data should be stored. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the X and Y ! coordinates of nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer nunk real ( kind = rk ) dudx real ( kind = rk ) dudy real ( kind = rk ), dimension(nunk) :: f integer, dimension(node_num) :: indx integer node real ( kind = rk ), dimension(2,node_num) :: node_xy character ( len = * ) :: output_filename integer output_status integer output_unit real ( kind = rk ) u real ( kind = rk ) x real ( kind = rk ) y call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, status = 'replace', & iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SOLUTION_WRITE - Warning!' write ( *, '(a)' ) ' Could not write the solution file.' return end if do node = 1, node_num x = node_xy(1,node) y = node_xy(2,node) if ( 0 < indx(node) ) then u = f(indx(node)) else call exact ( x, y, u, dudx, dudy ) end if write ( output_unit, '(g14.6)' ) u end do close ( unit = output_unit ) return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 August 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine triangulation_order6_plot ( file_name, node_num, node_xy, tri_num, & triangle_node, node_show, triangle_show ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER6_PLOT plots a 6-node triangulation of a pointset. ! ! Discussion: ! ! The triangulation is most usually a Delaunay triangulation, ! but this is not necessary. ! ! In a six node triangulation, it is assumed that nodes 1, 2, and 3 ! are the vertices of the triangles, and that nodes 4, 5, and 6 ! lie between 1 and 2, 2 and 3, and 3 and 1 respectively. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer NODE_NUM, the number of points. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the nodes. ! ! Input, integer TRI_NUM, the number of triangles. ! ! Input, integer TRIANGLE_NODE(6,TRI_NUM), lists, for ! each triangle, the indices of the points that form the vertices of ! the triangle. ! ! Input, integer NODE_SHOW, ! 0, do not show nodes; ! 1, show nodes; ! 2, show nodes and label them. ! ! Input, integer TRIANGLE_SHOW, ! 0, do not show triangles; ! 1, show triangles; ! 2, show triangles and label them. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer tri_num real ( kind = rk ) ave_x real ( kind = rk ) ave_y integer :: circle_size integer delta character ( len = * ) file_name integer file_unit integer i integer ios integer node integer node_show real ( kind = rk ) node_xy(2,node_num) character ( len = 40 ) string integer triangle integer triangle_node(6,tri_num) integer triangle_show real ( kind = rk ) x_max real ( kind = rk ) x_min integer x_ps integer :: x_ps_max = 576 integer :: x_ps_max_clip = 594 integer :: x_ps_min = 36 integer :: x_ps_min_clip = 18 real ( kind = rk ) x_scale real ( kind = rk ) y_max real ( kind = rk ) y_min integer y_ps integer :: y_ps_max = 666 integer :: y_ps_max_clip = 684 integer :: y_ps_min = 126 integer :: y_ps_min_clip = 108 real ( kind = rk ) y_scale ! ! We need to do some figuring here, so that we can determine ! the range of the data, and hence the height and width ! of the piece of paper. ! x_max = maxval ( node_xy(1,1:node_num) ) x_min = minval ( node_xy(1,1:node_num) ) x_scale = x_max - x_min x_max = x_max + 0.05D+00 * x_scale x_min = x_min - 0.05D+00 * x_scale x_scale = x_max - x_min y_max = maxval ( node_xy(2,1:node_num) ) y_min = minval ( node_xy(2,1:node_num) ) y_scale = y_max - y_min y_max = y_max + 0.05D+00 * y_scale y_min = y_min - 0.05D+00 * y_scale y_scale = y_max - y_min if ( x_scale < y_scale ) then delta = nint ( real ( x_ps_max - x_ps_min, kind = rk ) & * ( y_scale - x_scale ) / ( 2.0D+00 * y_scale ) ) x_ps_max = x_ps_max - delta x_ps_min = x_ps_min + delta x_ps_max_clip = x_ps_max_clip - delta x_ps_min_clip = x_ps_min_clip + delta x_scale = y_scale else if ( y_scale < x_scale ) then delta = nint ( real ( y_ps_max - y_ps_min, kind = rk ) & * ( x_scale - y_scale ) / ( 2.0D+00 * x_scale ) ) y_ps_max = y_ps_max - delta y_ps_min = y_ps_min + delta y_ps_max_clip = y_ps_max_clip - delta y_ps_min_clip = y_ps_min_clip + delta y_scale = x_scale end if ! ! Open the file. ! call get_unit ( file_unit ) open ( unit = file_unit, file = file_name, status = 'replace', & iostat = ios ) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGULATION_ORDER6_PLOT - Fatal error!' write ( *, '(a)' ) ' Can not open output file.' return end if ! ! Write the header. ! write ( file_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( file_unit, '(a)' ) '%%Creator: triangulation_order6_plot.f90' write ( file_unit, '(a)' ) '%%Title: ' // trim ( file_name ) write ( file_unit, '(a)' ) '%%Pages: 1' write ( file_unit, '(a,i3,2x,i3,2x,i3,2x,i3)' ) '%%BoundingBox: ', & x_ps_min, y_ps_min, x_ps_max, y_ps_max write ( file_unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( file_unit, '(a)' ) '%%LanguageLevel: 1' write ( file_unit, '(a)' ) '%%EndComments' write ( file_unit, '(a)' ) '%%BeginProlog' write ( file_unit, '(a)' ) '/inch {72 mul} def' write ( file_unit, '(a)' ) '%%EndProlog' write ( file_unit, '(a)' ) '%%Page: 1 1' write ( file_unit, '(a)' ) 'save' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB line color to very light gray.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.900 0.900 0.900 setrgbcolor' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Draw a gray border around the page.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) 'newpath' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_min, ' moveto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_max, y_ps_min, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_max, y_ps_max, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_max, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_min, ' lineto' write ( file_unit, '(a)' ) 'stroke' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB line color to black.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.000 0.000 0.000 setrgbcolor' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the font and its size.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '/Times-Roman findfont' write ( file_unit, '(a)' ) '0.50 inch scalefont' write ( file_unit, '(a)' ) 'setfont' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Print a title.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% 210 702 moveto' write ( file_unit, '(a)' ) '% (Triangulation) show' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Define a clipping polygon.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) 'newpath' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_min_clip, ' moveto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_max_clip, y_ps_min_clip, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_max_clip, y_ps_max_clip, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_max_clip, ' lineto' write ( file_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_min_clip, ' lineto' write ( file_unit, '(a)' ) 'clip newpath' ! ! Draw the nodes. ! if ( node_num <= 200 ) then circle_size = 5 else if ( node_num <= 500 ) then circle_size = 4 else if ( node_num <= 1000 ) then circle_size = 3 else if ( node_num <= 5000 ) then circle_size = 2 else circle_size = 1 end if if ( 1 <= node_show ) then write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Draw filled dots at the nodes.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB color to blue.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.000 0.150 0.750 setrgbcolor' write ( file_unit, '(a)' ) '%' do node = 1, node_num x_ps = int ( & ( ( x_max - node_xy(1,node) ) * real ( x_ps_min, kind = rk ) & + ( node_xy(1,node) - x_min ) * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( file_unit, '(a,i4,2x,i4,2x,i4,2x,a)' ) 'newpath ', x_ps, y_ps, & circle_size, '0 360 arc closepath fill' end do end if ! ! Label the nodes. ! if ( 2 <= node_show ) then write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Label the nodes:' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB color to darker blue.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.100 0.250 0.850 setrgbcolor' write ( file_unit, '(a)' ) '/Times-Roman findfont' write ( file_unit, '(a)' ) '0.20 inch scalefont' write ( file_unit, '(a)' ) 'setfont' do node = 1, node_num x_ps = int ( & ( ( x_max - node_xy(1,node) ) * real ( x_ps_min, kind = rk ) & + ( + node_xy(1,node) - x_min ) * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( string, '(i4)' ) node string = adjustl ( string ) write ( file_unit, '(i4,2x,i4,a)' ) x_ps, y_ps+5, & ' moveto (' // trim ( string ) // ') show' end do end if ! ! Draw the triangles. ! if ( 1 <= triangle_show ) then write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB color to red.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.900 0.200 0.100 setrgbcolor' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Draw the triangles.' write ( file_unit, '(a)' ) '%' do triangle = 1, tri_num write ( file_unit, '(a)' ) 'newpath' node = triangle_node(6,triangle) x_ps = int ( & ( ( x_max - node_xy(1,node) ) * real ( x_ps_min, kind = rk ) & + ( node_xy(1,node) - x_min ) * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( file_unit, '(i3,2x,i3,2x,a)' ) x_ps, y_ps, ' moveto' do i = 1, 3 node = triangle_node(i,triangle) x_ps = int ( & ( ( x_max - node_xy(1,node) ) & * real ( x_ps_min, kind = rk ) & + ( node_xy(1,node) - x_min ) & * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) & * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) & * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( file_unit, '(i3,2x,i3,2x,a)' ) x_ps, y_ps, ' lineto' node = triangle_node(i+3,triangle) x_ps = int ( & ( ( x_max - node_xy(1,node) ) & * real ( x_ps_min, kind = rk ) & + ( node_xy(1,node) - x_min ) & * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - node_xy(2,node) ) & * real ( y_ps_min, kind = rk ) & + ( node_xy(2,node) - y_min ) & * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( file_unit, '(i3,2x,i3,2x,a)' ) x_ps, y_ps, ' lineto' end do write ( file_unit, '(a)' ) 'stroke' end do end if ! ! Label the triangles. ! if ( 2 <= triangle_show ) then write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Label the triangles:' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% Set the RGB color to darker red.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '0.950 0.250 0.150 setrgbcolor' write ( file_unit, '(a)' ) '/Times-Roman findfont' write ( file_unit, '(a)' ) '0.20 inch scalefont' write ( file_unit, '(a)' ) 'setfont' do triangle = 1, tri_num ave_x = 0.0D+00 ave_y = 0.0D+00 do i = 1, 6 node = triangle_node(i,triangle) ave_x = ave_x + node_xy(1,node) ave_y = ave_y + node_xy(2,node) end do ave_x = ave_x / 6.0D+00 ave_y = ave_y / 6.0D+00 x_ps = int ( & ( ( x_max - ave_x ) * real ( x_ps_min, kind = rk ) & + ( + ave_x - x_min ) * real ( x_ps_max, kind = rk ) ) & / ( x_max - x_min ) ) y_ps = int ( & ( ( y_max - ave_y ) * real ( y_ps_min, kind = rk ) & + ( ave_y - y_min ) * real ( y_ps_max, kind = rk ) ) & / ( y_max - y_min ) ) write ( string, '(i4)' ) triangle string = adjustl ( string ) write ( file_unit, '(i4,2x,i4,a)' ) x_ps, y_ps, ' moveto (' & // trim ( string ) // ') show' end do end if write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) 'restore showpage' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '% End of page.' write ( file_unit, '(a)' ) '%' write ( file_unit, '(a)' ) '%%Trailer' write ( file_unit, '(a)' ) '%%EOF' close ( unit = file_unit ) return end subroutine xy_set ( nx, ny, node_num, xl, xr, yb, yt, node_xy ) !*****************************************************************************80 ! !! XY_SET sets the XY coordinates of the nodes. ! ! Discussion: ! ! The nodes are laid out in an evenly spaced grid, in the unit square. ! ! The first node is at the origin. More nodes are created to the ! right until the value of X = 1 is reached, at which point ! the next layer is generated starting back at X = 0, and an ! increased value of Y. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 September 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NX, NY, the number of elements in the X and ! Y direction. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) XL, XR, YB, YT, the X coordinates of ! the left and right sides of the rectangle, and the Y coordinates ! of the bottom and top of the rectangle. ! ! Output, real ( kind = rk ) NODE_XY(2,NODE_NUM), ! the coordinates of the nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer i integer j real ( kind = rk ) node_xy(2,node_num) integer nx integer ny real ( kind = rk ) xl real ( kind = rk ) xr real ( kind = rk ) yb real ( kind = rk ) yt do j = 1, 2 * ny - 1 do i = 1, 2 * nx - 1 node_xy(1,i+(j-1)*(2*nx-1)) = & ( real ( 2 * nx - i - 1, kind = rk ) * xl & + real ( i - 1, kind = rk ) * xr ) & / real ( 2 * nx - 2, kind = rk ) node_xy(2,i+(j-1)*(2*nx-1)) = & ( real ( 2 * ny - j - 1, kind = rk ) * yb & + real ( j - 1, kind = rk ) * yt ) & / real ( 2 * ny - 2, kind = rk ) end do end do return end