program main !*****************************************************************************80 ! !! fem2d_heat() solves the heat equation. ! ! Discussion: ! ! The heat equation being solved is: ! ! dUdT - Laplacian U(X,Y,T) + K(X,Y,T) * U(X,Y,T) = F(X,Y,T) ! ! in a triangulated region in the plane. ! ! Along the boundary of the region, Dirichlet conditions ! are imposed: ! ! U(X,Y,T) = G(X,Y,T) ! ! At the initial time T_INIT, the value of U is given ! at all points in the region: ! ! U(X,Y,T_INIT) = H(X,Y) ! ! The code uses continuous piecewise linear basis functions on ! triangles. ! ! The backward Euler approximation is used for the time derivatives. ! ! Problem specification: ! ! The user defines the geometry by supplying two data files ! which list the node coordinates, and list the nodes that make up ! each triangular element.. ! ! The user specifies the coefficient function K(X,Y,T) ! by supplying a routine of the form ! ! subroutine k_coef ( node_num, node_xy, time, k ) ! ! The user specifies the right hand side ! by supplying a routine of the form ! ! subroutine rhs ( node_num, node_xy, time, f ) ! ! The user specifies the right hand side of the Dirichlet boundary ! conditions by supplying a function ! ! subroutine dirichlet_condition ( node_num, node_xy, time, u ) ! ! The user specifies the initial condition by supplying a function ! ! subroutine initial_condition ( node_num, node_xy, time, u ) ! ! Usage: ! ! fem2d_heat node_file element_file ! ! invokes the program: ! ! * "node_file" contains the coordinates of the nodes; ! * "element_file" contains the indices of nodes that make up each ! triangular element. ! ! Files created include: ! ! * "nodes.eps", an image of the nodes; ! * "elements.eps", an image of the elements; ! * "u0000.txt", the value of the solution at the initial condition. ! * "u0001.txt" through "UNNNN.txt", the value of the solution at ! each time step; ! * "time.txt", the value of time at each step, from the initial to ! final times. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 25 December 2010 ! ! Author: ! ! John Burkardt ! ! Local: ! ! real ( kind = rk ) A(3*IB+1,NODE_NUM), the coefficient matrix. ! ! integer DIM_NUM, the spatial dimension, which is 2. ! ! character ( len = 255 ) ELEMENT_FILE_NAME, the name of the ! input file containing the element information. ! ! integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! integer ELEMENT_NUM, the number of elements. ! ! integer ELEMENT_ORDER, the order of each element. ! ! real ( kind = rk ) F(NODE_NUM), the right hand side. ! ! integer IB, the half-bandwidth of the matrix. ! ! logical NODE_BOUNDARY(NODE_NUM), is TRUE if a given ! node is on the boundary. ! ! integer NODE_CONDITION(NODE_NUM), indicates the type of ! boundary condition being applied to nodes on the boundary. ! 0, there is no condition (and no variable) at this node. ! 1, a finite element equation is used; ! 2, a Dirichlet condition is used. ! 3, a Neumann condition is used. ! ! character ( len = 255 ) NODE_FILE_NAME, the name of the ! input file containing the node coordinate information. ! ! integer NODE_NUM, the number of nodes. ! ! real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of nodes. ! ! integer QUAD_NUM, the number of quadrature points used ! for assembly. This is currently set to 3, the lowest reasonable value. ! Legal values are 1, 3, 4, 6, 7, 9, 13, and for some problems, a value of ! QUAD_NUM greater than 3 may be appropriate. ! ! real ( kind = rk ) U(NODE_NUM), the finite element coefficients ! defining the solution at the current time. ! ! real ( kind = rk ) U_OLD(NODE_NUM), the finite element coefficients ! defining the solution at the previous time. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ), allocatable, dimension (:,:) :: a logical, parameter :: debug = .false. integer dim_num character ( len = 255 ) element_file_name integer, allocatable, dimension(:,:) :: element_node integer element_num integer element_order real ( kind = rk ), allocatable, dimension (:) :: f integer ib integer ierr integer job integer node logical, allocatable, dimension(:) :: node_boundary integer, allocatable, dimension(:) :: node_condition character ( len = 255 ) node_file_name logical node_label integer node_num integer node_show real ( kind = rk ), allocatable, dimension(:,:) :: node_xy integer, allocatable, dimension (:) :: pivot integer, parameter :: quad_num = 7 character ( len = 255 ) solution_file_name real ( kind = rk ) time character ( len = 255 ) time_file_name real ( kind = rk ) time_final real ( kind = rk ) time_init real ( kind = rk ) time_old integer time_step integer time_step_num real ( kind = rk ) time_step_size integer time_unit integer triangle_show character ( len = 255 ) triangulation_eps_file_name real ( kind = rk ), allocatable, dimension ( : ) :: u real ( kind = rk ), allocatable, dimension ( : ) :: u_old call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_HEAT' write ( *, '(a)' ) ' FORTRAN90 version:' write ( *, '(a)' ) ' Solution of the time dependent heat equation' write ( *, '(a)' ) ' on an arbitrary triangulated region D in 2 dimensions.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Ut - Uxx - Uyy + K(x,y,t) * U = F(x,y,t) in D;' write ( *, '(a)' ) ' U = G(x,y,t) on boundary;' write ( *, '(a)' ) & ' U = H(x,y,t) at initial time.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The finite element method is used, with' write ( *, '(a)' ) ' 6 node quadratic triangular elements ("T6").' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The time derivative is approximated using the' write ( *, '(a)' ) ' backward Euler formula.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Current status:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' * Time step information currently set internally!' write ( *, '(a)' ) ' * Would be easy to do linear triangles as well.' write ( *, '(a)' ) ' * Do you want ability to compare to an exact solution?' write ( *, '(a)' ) ' ' ! ! Get the file names. ! call file_name_specification ( node_file_name, element_file_name ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Node file is "' // trim ( node_file_name ) // '".' write ( *, '(a)' ) ' Element file is "' // trim ( element_file_name ) & // '".' ! ! Read the node coordinate file. ! call r8mat_header_read ( node_file_name, dim_num, node_num ) write ( *, '(a,i8)' ) ' Number of nodes = ', node_num allocate ( node_boundary(node_num) ) allocate ( node_condition(node_num) ) allocate ( node_xy(dim_num,node_num) ) call r8mat_data_read ( node_file_name, dim_num, node_num, node_xy ) call r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, & dim_num, 10, ' First 10 nodes' ) ! ! Read the element description file. ! call i4mat_header_read ( element_file_name, element_order, element_num ) write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Element order = ', element_order write ( *, '(a,i8)' ) ' Number of elements = ', element_num if ( element_order /= 6 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_HEAT - Fatal error!' write ( *, '(a,i8)' ) ' The input element has order ', element_order write ( *, '(a)' ) ' However, elements of order 6 are required.' stop end if allocate ( element_node(element_order,element_num) ) call i4mat_data_read ( element_file_name, element_order, element_num, & element_node ) call i4mat_transpose_print_some ( element_order, element_num, & element_node, 1, 1, element_order, 10, ' First 10 elements' ) write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Quadrature order = ', quad_num ! ! Determine which nodes are boundary nodes and which have a ! finite element unknown. Then set the boundary values. ! call triangulation_order6_boundary_node ( node_num, element_num, & element_node, node_boundary ) if ( debug ) then call lvec_print ( node_num, node_boundary, ' Node Boundary?' ) end if ! ! Determine the node conditions. ! For now, we'll just assume all boundary nodes are Dirichlet. ! node_condition(1:node_num) = 1 do node = 1, node_num if ( node_boundary(node) ) then node_condition(node) = 2 end if end do ! ! Determine the bandwidth of the coefficient matrix. ! call bandwidth ( element_order, element_num, element_node, ib ) write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' The matrix half bandwidth is ', ib write ( *, '(a,i8)' ) ' The matrix bandwidth is ', 2 * ib + 1 write ( *, '(a,i8)' ) ' The storage bandwidth is ', 3 * ib + 1 ! ! Make a picture of the nodes. ! if ( node_num <= 100 ) then node_file_name = 'nodes.eps' node_label = .true. call points_plot ( node_file_name, node_num, node_xy, node_label ) end if ! ! Make a picture of the elements. ! if ( node_num <= 100 ) then triangulation_eps_file_name = '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 ) end if ! ! Set time stepping quantities. ! time_init = 0.0D+00 time_final = 0.5D+00 time_step_num = 10 time_step_size = ( time_final - time_init ) / real ( time_step_num, kind = rk ) write ( *, '(a)' ) ' ' write ( *, '(a,g14.6)' ) ' Initial time = ', time_init write ( *, '(a,g14.6)' ) ' Final time = ', time_final write ( *, '(a,g14.6)' ) ' Step size = ', time_step_size write ( *, '(a,i8)' ) ' Number of steps = ', time_step_num ! ! Allocate space for the coefficient matrix A and right hand side F. ! allocate ( a(3*ib+1,node_num) ) allocate ( f(node_num) ) allocate ( pivot(node_num) ) allocate ( u(node_num) ) allocate ( u_old(node_num) ) ! ! Initialize the names of the time and solution file. ! time_file_name = 'time.txt' solution_file_name = 'u0000.txt' ! ! Set the value of U at the initial time. ! time = time_init call initial_condition ( node_num, node_xy, time, u ) call get_unit ( time_unit ) open ( unit = time_unit, file = time_file_name, status = 'replace' ) write ( time_unit, '(g14.6)' ) time call solution_write ( node_num, u, solution_file_name, time ) ! ! Time looping. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'Time L2 Error H1 Error' write ( *, '(a)' ) ' ' do time_step = 1, time_step_num time_old = time u_old(1:node_num) = u(1:node_num) time = ( real ( time_step_num - time_step, kind = rk ) * time_init & + real ( time_step, kind = rk ) * time_final ) & / real ( time_step_num, kind = rk ) ! ! Assemble the finite element coefficient matrix A and the right-hand side F. ! call assemble_heat ( node_num, node_xy, node_condition, & element_order, element_num, element_node, quad_num, ib, time, a, f ) if ( debug ) then call dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, & 10, 10, ' Initial block of matrix A:' ) call r8vec_print_some ( node_num, f, 1, 10, & ' Part of right hand side F:' ) end if ! ! Adjust the linear system for the dU/dT term, which we are treating ! using the backward Euler formula. ! call assemble_backward_euler ( node_num, node_xy, element_order, & element_num, element_node, quad_num, ib, time, time_step_size, & u_old, a, f ) if ( debug ) then call dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, & 10, 10, ' A after DT adjustment:' ) call r8vec_print_some ( node_num, f, 1, 10, & ' F after DT adjustment:' ) end if ! ! Adjust the linear system to account for boundary conditions. ! call assemble_boundary ( node_num, node_xy, node_condition, ib, time, a, f ) if ( debug ) then call dgb_print_some ( node_num, node_num, ib, ib, a, 1, 1, & 10, 10, ' A after BC adjustment:' ) call r8vec_print_some ( node_num, f, 1, 10, & ' F after BC adjustment:' ) end if ! ! Solve the linear system using a banded solver. ! call dgb_fa ( node_num, ib, ib, a, pivot, ierr ) if ( ierr /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_HEAT - 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 u(1:node_num) = f(1:node_num) call dgb_sl ( node_num, ib, ib, a, pivot, u, job ) if ( debug ) then call r8vec_print_some ( node_num, u, 1, 10, & ' Part of the solution vector U:' ) end if ! ! Increment the file name, and write the new solution. ! write ( time_unit, '(g14.6)' ) time call file_name_inc ( solution_file_name ) call solution_write ( node_num, u, solution_file_name, time ) end do close ( unit = time_unit ) ! ! Free memory. ! deallocate ( a ) deallocate ( f ) deallocate ( node_boundary ) deallocate ( node_condition ) deallocate ( node_xy ) deallocate ( pivot ) deallocate ( element_node ) deallocate ( u ) deallocate ( u_old ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FEM2D_HEAT:' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop end subroutine assemble_backward_euler ( node_num, node_xy, element_order, & element_num, element_node, quad_num, ib, time, time_step_size, & u_old, a, f ) !*****************************************************************************80 ! !! ASSEMBLE_BACKWARD_EULER adjusts the system for the backward Euler term. ! ! Discussion: ! ! The input linear system ! ! A * U = F ! ! is appropriate for the equation ! ! -Uxx - Uyy - K * U = RHS ! ! We need to modify the matrix A and the right hand side F to ! account for the approximation of the time derivative in ! ! Ut - Uxx - Uyy - K * U = RHS ! ! by the backward Euler approximation: ! ! Ut approximately equal to ( U - Uold ) / dT ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 July 2007 ! ! 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 ELEMENT_ORDER, the number of nodes used ! to form one element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer QUAD_NUM, the number of quadrature points ! used in assembly. ! ! Input, integer IB, the half-bandwidth of the matrix. ! ! Input, real ( kind = rk ) TIME, the current time. ! ! Input, real ( kind = rk ) TIME_STEP_SIZE, the size of the time step. ! ! Input, real ( kind = rk ) U_OLD(NODE_NUM), the finite element ! coefficients for the solution at the previous time. ! ! Input/output, real ( kind = rk ) A(3*IB+1,NODE_NUM), the NODE_NUM ! by NODE_NUM coefficient matrix, stored in a compressed format. ! ! Input/output, real ( kind = rk ) F(NODE_NUM), the right hand side. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ib integer element_num integer element_order integer node_num integer quad_num real ( kind = rk ), dimension(3*ib+1,node_num) :: a real ( kind = rk ) area 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 integer, dimension(element_order,element_num) :: element_node real ( kind = rk ), dimension(node_num) :: f integer j integer node real ( kind = rk ), dimension(2,node_num) :: node_xy real ( kind = rk ), dimension(2,quad_num) :: phys_xy integer quad real ( kind = rk ), dimension (quad_num) :: quad_w real ( kind = rk ), dimension (2,quad_num) :: quad_xy real ( kind = rk ) t3(2,3) real ( kind = rk ) t6(2,6) integer test real ( kind = rk ) time real ( kind = rk ) time_step_size real ( kind = rk ) triangle_area_2d real ( kind = rk ), dimension(node_num) :: u_old real ( kind = rk ), dimension(quad_num) :: w call r8_fake_use ( time ) ! ! Get the quadrature rule weights and nodes. ! call quad_rule ( quad_num, quad_w, quad_xy ) do element = 1, element_num ! ! Make two copies of the triangle. ! t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element)) t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element)) ! ! Map the quadrature points QUAD_XY to points PHYS_XY in the physical triangle. ! call reference_to_physical_t3 ( t3, quad_num, quad_xy, phys_xy ) area = abs ( triangle_area_2d ( t3 ) ) w(1:quad_num) = area * quad_w(1:quad_num) do quad = 1, quad_num do test = 1, element_order node = element_node(test,element) call basis_11_t6 ( t6, test, phys_xy(1:2,quad), bi, dbidx, dbidy ) ! ! Carry the U_OLD term to the right hand side. ! f(node) = f(node) + w(quad) * bi * u_old(node) / time_step_size ! ! Modify the diagonal entries of A. ! do basis = 1, element_order j = element_node(basis,element) call basis_11_t6 ( t6, basis, phys_xy(1:2,quad), bj, dbjdx, dbjdy ) a(node-j+2*ib+1,j) = a(node-j+2*ib+1,j) & + w(quad) * bi * bj / time_step_size end do end do end do end do return end subroutine assemble_boundary ( node_num, node_xy, node_condition, ib, time, & a, f ) !*****************************************************************************80 ! !! ASSEMBLE_BOUNDARY modifies the linear system for the boundary conditions. ! ! Discussion: ! ! For now, we are only working with Dirichlet boundary conditions. ! ! Neumann conditions will be OK if the specified flux value is zero. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 January 2007 ! ! 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 NODE_CONDITION(NODE_NUM), reports the condition ! used to set the unknown associated with the node. ! 0, unknown. ! 1, finite element equation. ! 2, Dirichlet condition; ! 3, Neumann condition. ! ! Input, integer IB, the half-bandwidth of the matrix. ! ! Input, real ( kind = rk ) TIME, the current time. ! ! Input/output, real ( kind = rk ) A(3*IB+1,NODE_NUM), the NODE_NUM by ! NODE_NUM coefficient matrix, stored in a compressed format; on output, ! the matrix has been adjusted for Dirichlet boundary conditions. ! ! Input/output, real ( kind = rk ) F(NODE_NUM), the right hand side. ! On output, the right hand side has been adjusted for Dirichlet ! boundary conditions. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ib integer node_num real ( kind = rk ), dimension(3*ib+1,node_num) :: a real ( kind = rk ) bc_value(node_num) integer column integer column_high integer column_low integer, parameter :: DIRICHLET = 2 real ( kind = rk ), dimension(node_num) :: f integer node integer node_condition(node_num) real ( kind = rk ), dimension(2,node_num) :: node_xy real ( kind = rk ) time call dirichlet_condition ( node_num, node_xy, time, bc_value ) do node = 1, node_num if ( node_condition(node) == DIRICHLET ) then column_low = max ( node - ib, 1 ) column_high = min ( node + ib, node_num ) do column = column_low, column_high a(node-column+2*ib+1,column) = 0.0D+00 end do a(2*ib+1,node) = 1.0D+00 f(node) = bc_value(node) end if end do return end subroutine assemble_heat ( node_num, node_xy, node_condition, & element_order, element_num, element_node, quad_num, ib, time, a, f ) !*****************************************************************************80 ! !! ASSEMBLE_HEAT assembles the finite element system for the heat equation. ! ! 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: ! ! 30 August 2006 ! ! 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 NODE_CONDITION(NODE_NUM), reports the condition ! used to set the unknown associated with the node. ! 0, unknown. ! 1, finite element equation. ! 2, Dirichlet condition; ! 3, Neumann condition. ! ! Input, integer ELEMENT_ORDER, the number of nodes used to ! form one element. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer QUAD_NUM, the number of quadrature points ! used in assembly. ! ! Input, integer IB, the half-bandwidth of the matrix. ! ! Input, real ( kind = rk ) TIME, the current time. ! ! Output, real ( kind = rk ) A(3*IB+1,NODE_NUM), the NODE_NUM by NODE_NUM ! coefficient matrix, stored in a compressed format. ! ! Output, real ( kind = rk ) F(NODE_NUM), the right hand side. ! ! Local parameters: ! ! Local, real ( kind = rk ) BI, DBIDX, DBIDY, the value of some basis function ! and its first derivatives at a quadrature point. ! ! Local, real ( kind = rk ) BJ, DBJDX, DBJDY, the value of another basis ! function and its first derivatives at a quadrature point. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer element_order integer ib integer node_num integer quad_num real ( kind = rk ), dimension(3*ib+1,node_num) :: a real ( kind = rk ) area 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, dimension(6,element_num) :: element_node real ( kind = rk ), dimension(node_num) :: f integer i integer j real ( kind = rk ) k_value integer node_condition(node_num) real ( kind = rk ), dimension(2,node_num) :: node_xy real ( kind = rk ) phys_xy(2,quad_num) integer quad real ( kind = rk ), dimension(quad_num) :: quad_w real ( kind = rk ), dimension(2,quad_num) :: quad_xy real ( kind = rk ) :: rhs_value real ( kind = rk ), dimension (2,3) :: t3 real ( kind = rk ), dimension (2,6) :: t6 integer test real ( kind = rk ) time integer element real ( kind = rk ) triangle_area_2d real ( kind = rk ) w(quad_num) call i4_fake_use ( node_condition(1) ) ! ! Initialize the arrays to zero. ! f(1:node_num) = 0.0D+00 a(1:3*ib+1,1:node_num) = 0.0D+00 ! ! Get the quadrature weights and nodes. ! call quad_rule ( quad_num, quad_w, quad_xy ) ! ! Add up all quantities associated with the ELEMENT-th element. ! do element = 1, element_num ! ! Make two copies of the triangle. ! t3(1:2,1:3) = node_xy(1:2,element_node(1:3,element)) t6(1:2,1:6) = node_xy(1:2,element_node(1:6,element)) ! ! Map the quadrature points QUAD_XY to points PHYS_XY in the physical triangle. ! call reference_to_physical_t3 ( t3, quad_num, quad_xy, phys_xy ) area = abs ( triangle_area_2d ( t3 ) ) w(1:quad_num) = area * quad_w(1:quad_num) ! ! Consider the QUAD-th quadrature point. ! do quad = 1, quad_num call k_coef ( 1, phys_xy(1:2,quad), time, k_value ) call rhs ( 1, phys_xy(1:2,quad), time, rhs_value ) ! ! Consider the TEST-th test function. ! ! We generate an integral for every node associated with an unknown. ! But if a node is associated with a boundary condition, we do nothing. ! do test = 1, element_order i = element_node(test,element) call basis_11_t6 ( t6, test, phys_xy(1:2,quad), bi, dbidx, dbidy ) f(i) = f(i) + w(quad) * rhs_value * bi ! ! Consider the BASIS-th basis function, which is used to form the ! value of the solution function. ! do basis = 1, element_order j = element_node(basis,element) call basis_11_t6 ( t6, basis, phys_xy(1:2,quad), bj, dbjdx, dbjdy ) a(i-j+2*ib+1,j) = a(i-j+2*ib+1,j) & + w(quad) * ( dbidx * dbjdx + dbidy * dbjdy + k_value * bj * bi ) end do end do end do end do return end subroutine bandwidth ( element_order, element_num, element_node, nhba ) !*****************************************************************************80 ! !! BANDWIDTH determines the bandwidth of the coefficient matrix. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 August 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ELEMENT_ORDER, the order of the elements. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Output, integer NHBA, the half bandwidth of the matrix. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer element_order integer element integer element_node(element_order,element_num) integer global_i integer global_j integer local_i integer local_j integer nhba nhba = 0 do element = 1, element_num do local_i = 1, element_order global_i = element_node(local_i,element) do local_j = 1, element_order global_j = element_node(local_j,element) nhba = max ( nhba, abs ( global_j - global_i ) ) end do end do end do return end subroutine basis_11_t6 ( t, i, p, bi, dbidx, dbidy ) !*****************************************************************************80 ! !! BASIS_11_T6: one basis at one point for the T6 6 node triangular element. ! ! Discussion: ! ! The routine is given the coordinates of the nodes of a triangle. ! ! 3 ! / \ ! 6 5 ! / \ ! 1---4---2 ! ! It evaluates the quadratic basis function B(I)(X,Y) associated with ! node I, which has the property that it is a quadratic function ! which is 1 at node I and zero at the other five nodes. ! ! This routine assumes that the sides of the triangle are straight, ! so that the midside nodes fall on the line between two vertices. ! ! This routine relies on the fact that each basis function can be ! written as the product of two linear factors, which are easily ! computed and normalized. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) T(2,6), the coordinates of the nodes. ! ! Input, integer I, the index of the desired basis function. ! I should be between 1 and 6. ! ! Input, real ( kind = rk ) P(2), the coordinates of a point at which the ! basis function is to be evaluated. ! ! Output, real ( kind = rk ) BI, DBIDX, DBIDY, the values of the basis ! function and its X and Y derivatives. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) bi real ( kind = rk ) dbidx real ( kind = rk ) dbidy real ( kind = rk ) gf real ( kind = rk ) gn real ( kind = rk ) hf real ( kind = rk ) hn integer i integer i4_wrap integer j1 integer j2 integer k1 integer k2 real ( kind = rk ) p(2) real ( kind = rk ) t(2,6) if ( i < 1 .or. 6 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_11_T6 - Fatal error!' write ( *, '(a)' ) ' Basis index I is not between 1 and 6.' write ( *, '(a,i8)' ) ' I = ', i stop end if ! ! Determine the pairs of nodes. ! if ( i <= 3 ) then j1 = i4_wrap ( i + 1, 1, 3 ) j2 = i4_wrap ( i + 2, 1, 3 ) k1 = i + 3 k2 = i4_wrap ( i + 5, 4, 6 ) else j1 = i - 3 j2 = i4_wrap ( i - 3 + 2, 1, 3 ) k1 = i4_wrap ( i - 3 + 1, 1, 3 ) k2 = i4_wrap ( i - 3 + 2, 1, 3 ) end if ! ! Evaluate the two linear factors GF and HF, ! and their normalizers GN and HN. ! gf = ( p(1) - t(1,j1) ) * ( t(2,j2) - t(2,j1) ) & - ( t(1,j2) - t(1,j1) ) * ( p(2) - t(2,j1) ) gn = ( t(1,i) - t(1,j1) ) * ( t(2,j2) - t(2,j1) ) & - ( t(1,j2) - t(1,j1) ) * ( t(2,i) - t(2,j1) ) hf = ( p(1) - t(1,k1) ) * ( t(2,k2) - t(2,k1) ) & - ( t(1,k2) - t(1,k1) ) * ( p(2) - t(2,k1) ) hn = ( t(1,i) - t(1,k1) ) * ( t(2,k2) - t(2,k1) ) & - ( t(1,k2) - t(1,k1) ) * ( t(2,i) - t(2,k1) ) ! ! Construct the basis function and its derivatives. ! bi = ( gf / gn ) * ( hf / hn ) dbidx = ( ( t(2,j2) - t(2,j1) ) / gn ) * ( hf / hn ) & + ( gf / gn ) * ( ( t(2,k2) - t(2,k1) ) / hn ) dbidy = - ( ( t(1,j2) - t(1,j1) ) / gn ) * ( hf / hn ) & - ( gf / gn ) * ( ( t(1,k2) - t(1,k1) ) / hn ) return end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) !*****************************************************************************80 ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! Example: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C1, C2, the characters to compare. ! ! Output, logical CH_EQI, the result of the comparison. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) logical ch_eqi character c1 character c1_cap character c2 character c2_cap c1_cap = c1 c2_cap = c2 call ch_cap ( c1_cap ) call ch_cap ( c2_cap ) if ( c1_cap == c2_cap ) then ch_eqi = .true. else ch_eqi = .false. end if return end subroutine ch_to_digit ( c, digit ) !*****************************************************************************80 ! !! CH_TO_DIGIT returns the integer value of a base 10 digit. ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer DIGIT, the corresponding integer value. ! If C was 'illegal', then DIGIT is -1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer digit if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c == ' ' ) then digit = 0 else digit = -1 end if return end subroutine 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: ! ! 04 March 1999 ! ! Author: ! ! FORTRAN77 original version by Dongarra, Bunch, Moler, Stewart. ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Dongarra, Bunch, Moler, 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 ( l /= k ) 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_mxv ( m, n, ml, mu, a, x, b ) !*****************************************************************************80 ! !! DGB_MXV multiplies a DGB matrix times a vector. ! ! 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. ! ! LINPACK and LAPACK storage of general band matrices requires ! an extra ML upper diagonals for possible fill in entries during ! Gauss elimination. This routine does not access any entries ! in the fill in diagonals, because it assumes that the matrix ! has NOT had Gauss elimination applied to it. If the matrix ! has been Gauss eliminated, then the routine DGB_MU must be ! used instead. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 January 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Dongarra, Bunch, Moler, Stewart, ! LINPACK User's Guide, ! SIAM, 1979. ! ! 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, real ( kind = rk ) X(N), the vector to be multiplied by A. ! ! Output, real ( kind = rk ) B(M), the product A * x. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer ml integer mu integer m integer n real ( kind = rk ) a(2*ml+mu+1,n) real ( kind = rk ) b(m) integer i integer j integer jhi integer jlo real ( kind = rk ) x(n) do i = 1, m b(i) = 0.0D+00 jlo = max ( 1, i - ml ) jhi = min ( n, i + mu ) do j = jlo, jhi b(i) = b(i) + a(i-j+ml+mu+1,j) * x(j) end do end do 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. ! ! DGB storage is used by LINPACK and LAPACK. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 April 2006 ! ! 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 write ( *, '(a)' ) ' ' 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: ! ! 04 March 1999 ! ! Author: ! ! FORTRAN77 original version by Dongarra, Bunch, Moler, Stewart. ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! Dongarra, Bunch, Moler, 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 file_column_count ( input_file_name, column_num ) !*****************************************************************************80 ! !! FILE_COLUMN_COUNT counts the number of columns in the first line of a file. ! ! Discussion: ! ! The file is assumed to be a simple text file. ! ! Most lines of the file is presumed to consist of COLUMN_NUM words, ! separated by spaces. There may also be some blank lines, and some ! comment lines, ! which have a "#" in column 1. ! ! The routine tries to find the first non-comment non-blank line and ! counts the number of words in that line. ! ! If all lines are blanks or comments, it goes back and tries to analyze ! a comment line. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILE_NAME, the name of the file. ! ! Output, integer COLUMN_NUM, the number of columns ! in the file. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer column_num logical got_one character ( len = * ) input_file_name integer input_status integer input_unit character ( len = 256 ) line ! ! Open the file. ! call get_unit ( input_unit ) open ( unit = input_unit, file = input_file_name, status = 'old', & form = 'formatted', access = 'sequential', iostat = input_status ) if ( input_status /= 0 ) then column_num = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' & // trim ( input_file_name ) // '" on unit ', input_unit return end if ! ! Read one line, but skip blank lines and comment lines. ! got_one = .false. do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if if ( line(1:1) == '#' ) then cycle end if got_one = .true. exit end do if ( .not. got_one ) then rewind ( input_unit ) do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if got_one = .true. exit end do end if close ( unit = input_unit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!' write ( *, '(a)' ) ' The file does not seem to contain any data.' column_num = -1 return end if call s_word_count ( line, column_num ) return end subroutine file_name_inc ( file_name ) !*****************************************************************************80 ! !! FILE_NAME_INC increments a partially numeric filename. ! ! Discussion: ! ! It is assumed that the digits in the name, whether scattered or ! connected, represent a number that is to be increased by 1 on ! each call. If this number is all 9's on input, the output number ! is all 0's. Non-numeric letters of the name are unaffected. ! ! If the name is empty, then the routine stops. ! ! If the name contains no digits, the empty string is returned. ! ! Example: ! ! Input Output ! ----- ------ ! 'a7to11.txt' 'a7to12.txt' ! 'a7to99.txt' 'a8to00.txt' ! 'a9to99.txt' 'a0to00.txt' ! 'cat.txt' ' ' ! ' ' STOP! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) FILE_NAME. ! On input, a character string to be incremented. ! On output, the incremented string. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer change integer digit character ( len = * ) file_name integer i integer lens lens = len_trim ( file_name ) if ( lens <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_NAME_INC - Fatal error!' write ( *, '(a)' ) ' The input string is empty.' stop end if change = 0 do i = lens, 1, -1 c = file_name(i:i) if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then change = change + 1 digit = ichar ( c ) - 48 digit = digit + 1 if ( digit == 10 ) then digit = 0 end if c = char ( digit + 48 ) file_name(i:i) = c if ( c /= '0' ) then return end if end if end do if ( change == 0 ) then file_name = ' ' return end if return end subroutine file_name_specification ( node_file_name, element_file_name ) !*****************************************************************************80 ! !! FILE_NAME_SPECIFICATION determines the names of the input files. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 October 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, character ( len = * ) NODE_FILE_NAME, the name of the node file. ! ! Output, character ( len = * ) ELEMENT_FILE_NAME, the name ! of the element file. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer arg_num character ( len = * ) :: element_file_name integer iarg integer iargc character ( len = * ) :: node_file_name ! ! Get the number of command line arguments. ! arg_num = iargc ( ) ! ! If at least one command line argument, it's the node file name. ! if ( 1 <= arg_num ) then iarg = 1 call getarg ( iarg, node_file_name ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_NAME_SPECIFICATION:' write ( *, '(a)' ) ' Please enter the name of the node file.' read ( *, '(a)' ) node_file_name end if ! ! If at least two command line arguments, the second is the triangulation file. ! if ( 2 <= arg_num ) then iarg = 2 call getarg ( iarg, element_file_name ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_NAME_SPECIFICATION:' write ( *, '(a)' ) ' Please enter the name of the triangulation file.' read ( *, '(a)' ) element_file_name end if return end subroutine file_row_count ( input_file_name, row_num ) !*****************************************************************************80 ! !! FILE_ROW_COUNT counts the number of row records in a file. ! ! Discussion: ! ! It does not count lines that are blank, or that begin with a ! comment symbol '#'. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILE_NAME, the name of the input file. ! ! Output, integer ROW_NUM, the number of rows found. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer bad_num integer comment_num integer ierror character ( len = * ) input_file_name integer input_status integer input_unit character ( len = 100 ) line integer record_num integer row_num call get_unit ( input_unit ) open ( unit = input_unit, file = input_file_name, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then row_num = -1; ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_ROW_COUNT - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_file_name ) // '" on unit ', input_unit stop end if comment_num = 0 row_num = 0 record_num = 0 bad_num = 0 do read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = record_num exit end if record_num = record_num + 1 if ( line(1:1) == '#' ) then comment_num = comment_num + 1 cycle end if if ( len_trim ( line ) == 0 ) then comment_num = comment_num + 1 cycle end if row_num = row_num + 1 end do close ( unit = input_unit ) return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer 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 an integer 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: ! ! 18 September 2005 ! ! 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 i4_fake_use ( n ) !*****************************************************************************80 ! !! i4_fake_use() pretends to use a variable. ! ! Discussion: ! ! Some compilers will issue a warning if a variable is unused. ! Sometimes there's a good reason to include a variable in a program, ! but not to use it. Calling this function with that variable as ! the argument will shut the compiler up. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 April 2020 ! ! Author: ! ! John Burkardt ! ! Input: ! ! integer N, the variable to be "used". ! implicit none integer n if ( n /= n ) then write ( *, '(a)' ) ' i4_fake_use: variable is NAN.' end if return end function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of integer division. ! ! Formula: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! Discussion: ! ! The MOD function computes a result with the same sign as the ! quantity being divided. Thus, suppose you had an angle A, ! and you wanted to ensure that it was between 0 and 360. ! Then mod(A,360) would do, if A was positive, but if A ! was negative, your result would be between -360 and 0. ! ! On the other hand, I4_MODP(A,360) is between 0 and 360, always. ! ! Example: ! ! I J MOD I4_MODP Factorization ! ! 107 50 7 7 107 = 2 * 50 + 7 ! 107 -50 7 7 107 = -2 * -50 + 7 ! -107 50 -7 43 -107 = -3 * 50 + 43 ! -107 -50 -7 43 -107 = 3 * -50 + 43 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number to be divided. ! ! Input, integer J, the number that divides I. ! ! Output, integer I4_MODP, the nonnegative remainder when I is ! divided by J. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer i4_modp integer j if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i8)' ) ' I4_MODP ( I, J ) called with J = ', j stop end if i4_modp = mod ( i, j ) if ( i4_modp < 0 ) then i4_modp = i4_modp + abs ( j ) end if return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an integer to lie between given limits by wrapping. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I I4_WRAP ! ! -2 8 ! -1 4 ! 0 5 ! 1 6 ! 2 7 ! 3 8 ! 4 4 ! 5 5 ! 6 6 ! 7 7 ! 8 8 ! 9 4 ! 10 5 ! 11 6 ! 12 7 ! 13 8 ! 14 4 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 August 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer IVAL, a value. ! ! Input, integer ILO, IHI, the desired bounds for the value. ! ! Output, integer I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i4_modp integer i4_wrap integer ihi integer ilo integer ival integer jhi integer jlo integer wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then i4_wrap = jlo else i4_wrap = jlo + i4_modp ( ival - jlo, wide ) end if return end subroutine i4col_compare ( m, n, a, i, j, isgn ) !*****************************************************************************80 ! !! I4COL_COMPARE compares columns I and J of a integer array. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! ISGN = -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an array of N columns of vectors ! of length M. ! ! Input, integer I, J, the columns to be compared. ! I and J must be between 1 and N. ! ! Output, integer ISGN, the results of the comparison: ! -1, column I < column J, ! 0, column I = column J, ! +1, column J < column I. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer a(m,n) integer i integer isgn integer j integer k ! ! Check. ! if ( i < 1 .or. n < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index I is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i stop end if if ( j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_COMPARE - Fatal error!' write ( *, '(a)' ) ' Column index J is out of bounds.' write ( *, '(a,i8)' ) ' J = ', j stop end if isgn = 0 if ( i == j ) then return end if k = 1 do while ( k <= m ) if ( a(k,i) < a(k,j) ) then isgn = -1 return else if ( a(k,j) < a(k,i) ) then isgn = +1 return end if k = k + 1 end do return end subroutine i4col_sort_a ( m, n, a ) !*****************************************************************************80 ! !! I4COL_SORT_A ascending sorts an integer array of columns. ! ! Discussion: ! ! In lexicographic order, the statement "X < Y", applied to two real ! vectors X and Y of length M, means that there is some index I, with ! 1 <= I <= M, with the property that ! ! X(J) = Y(J) for J < I, ! and ! X(I) < Y(I). ! ! In other words, the first time they differ, X is smaller. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 25 September 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of rows of A, and the length of ! a vector of data. ! ! Input, integer N, the number of columns of A. ! ! Input/output, integer A(M,N). ! On input, the array of N columns of M-vectors. ! On output, the columns of A have been sorted in ascending ! lexicographic order. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer a(m,n) integer i integer indx integer isgn integer j if ( m <= 0 ) then return end if if ( n <= 1 ) then return end if ! ! Initialize. ! i = 0 indx = 0 isgn = 0 j = 0 ! ! Call the external heap sorter. ! do call sort_heap_external ( n, indx, i, j, isgn ) ! ! Interchange the I and J objects. ! if ( 0 < indx ) then call i4col_swap ( m, n, a, i, j ) ! ! Compare the I and J objects. ! else if ( indx < 0 ) then call i4col_compare ( m, n, a, i, j, isgn ) else if ( indx == 0 ) then exit end if end do return end subroutine i4col_swap ( m, n, a, i, j ) !*****************************************************************************80 ! !! I4COL_SWAP swaps columns I and J of a integer array of column data. ! ! Example: ! ! Input: ! ! M = 3, N = 4, I = 2, J = 4 ! ! A = ( ! 1 2 3 4 ! 5 6 7 8 ! 9 10 11 12 ) ! ! Output: ! ! A = ( ! 1 4 3 2 ! 5 8 7 6 ! 9 12 11 10 ) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 April 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns ! in the array. ! ! Input/output, integer A(M,N), an array of N columns ! of length M. ! ! Input, integer I, J, the columns to be swapped. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer a(m,n) integer col(m) integer i integer j if ( i < 1 .or. n < i .or. j < 1 .or. n < j ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4COL_SWAP - Fatal error!' write ( *, '(a)' ) ' I or J is out of bounds.' write ( *, '(a,i8)' ) ' I = ', i write ( *, '(a,i8)' ) ' J = ', j write ( *, '(a,i8)' ) ' N = ', n stop end if if ( i == j ) then return end if col(1:m) = a(1:m,i) a(1:m,i) = a(1:m,j) a(1:m,j) = col(1:m) return end subroutine i4mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! I4MAT_DATA_READ reads data from an I4MAT file. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! The file may contain more than N points, but this routine ! will return after reading N points. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 January 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Output, integer TABLE(M,N), the data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer ierror character ( len = * ) input_filename integer input_status integer input_unit integer j character ( len = 255 ) line integer table(m,n) integer x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_i4vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine i4mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! I4MAT_HEADER_READ reads the header from an I4MAT. ! ! Discussion: ! ! An I4MAT is an array of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer M, spatial dimension. ! ! Output, integer N, the number of points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = * ) input_filename integer m integer n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 7 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine lvec_print ( n, a, title ) !*****************************************************************************80 ! !! LVEC_PRINT prints a logical vector. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, logical A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n logical a(n) integer i character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,l1)' ) i, a(i) end do return end subroutine points_plot ( file_name, node_num, node_xy, node_label ) !*****************************************************************************80 ! !! POINTS_PLOT plots a pointset. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 September 2006 ! ! 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_status integer file_unit 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 = file_status ) if ( file_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POINTS_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: points_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 quad_rule ( quad_num, quad_w, quad_xy ) !*****************************************************************************80 ! !! QUAD_RULE sets the quadrature rule for assembly. ! ! Discussion: ! ! The quadrature rule is given for a reference element. ! ! 0 <= X, ! 0 <= Y, and ! X + Y <= 1. ! ! ^ ! 1 | * ! | |\ ! Y | | \ ! | | \ ! 0 | *---* ! +-------> ! 0 X 1 ! ! The rules have the following precision: ! ! QUAD_NUM Precision ! ! 1 1 ! 3 2 ! 4 3 ! 6 4 ! 7 5 ! 9 6 ! 13 7 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 July 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer QUAD_NUM, the number of quadrature nodes. ! ! Output, real ( kind = rk ) QUAD_W(QUAD_NUM), the quadrature weights. ! ! Output, real ( kind = rk ) QUAD_XY(2,QUAD_NUM), ! the coordinates of the quadrature nodes. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 integer quad_num real ( kind = rk ) a real ( kind = rk ) b real ( kind = rk ) c real ( kind = rk ) d real ( kind = rk ) e real ( kind = rk ) f real ( kind = rk ) g real ( kind = rk ) h real ( kind = rk ), dimension(quad_num) :: quad_w real ( kind = rk ), dimension(dim_num,quad_num) :: quad_xy real ( kind = rk ) t real ( kind = rk ) u real ( kind = rk ) v real ( kind = rk ) w if ( quad_num == 1 ) then quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & 1.0D+00 / 3.0D+00, 1.0D+00 / 3.0D+00 /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = 1.0D+00 else if ( quad_num == 3 ) then quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & 0.5D+00, 0.0D+00, & 0.5D+00, 0.5D+00, & 0.0D+00, 0.5D+00 /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = 1.0D+00 / 3.0D+00 else if ( quad_num == 4 ) then a = 6.0D+00 / 30.0D+00 b = 10.0D+00 / 30.0D+00 c = 18.0D+00 / 30.0D+00 d = 25.0D+00 / 48.0D+00 e = -27.0D+00 / 48.0D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & b, b, & c, a, & a, c, & a, a /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ e, d, d, d /) else if ( quad_num == 6 ) then a = 0.816847572980459D+00 b = 0.091576213509771D+00 c = 0.108103018168070D+00 d = 0.445948490915965D+00 v = 0.109951743655322D+00 w = 0.223381589678011D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & a, b, & b, a, & b, b, & c, d, & d, c, & d, d /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ v, v, v, w, w, w /) else if ( quad_num == 7 ) then a = 1.0D+00 / 3.0D+00 b = ( 9.0D+00 + 2.0D+00 * sqrt ( 15.0D+00 ) ) / 21.0D+00 c = ( 6.0D+00 - sqrt ( 15.0D+00 ) ) / 21.0D+00 d = ( 9.0D+00 - 2.0D+00 * sqrt ( 15.0D+00 ) ) / 21.0D+00 e = ( 6.0D+00 + sqrt ( 15.0D+00 ) ) / 21.0D+00 u = 0.225D+00 v = ( 155.0D+00 - sqrt ( 15.0D+00 ) ) / 1200.0D+00 w = ( 155.0D+00 + sqrt ( 15.0D+00 ) ) / 1200.0D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & a, a, & b, c, & c, b, & c, c, & d, e, & e, d, & e, e /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ u, v, v, v, w, w, w /) else if ( quad_num == 9 ) then a = 0.124949503233232D+00 b = 0.437525248383384D+00 c = 0.797112651860071D+00 d = 0.165409927389841D+00 e = 0.037477420750088D+00 u = 0.205950504760887D+00 v = 0.063691414286223D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & a, b, & b, a, & b, b, & c, d, & c, e, & d, c, & d, e, & e, c, & e, d /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ u, u, u, v, v, v, v, v, v /) else if ( quad_num == 13 ) then h = 1.0D+00 / 3.0D+00 a = 0.479308067841923D+00 b = 0.260345966079038D+00 c = 0.869739794195568D+00 d = 0.065130102902216D+00 e = 0.638444188569809D+00 f = 0.312865496004875D+00 g = 0.048690315425316D+00 w = -0.149570044467670D+00 t = 0.175615257433204D+00 u = 0.053347235608839D+00 v = 0.077113760890257D+00 quad_xy(1:dim_num,1:quad_num) = reshape ( (/ & h, h, & a, b, & b, a, & b, b, & c, d, & d, c, & d, d, & e, f, & e, g, & f, e, & f, g, & g, e, & g, f /), (/ dim_num, quad_num /) ) quad_w(1:quad_num) = (/ w, t, t, t, u, u, u, v, v, v, v, v, v /) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'QUAD_RULE - Fatal error!' write ( *, '(a,i8)' ) ' No rule is available of order QUAD_NUM = ', & quad_num stop end if return end subroutine r8_fake_use ( x ) !*****************************************************************************80 ! !! r8_fake_use() pretends to use an R8 variable. ! ! Discussion: ! ! Some compilers will issue a warning if a variable is unused. ! Sometimes there's a good reason to include a variable in a program, ! but not to use it. Calling this function with that variable as ! the argument will shut the compiler up. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 21 April 2020 ! ! Author: ! ! John Burkardt ! ! Input: ! ! real ( kind = rk ) X, the variable to be "used". ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) x if ( x /= x ) then write ( *, '(a)' ) ' r8_fake_use: variable is NAN.' end if return end subroutine r8mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_DATA_READ reads data from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Discussion: ! ! The file may contain more than N points, but this routine will ! return after reading N of them. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 October 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Output, real ( kind = rk ) TABLE(M,N), the data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer ierror character ( len = * ) input_filename integer input_status integer input_unit integer j character ( len = 255 ) line real ( kind = rk ) table(m,n) real ( kind = rk ) x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_r8vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine r8mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! R8MAT_HEADER_READ reads the header from an R8MAT file. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer M, spatial dimension. ! ! Output, integer N, the number of points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = * ) input_filename integer m integer n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, real ( kind = rk ) A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 5 integer m integer n real ( kind = rk ) a(m,n) character ( len = 14 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7,7x)') i end do write ( *, '('' Row '',5a14)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(g14.6)' ) a(i,j) end do write ( *, '(i5,1x,5a14)' ) j, ( ctemp(i), i = 1, inc ) end do end do 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: ! ! 16 October 2006 ! ! 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 subroutine reference_to_physical_t3 ( t, n, ref, phy ) !*****************************************************************************80 ! !! REFERENCE_TO_PHYSICAL_T3 maps reference points to physical points. ! ! Discussion: ! ! Given the vertices of an order 3 physical triangle and a point ! (XSI,ETA) in the reference triangle, the routine computes the value ! of the corresponding image point (X,Y) in physical space. ! ! Note that this routine may also be appropriate for an order 6 ! triangle, if the mapping between reference and physical space ! is linear. This implies, in particular, that the sides of the ! image triangle are straight and that the "midside" nodes in the ! physical triangle are literally halfway along the sides of ! the physical triangle. ! ! Reference Element T3: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 24 June 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) T(2,3), the coordinates of the vertices. ! The vertices are assumed to be the images of (0,0), (1,0) and ! (0,1) respectively. ! ! Input, integer N, the number of objects to transform. ! ! Input, real ( kind = rk ) REF(2,N), points in the reference triangle. ! ! Output, real ( kind = rk ) PHY(2,N), corresponding points in the ! physical triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i real ( kind = rk ) phy(2,n) real ( kind = rk ) ref(2,n) real ( kind = rk ) t(2,3) do i = 1, 2 phy(i,1:n) = t(i,1) * ( 1.0D+00 - ref(1,1:n) - ref(2,1:n) ) & + t(i,2) * ref(1,1:n) & + t(i,3) * ref(2,1:n) end do return end subroutine s_to_i4 ( s, ival, ierror, length ) !*****************************************************************************80 ! !! S_TO_I4 reads an I4 from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, a string to be examined. ! ! Output, integer IVAL, the integer value read from the string. ! If the string is blank, then IVAL will be returned 0. ! ! Output, integer IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, integer LENGTH, the number of characters of S ! used to make IVAL. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer i integer ierror integer isgn integer istate integer ival integer length character ( len = * ) s ierror = 0 istate = 0 isgn = 1 ival = 0 do i = 1, len_trim ( s ) c = s(i:i) ! ! Haven't read anything. ! if ( istate == 0 ) then if ( c == ' ' ) then else if ( c == '-' ) then istate = 1 isgn = -1 else if ( c == '+' ) then istate = 1 isgn = + 1 else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read the sign, expecting digits. ! else if ( istate == 1 ) then if ( c == ' ' ) then else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read at least one digit, expecting more. ! else if ( istate == 2 ) then if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then ival = 10 * ival + ichar ( c ) - ichar ( '0' ) else ival = isgn * ival length = i - 1 return end if end if end do ! ! If we read all the characters in the string, see if we're OK. ! if ( istate == 2 ) then ival = isgn * ival length = len_trim ( s ) else ierror = 1 length = 0 end if return end subroutine s_to_i4vec ( s, n, ivec, ierror ) !*****************************************************************************80 ! !! S_TO_I4VEC reads an integer vector from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 October 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer N, the number of values expected. ! ! Output, integer IVEC(N), the values read from the string. ! ! Output, integer IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer ierror integer ilo integer ivec(n) integer length character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_i4 ( s(ilo:), ivec(i), ierror, length ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + length end do return end subroutine s_to_r8 ( s, dval, ierror, length ) !*****************************************************************************80 ! !! S_TO_R8 reads an R8 from a string. ! ! Discussion: ! ! The routine will read as many characters as possible until it reaches ! the end of the string, or encounters a character which cannot be ! part of the number. ! ! Legal input is: ! ! 1 blanks, ! 2 '+' or '-' sign, ! 2.5 blanks ! 3 integer part, ! 4 decimal point, ! 5 fraction part, ! 6 'E' or 'e' or 'D' or 'd', exponent marker, ! 7 exponent sign, ! 8 exponent integer part, ! 9 exponent decimal point, ! 10 exponent fraction part, ! 11 blanks, ! 12 final comma or semicolon, ! ! with most quantities optional. ! ! Example: ! ! S DVAL ! ! '1' 1.0 ! ' 1 ' 1.0 ! '1A' 1.0 ! '12,34,56' 12.0 ! ' 34 7' 34.0 ! '-1E2ABCD' -100.0 ! '-1X2ABCD' -1.0 ! ' 2E-1' 0.2 ! '23.45' 23.45 ! '-4.2E+2' -420.0 ! '17d2' 1700.0 ! '-14e-2' -0.14 ! 'e2' 100.0 ! '-12.73e-9.23' -12.73 * 10.0**(-9.23) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string containing the ! data to be read. Reading will begin at position 1 and ! terminate at the end of the string, or when no more ! characters can be read to form a legal real. Blanks, ! commas, or other nonnumeric data will, in particular, ! cause the conversion to halt. ! ! Output, real ( kind = rk ) DVAL, the value read from the string. ! ! Output, integer IERROR, error flag. ! 0, no errors occurred. ! 1, 2, 6 or 7, the input number was garbled. The ! value of IERROR is the last type of input successfully ! read. For instance, 1 means initial blanks, 2 means ! a plus or minus sign, and so on. ! ! Output, integer LENGTH, the number of characters read ! to form the number, including any terminating ! characters such as a trailing comma or blanks. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) logical ch_eqi character c real ( kind = rk ) dval integer ierror integer ihave integer isgn integer iterm integer jbot integer jsgn integer jtop integer length integer nchar integer ndig real ( kind = rk ) rbot real ( kind = rk ) rexp real ( kind = rk ) rtop character ( len = * ) s nchar = len_trim ( s ) ierror = 0 dval = 0.0D+00 length = -1 isgn = 1 rtop = 0 rbot = 1 jsgn = 1 jtop = 0 jbot = 1 ihave = 1 iterm = 0 do length = length + 1 if ( nchar < length+1 ) then exit end if c = s(length+1:length+1) ! ! Blank character. ! if ( c == ' ' ) then if ( ihave == 2 ) then else if ( ihave == 6 .or. ihave == 7 ) then iterm = 1 else if ( 1 < ihave ) then ihave = 11 end if ! ! Comma. ! else if ( c == ',' .or. c == ';' ) then if ( ihave /= 1 ) then iterm = 1 ihave = 12 length = length + 1 end if ! ! Minus sign. ! else if ( c == '-' ) then if ( ihave == 1 ) then ihave = 2 isgn = -1 else if ( ihave == 6 ) then ihave = 7 jsgn = -1 else iterm = 1 end if ! ! Plus sign. ! else if ( c == '+' ) then if ( ihave == 1 ) then ihave = 2 else if ( ihave == 6 ) then ihave = 7 else iterm = 1 end if ! ! Decimal point. ! else if ( c == '.' ) then if ( ihave < 4 ) then ihave = 4 else if ( 6 <= ihave .and. ihave <= 8 ) then ihave = 9 else iterm = 1 end if ! ! Scientific notation exponent marker. ! else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then if ( ihave < 6 ) then ihave = 6 else iterm = 1 end if ! ! Digit. ! else if ( ihave < 11 .and. lle ( '0', c ) .and. lle ( c, '9' ) ) then if ( ihave <= 2 ) then ihave = 3 else if ( ihave == 4 ) then ihave = 5 else if ( ihave == 6 .or. ihave == 7 ) then ihave = 8 else if ( ihave == 9 ) then ihave = 10 end if call ch_to_digit ( c, ndig ) if ( ihave == 3 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = rk ) else if ( ihave == 5 ) then rtop = 10.0D+00 * rtop + real ( ndig, kind = rk ) rbot = 10.0D+00 * rbot else if ( ihave == 8 ) then jtop = 10 * jtop + ndig else if ( ihave == 10 ) then jtop = 10 * jtop + ndig jbot = 10 * jbot end if ! ! Anything else is regarded as a terminator. ! else iterm = 1 end if ! ! If we haven't seen a terminator, and we haven't examined the ! entire string, go get the next character. ! if ( iterm == 1 ) then exit end if end do ! ! If we haven't seen a terminator, and we have examined the ! entire string, then we're done, and LENGTH is equal to NCHAR. ! if ( iterm /= 1 .and. length+1 == nchar ) then length = nchar end if ! ! Number seems to have terminated. Have we got a legal number? ! Not if we terminated in states 1, 2, 6 or 7! ! if ( ihave == 1 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then ierror = ihave write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'S_TO_R8 - Serious error!' write ( *, '(a)' ) ' Illegal or nonnumeric input:' write ( *, '(a)' ) ' ' // trim ( s ) return end if ! ! Number seems OK. Form it. ! if ( jtop == 0 ) then rexp = 1.0D+00 else if ( jbot == 1 ) then rexp = 10.0D+00 ** ( jsgn * jtop ) else rexp = 10.0D+00 ** ( real ( jsgn * jtop, kind = rk ) & / real ( jbot, kind = rk ) ) end if end if dval = real ( isgn, kind = rk ) * rexp * rtop / rbot return end subroutine s_to_r8vec ( s, n, rvec, ierror ) !*****************************************************************************80 ! !! S_TO_R8VEC reads an R8VEC from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 07 September 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer N, the number of values expected. ! ! Output, real ( kind = rk ) RVEC(N), the values read from the string. ! ! Output, integer IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer ierror integer ilo integer lchar real ( kind = rk ) rvec(n) character ( len = * ) s i = 0 ierror = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_r8 ( s(ilo:), rvec(i), ierror, lchar ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + lchar end do return end subroutine s_word_count ( s, nword ) !*****************************************************************************80 ! !! S_WORD_COUNT counts the number of "words" in a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be examined. ! ! Output, integer NWORD, the number of "words" in the string. ! Words are presumed to be separated by one or more blanks. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) logical blank integer i integer lens integer nword character ( len = * ) s nword = 0 lens = len ( s ) if ( lens <= 0 ) then return end if blank = .true. do i = 1, lens if ( s(i:i) == ' ' ) then blank = .true. else if ( blank ) then nword = nword + 1 blank = .false. end if end do return end subroutine solution_write ( node_num, u, solution_file_name, time ) !*****************************************************************************80 ! !! SOLUTION_WRITE writes the solution to a file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 August 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) U(NODE_NUM), the coefficients of ! the solution. ! ! Input, character ( len = * ) SOLUTION_FILE_NAME, the name of the file ! in which the data should be stored. ! ! Input, real ( kind = rk ) TIME, the current time. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num logical, parameter :: debug = .true. integer node real ( kind = rk ), dimension(node_num) :: u character ( len = * ) :: solution_file_name integer solution_file_status integer solution_file_unit real ( kind = rk ) time call get_unit ( solution_file_unit ) open ( unit = solution_file_unit, file = solution_file_name, & status = 'replace', iostat = solution_file_status ) if ( solution_file_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SOLUTION_WRITE - Warning!' write ( *, '(a)' ) ' Could not write solution file "' & // trim ( solution_file_name ) // '" for time T = ', time return end if do node = 1, node_num write ( solution_file_unit, '(g14.6)' ) u(node) end do close ( unit = solution_file_unit ) if ( debug ) then write ( *, '(a,g14.6)' ) ' Wrote solution file "' & // trim ( solution_file_name ) // '" for time T = ', time end if return end subroutine sort_heap_external ( n, indx, i, j, isgn ) !*****************************************************************************80 ! !! SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. ! ! Discussion: ! ! The actual list of data is not passed to the routine. Hence this ! routine may be used to sort integers, reals, numbers, names, ! dates, shoe sizes, and so on. After each call, the routine asks ! the user to compare or interchange two items, until a special ! return value signals that the sorting is completed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 05 February 2004 ! ! Author: ! ! Original FORTRAN77 version by Nijenhuis and Wilf, ! FORTRAN90 version by John Burkardt ! ! Reference: ! ! A Nijenhuis and H Wilf, ! Combinatorial Algorithms, ! Academic Press, 1978, second edition, ! ISBN 0-12-519260-6. ! ! Parameters: ! ! Input, integer N, the number of items to be sorted. ! ! Input/output, integer INDX, the main communication signal. ! ! The user must set INDX to 0 before the first call. ! Thereafter, the user should not change the value of INDX until ! the sorting is done. ! ! On return, if INDX is ! ! greater than 0, ! * interchange items I and J; ! * call again. ! ! less than 0, ! * compare items I and J; ! * set ISGN = -1 if I < J, ISGN = +1 if J < I; ! * call again. ! ! equal to 0, the sorting is done. ! ! Output, integer I, J, the indices of two items. ! On return with INDX positive, elements I and J should be interchanged. ! On return with INDX negative, elements I and J should be compared, and ! the result reported in ISGN on the next call. ! ! Input, integer ISGN, results of comparison of elements ! I and J. (Used only when the previous call returned INDX less than 0). ! ISGN <= 0 means I is less than or equal to J; ! 0 <= ISGN means I is greater than or equal to J. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer, save :: i_save = 0 integer indx integer isgn integer j integer, save :: j_save = 0 integer, save :: k = 0 integer, save :: k1 = 0 integer n integer, save :: n1 = 0 ! ! INDX = 0: This is the first call. ! if ( indx == 0 ) then i_save = 0 j_save = 0 k = n / 2 k1 = k n1 = n ! ! INDX < 0: The user is returning the results of a comparison. ! else if ( indx < 0 ) then if ( indx == -2 ) then if ( isgn < 0 ) then i_save = i_save + 1 end if j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return end if if ( 0 < isgn ) then indx = 2 i = i_save j = j_save return end if if ( k <= 1 ) then if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 end if i = i_save j = j_save return end if k = k - 1 k1 = k ! ! 0 < INDX, the user was asked to make an interchange. ! else if ( indx == 1 ) then k1 = k end if do i_save = 2 * k1 if ( i_save == n1 ) then j_save = k1 k1 = i_save indx = -1 i = i_save j = j_save return else if ( i_save <= n1 ) then j_save = i_save + 1 indx = -2 i = i_save j = j_save return end if if ( k <= 1 ) then exit end if k = k - 1 k1 = k end do if ( n1 == 1 ) then i_save = 0 j_save = 0 indx = 0 i = i_save j = j_save else i_save = n1 n1 = n1 - 1 j_save = 1 indx = 1 i = i_save j = j_save end if 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 function triangle_area_2d ( t ) !*****************************************************************************80 ! !! TRIANGLE_AREA_2D computes the area of a triangle in 2D. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) T(2,3), the triangle vertices. ! ! Output, real ( kind = rk ) TRIANGLE_AREA_2D, the absolute area of ! the triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 real ( kind = rk ) t(dim_num,3) real ( kind = rk ) triangle_area_2d triangle_area_2d = 0.5D+00 * abs ( & t(1,1) * ( t(2,2) - t(2,3) ) & + t(1,2) * ( t(2,3) - t(2,1) ) & + t(1,3) * ( t(2,1) - t(2,2) ) ) return end subroutine triangulation_order6_boundary_node ( node_num, triangle_num, & triangle_node, node_boundary ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER6_BOUNDARY_NODE indicates which nodes are on the boundary. ! ! Discussion: ! ! This routine is given a triangulation, an abstract list of triples ! of nodes. It is assumed that the nodes in each triangle are listed ! in a counterclockwise order, although the routine should work ! if the nodes are consistently listed in a clockwise order as well. ! ! It is assumed that each edge of the triangulation is either ! * an INTERIOR edge, which is listed twice, once with positive ! orientation and once with negative orientation, or; ! * a BOUNDARY edge, which will occur only once. ! ! This routine should work even if the region has holes - as long ! as the boundary of the hole comprises more than 3 edges! ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 13 June 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer TRIANGLE_NUM, the number of triangles. ! ! Input, integer TRIANGLE_NODE(6,TRIANGLE_NUM), the nodes ! that make up the triangles. ! ! Output, logical NODE_BOUNDARY(NODE_NUM), is TRUE if the node ! is on a boundary edge. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer triangle_num integer e1(3*triangle_num) integer e2(3*triangle_num) integer edge(3,3*triangle_num) integer i integer m integer n logical node_boundary(node_num) integer triangle_node(6,triangle_num) m = 3 n = 3 * triangle_num ! ! Set up the edge array. The midside node is listed last, as ! it is not needed for the sorting process. ! edge(1, 1: triangle_num) = triangle_node(1,1:triangle_num) edge(2, 1: triangle_num) = triangle_node(4,1:triangle_num) edge(3, 1: triangle_num) = triangle_node(2,1:triangle_num) edge(1, triangle_num+1:2*triangle_num) = triangle_node(2,1:triangle_num) edge(2, triangle_num+1:2*triangle_num) = triangle_node(5,1:triangle_num) edge(3, triangle_num+1:2*triangle_num) = triangle_node(3,1:triangle_num) edge(1,2*triangle_num+1:3*triangle_num) = triangle_node(3,1:triangle_num) edge(2,2*triangle_num+1:3*triangle_num) = triangle_node(6,1:triangle_num) edge(3,2*triangle_num+1:3*triangle_num) = triangle_node(1,1:triangle_num) ! ! In each column, force the smaller of the two vertices to appear first. ! e1(1:n) = minval ( edge(1:3:2,1:n), dim = 1 ) e2(1:n) = maxval ( edge(1:3:2,1:n), dim = 1 ) edge(1,1:n) = e1(1:n) edge(3,1:n) = e2(1:n) ! ! Ascending sort the column array. ! call i4col_sort_a ( m, n, edge ) ! ! Records which appear twice are internal edges and can be ignored. ! node_boundary(1:node_num) = .false. i = 0 do while ( i < 3 * triangle_num ) i = i + 1 if ( i == 3 * triangle_num ) then node_boundary(edge(1:m,i)) = .true. else if ( all ( edge(1:m,i) == edge(1:m,i+1) ) ) then i = i + 1 else node_boundary(edge(1:m,i)) = .true. end if end do return end subroutine triangulation_order6_plot ( file_name, node_num, node_xy, & triangle_num, triangle_node, node_show, triangle_show ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER6_PLOT plots a 6-node triangulation of a set of nodes. ! ! 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: ! ! 27 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the output file. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. ! ! Input, integer TRIANGLE_NUM, the number of triangles. ! ! Input, integer TRIANGLE_NODE(6,TRIANGLE_NUM), lists, ! for each triangle, the indices of the nodes 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 triangle_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,triangle_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, triangle_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, triangle_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