subroutine adj_set_q4_mesh ( node_num, element_num, element_node, & element_neighbor, adj_num, adj_col, adj ) !*****************************************************************************80 ! !! adj_set_q4_mesh() sets adjacencies in a Q4 mesh. ! ! Discussion: ! ! This routine is called to set the adjacency values, after the ! appropriate amount of memory has been set aside for storage. ! ! The mesh is assumed to involve 4-node quadrilaterals. ! ! Two nodes are "adjacent" if they are both nodes in some element. ! Also, a node is considered to be adjacent to itself. ! ! This routine can be used to create the compressed column storage ! for a linear element finite element discretization of ! Poisson's equation in two dimensions. ! ! Diagram: ! ! side 3 ! 4-------3 ! s | | s ! i | | i ! d | | d ! e | | e ! | | ! 4 | | 2 ! | | ! 1-------2 ! ! side 1 ! ! The local node numbering ! ! ! 20-21-22-23-24 ! | | | | | ! | | | | | ! 15-16-17-18-19 ! | | | | | ! | | | | | ! 10-11-12-13-14 ! | | | | | ! | | | | | ! 5--6--7--8--9 ! | | | | | ! | | | | | ! 0--1--2--3--4 ! ! A sample grid. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 29 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), lists the nodes ! that make up each element in counterclockwise order. ! ! Input, integer ELEMENT_NEIGHBOR(4,ELEMENT_NUM), for each ! side of an element, lists the neighboring element, or -1 if there is ! no neighbor. ! ! Input, integer ADJ_NUM, the number of adjacencies. ! ! Input, integer ADJ_COL(NODE_NUM+1). Information about ! column J is stored in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ. ! ! Output, integer ADJ(ADJ_NUM), the adjacency information. ! implicit none integer adj_num integer node_num integer element_num integer, parameter :: element_order = 4 integer adj(adj_num) integer adj_col(node_num+1) integer adj_copy(node_num) integer k1 integer k2 integer n1 integer n2 integer n3 integer n4 integer node integer number integer element integer element2 integer element_neighbor(4,element_num) integer element_node(element_order,element_num) adj(1:adj_num) = -1 adj_copy(1:node_num) = adj_col(1:node_num) ! ! Set every node to be adjacent to itself. ! do node = 1, node_num adj(adj_copy(node)) = node adj_copy(node) = adj_copy(node) + 1 end do ! ! Examine each element. ! do element = 1, element_num n1 = element_node(1,element) n2 = element_node(2,element) n3 = element_node(3,element) n4 = element_node(4,element) ! ! Add edges (1,3) and (2,4). There is no need to check for redundancy, ! since this is the only case when these nodes can share an element. ! adj(adj_copy(n1)) = n3 adj_copy(n1) = adj_copy(n1) + 1 adj(adj_copy(n3)) = n1 adj_copy(n3) = adj_copy(n3) + 1 adj(adj_copy(n2)) = n4 adj_copy(n2) = adj_copy(n2) + 1 adj(adj_copy(n4)) = n2 adj_copy(n4) = adj_copy(n4) + 1 ! ! Add edge (1,2) if this is the first occurrence, ! that is, if the edge (1,2) is on a boundary (ELEMENT2 <= 0) ! or if this element is the first of the pair in which the edge ! occurs (ELEMENT < ELEMENT2). ! element2 = element_neighbor(1,element) if ( element2 < 0 .or. element < element2 ) then adj(adj_copy(n1)) = n2 adj_copy(n1) = adj_copy(n1) + 1 adj(adj_copy(n2)) = n1 adj_copy(n2) = adj_copy(n2) + 1 end if ! ! Add edge (2,3). ! element2 = element_neighbor(2,element) if ( element2 < 0 .or. element < element2 ) then adj(adj_copy(n2)) = n3 adj_copy(n2) = adj_copy(n2) + 1 adj(adj_copy(n3)) = n2 adj_copy(n3) = adj_copy(n3) + 1 end if ! ! Add edge (3,4). ! element2 = element_neighbor(3,element) if ( element2 < 0 .or. element < element2 ) then adj(adj_copy(n4)) = n3 adj_copy(n4) = adj_copy(n4) + 1 adj(adj_copy(n3)) = n4 adj_copy(n3) = adj_copy(n3) + 1 end if ! ! Add edge (4,1). ! element2 = element_neighbor(4,element) if ( element2 < 0 .or. element < element2 ) then adj(adj_copy(n1)) = n4 adj_copy(n1) = adj_copy(n1) + 1 adj(adj_copy(n4)) = n1 adj_copy(n4) = adj_copy(n4) + 1 end if end do ! ! Ascending sort the entries for each node. ! do node = 1, node_num k1 = adj_col(node) k2 = adj_col(node+1)-1 number = k2 + 1 - k1 call i4vec_sort_heap_a ( number, adj(k1:k2) ) end do return end subroutine adj_size_q4_mesh ( node_num, element_num, element_node, & element_neighbor, adj_num, adj_col ) !*****************************************************************************80 ! !! ADJ_SIZE_Q4_MESH counts adjacencies in a Q4 mesh. ! ! Discussion: ! ! This routine is called to count the adjacencies, so that the ! appropriate amount of memory can be set aside for storage when ! the adjacency structure is created. ! ! The mesh is assumed to involve 4-node quadrilaterals. ! ! Two nodes are "adjacent" if they are both nodes in some quadrilateral. ! Also, a node is considered to be adjacent to itself. ! ! Diagram: ! ! side 3 ! 4-------3 ! s | | s ! i | | i ! d | | d ! e | | e ! | | ! 4 | | 2 ! | | ! 1-------2 ! ! side 1 ! ! The local node numbering ! ! ! 20-21-22-23-24 ! | | | | | ! | | | | | ! 15-16-17-18-19 ! | | | | | ! | | | | | ! 10-11-12-13-14 ! | | | | | ! | | | | | ! 5--6--7--8--9 ! | | | | | ! | | | | | ! 0--1--2--3--4 ! ! A sample grid. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 29 September 2009 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), lists the ! nodes that make up each element, in counterclockwise order. ! ! Input, integer ELEMENT_NEIGHBOR(4,ELEMENT_NUM), for each ! side of a element, lists the neighboring elment, or -1 if there is ! no neighbor. ! ! Output, integer ADJ_NUM, the number of adjacencies. ! ! Output, integer ADJ_COL(NODE_NUM+1), Information about ! column J is stored in entries ADJ_COL(J) through ADJ_COL(J+1)-1 of ADJ. ! implicit none integer element_num integer, parameter :: element_order = 4 integer node_num integer adj_col(node_num+1) integer adj_num integer element integer element_neighbor(4,element_num) integer element_node(element_order,element_num) integer element2 integer n1 integer n2 integer n3 integer n4 integer node adj_num = 0 ! ! Set every node to be adjacent to itself. ! adj_col(1:node_num) = 1 ! ! Examine each element. ! do element = 1, element_num n1 = element_node(1,element) n2 = element_node(2,element) n3 = element_node(3,element) n4 = element_node(4,element) ! ! Add edge (1,3). ! adj_col(n1) = adj_col(n1) + 1 adj_col(n3) = adj_col(n3) + 1 ! ! Add edge (2,4). ! adj_col(n2) = adj_col(n2) + 1 adj_col(n4) = adj_col(n4) + 1 ! ! Add edge (1,2) if this is the first occurrence, ! that is, if the edge (1,2) is on a boundary (ELEMENT2 <= 0) ! or if this element is the first of the pair in which the edge ! occurs (ELEMENT < ELEMENT2). ! element2 = element_neighbor(1,element) if ( element2 < 0 .or. element < element2 ) then adj_col(n1) = adj_col(n1) + 1 adj_col(n2) = adj_col(n2) + 1 end if ! ! Add edge (2,3). ! element2 = element_neighbor(2,element) if ( element2 < 0 .or. element < element2 ) then adj_col(n2) = adj_col(n2) + 1 adj_col(n3) = adj_col(n3) + 1 end if ! ! Add edge (3,4). ! element2 = element_neighbor(3,element) if ( element2 < 0 .or. element < element2 ) then adj_col(n3) = adj_col(n3) + 1 adj_col(n4) = adj_col(n4) + 1 end if ! ! Add edge (4,1). ! element2 = element_neighbor(4,element) if ( element2 < 0 .or. element < element2 ) then adj_col(n4) = adj_col(n4) + 1 adj_col(n1) = adj_col(n1) + 1 end if end do ! ! We used ADJ_COL to count the number of entries in each column. ! Convert it to pointers into the ADJ array. ! do node = node_num + 1, 2, -1 adj_col(node) = adj_col(node-1) end do adj_col(1) = 1 do node = 2, node_num + 1 adj_col(node) = adj_col(node) + adj_col(node-1) end do ! ! Finally, record the total number of adjacencies. ! adj_num = adj_col(node_num+1) - 1 return end subroutine area_q4_mesh ( node_num, element_num, node_xy, element_node, & element_area, mesh_area ) !*****************************************************************************80 ! !! AREA_Q4_MESH computes areas of elements in a Q4 mesh. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the node coordinates. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), lists the ! nodes that make up each element, in counterclockwise order. ! ! Output, real ( kind = rk ) ELEMENT_AREA(ELEMENT_NUM), the element areas. ! ! Output, real ( kind = rk ) MESH_AREA, the mesh area. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer node_num integer element real ( kind = rk ) element_area(element_num) integer element_node(4,element_num) real ( kind = rk ) mesh_area integer node real ( kind = rk ) node_xy(2,node_num) real ( kind = rk ) q4(2,4) do element = 1, element_num do node = 1, 4 q4(1:2,node) = node_xy(1:2,element_node(node,element)) end do call area_quad ( q4, element_area(element) ) end do mesh_area = sum ( element_area(1:element_num) ) return end subroutine area_quad ( quad_xy, area ) !*****************************************************************************80 ! !! AREA_QUAD returns the area of a quadrilateral. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 23 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) QUAD_XY(2,4), the coordinates of the nodes. ! ! Output, real ( kind = rk ) AREA, the area of the quadrilateral. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) real ( kind = rk ) area real ( kind = rk ) area1 real ( kind = rk ) area2 real ( kind = rk ) quad_xy(2,4) real ( kind = rk ) t1(2,3) real ( kind = rk ) t2(2,3) t1(1:2,1) = quad_xy(1:2,1) t1(1:2,2) = quad_xy(1:2,2) t1(1:2,3) = quad_xy(1:2,3) call triangle_area ( t1, area1 ) t2(1:2,1) = quad_xy(1:2,3) t2(1:2,2) = quad_xy(1:2,4) t2(1:2,3) = quad_xy(1:2,1) call triangle_area ( t2, area2 ) if ( area1 < 0.0D+00 .or. area2 < 0.0D+00 ) then t1(1:2,1) = quad_xy(1:2,2) t1(1:2,2) = quad_xy(1:2,3) t1(1:2,3) = quad_xy(1:2,4) call triangle_area ( t1, area1 ) t2(1:2,1) = quad_xy(1:2,4) t2(1:2,2) = quad_xy(1:2,1) t2(1:2,3) = quad_xy(1:2,2) call triangle_area ( t2, area2 ) if ( area1 < 0.0D+00 .or. area2 < 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AREA_QUAD - Fatal error!' write ( *, '(a)' ) ' The quadrilateral nodes seem to be listed in' write ( *, '(a)' ) ' the wrong order, or the quadrilateral is' write ( *, '(a)' ) ' degenerate.' stop end if end if area = area1 + area2 return end subroutine bandwidth ( element_order, element_num, element_node, ml, mu, m ) !*****************************************************************************80 ! !! BANDWIDTH determines the bandwidth associated with a finite element mesh. ! ! Discussion: ! ! The quantity computed here is the "geometric" bandwidth determined ! by the finite element mesh alone. ! ! If a single finite element variable is associated with each node ! of the mesh, and if the nodes and variables are numbered in the ! same way, then the geometric bandwidth is the same as the bandwidth ! of a typical finite element matrix. ! ! The bandwidth M is defined in terms of the lower and upper bandwidths: ! ! M = ML + 1 + MU ! ! where ! ! ML = maximum distance from any diagonal entry to a nonzero ! entry in the same row, but earlier column, ! ! MU = maximum distance from any diagonal entry to a nonzero ! entry in the same row, but later column. ! ! Because the finite element node adjacency relationship is symmetric, ! we are guaranteed that ML = MU. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 September 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 ML, MU, the lower and upper bandwidths ! of the matrix. ! ! Output, integer M, the bandwidth of the matrix. ! implicit none 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 m integer ml integer mu ml = 0 mu = 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) mu = max ( mu, global_j - global_i ) ml = max ( ml, global_i - global_j ) end do end do end do m = ml + 1 + mu return end subroutine boundary_edge_count_q4_mesh ( element_num, element_node, & boundary_edge_num ) !*****************************************************************************80 ! !! BOUNDARY_EDGE_COUNT_Q4_MESH counts the boundary edges. ! ! Discussion: ! ! This routine is given a Q4 mesh, an abstract list of sets of 4 nodes. ! It is assumed that the nodes in each Q4 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 mesh 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. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes ! that make up the elements. These should be listed in counterclockwise ! order. ! ! Output, integer BOUNDARY_EDGE_NUM, the number of boundary ! edges. ! implicit none integer element_num integer, parameter :: element_order = 4 integer boundary_edge_num integer e1(4*element_num) integer e2(4*element_num) integer edge(2,4*element_num) integer interior_edge_num integer m integer n integer element_node(element_order,element_num) integer unique_num m = 2 n = 4 * element_num ! ! Set up the edge array. ! edge(1:2, 1: element_num) = element_node(1:2,1:element_num) edge(1:2, element_num+1:2*element_num) = element_node(2:3,1:element_num) edge(1:2,2*element_num+1:3*element_num) = element_node(3:4,1:element_num) edge(1 ,3*element_num+1:4*element_num) = element_node(4, 1:element_num) edge(2 ,3*element_num+1:4*element_num) = element_node(1, 1:element_num) ! ! In each column, force the smaller entry to appear first. ! e1(1:n) = minval ( edge(1:2,1:n), dim = 1 ) e2(1:n) = maxval ( edge(1:2,1:n), dim = 1 ) edge(1,1:n) = e1(1:n) edge(2,1:n) = e2(1:n) ! ! Ascending sort the column array. ! call i4col_sort_a ( m, n, edge ) ! ! Get the number of unique columns in EDGE. ! call i4col_sorted_unique_count ( m, n, edge, unique_num ) interior_edge_num = 4 * element_num - unique_num boundary_edge_num = 4 * element_num - 2 * interior_edge_num return end subroutine boundary_edge_count_euler_q4_mesh ( node_num, element_num, & hole_num, boundary_num ) !*****************************************************************************80 ! !! BOUNDARY_EDGE_COUNT_EULER_Q4_MESH counts boundary edges by Euler's formula. ! ! Discussion: ! ! We assume we are given information about a quadrilateral mesh ! of a set of nodes in the plane. ! ! Given the number of nodes, elements and holes, we are going to apply ! Euler's formula to determine the number of edges that lie on the ! boundary of the set of nodes. ! ! The number of faces, including the infinite face and internal holes, ! is ELEMENT_NUM + HOLE_NUM + 1. ! ! Let BOUNDARY_NUM denote the number of edges on the boundary. ! Each of the ELEMENT_NUM quadrilaterals uses four edges. Every edge ! occurs in two different elements, so the number of edges must be ! ( 4 * ELEMENT_NUM + BOUNDARY_NUM ) / 2. ! ! The number of nodes used in the mesh is NODE_NUM. ! ! Euler's formula asserts that, for a simple connected figure in the ! plane with no edge crossings, NODE_NUM nodes, EDGE_NUM edges and ! FACE_NUM faces: ! ! NODE_NUM - EDGE_NUM + FACE_NUM = 2 ! ! In our context, this becomes ! ! NODE_NUM - ( 4 * ELEMENT_NUM + BOUNDARY_NUM ) / 2 ! + ELEMENT_NUM + HOLE_NUM + 1 = 2 ! ! or ! ! BOUNDARY_NUM = 2 * NODE_NUM + 2 * HOLE_NUM - 2 * ELEMENT_NUM - 2 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 February 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Marc de Berg, Marc Krevald, Mark Overmars, Otfried Schwarzkopf, ! Computational Geometry, Section 9.1, ! Springer, 2000. ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer HOLE_NUM, the number of internal holes. ! ! Output, integer BOUNDARY_NUM, the number of edges that ! lie on the boundary of the mesh. ! implicit none integer boundary_num integer element_num integer hole_num integer node_num boundary_num = 2 * node_num + 2 * hole_num - 2 * element_num - 2 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 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 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 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 example1_q4_mesh ( node_num, element_num, node_xy, element_node, & element_neighbor ) !*****************************************************************************80 ! !! EXAMPLE1_Q4_MESH sets up example #1 Q4 mesh. ! ! Discussion: ! ! The appropriate values of NODE_NUM and ELEMENT_NUM can be found by ! calling EXAMPLE1_Q4_MESH_SIZE first. ! ! 24---25---26---27---28 ! | 14 | 15 | 16 | 17 | ! 18---19---20---21---22---23 ! | 10 | -2 | 11 | 12 | 13 | ! 12---13---14---15---16---17 ! | 5 | 6 | 7 | 8 | 9 | ! 6----7----8----9---10---11 ! | 1 | 2 | 3 | 4 | ! 1----2----3----4----5 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Output, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of the ! nodes. ! ! Output, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes ! that make up the elements. ! ! Output, integer ELEMENT_NEIGHBOR(4,ELEMENT_NUM), the ! element neighbors on each side. Negative values indicate edges that ! lie on the exterior. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 integer element_num integer, parameter :: element_order = 4 integer node_num real ( kind = rk ) node_xy(dim_num,node_num) integer element_node(element_order,element_num) integer element_neighbor(4,element_num) node_xy = reshape ( (/ & 0.0D+00, 0.0D+00, & 1.0D+00, 0.0D+00, & 2.0D+00, 0.0D+00, & 3.0D+00, 0.0D+00, & 4.0D+00, 0.0D+00, & 0.0D+00, 1.0D+00, & 1.0D+00, 1.0D+00, & 2.0D+00, 1.0D+00, & 3.0D+00, 1.0D+00, & 4.0D+00, 1.0D+00, & 5.0D+00, 1.0D+00, & 0.0D+00, 2.0D+00, & 1.0D+00, 2.0D+00, & 2.0D+00, 2.0D+00, & 3.0D+00, 2.0D+00, & 4.0D+00, 2.0D+00, & 5.0D+00, 2.0D+00, & 0.0D+00, 3.0D+00, & 1.0D+00, 3.0D+00, & 2.0D+00, 3.0D+00, & 3.0D+00, 3.0D+00, & 4.0D+00, 3.0D+00, & 5.0D+00, 3.0D+00, & 0.0D+00, 4.0D+00, & 1.0D+00, 4.0D+00, & 2.0D+00, 4.0D+00, & 3.0D+00, 4.0D+00, & 4.0D+00, 4.0D+00 /), (/ dim_num, node_num /) ) element_node(1:element_order,1:element_num ) = reshape ( (/ & 1, 2, 7, 6, & 2, 3, 8, 7, & 3, 4, 9, 8, & 4, 5, 10, 9, & 6, 7, 13, 12, & 7, 8, 14, 13, & 8, 9, 15, 14, & 9, 10, 16, 15, & 10, 11, 17, 16, & 12, 13, 19, 18, & 14, 15, 21, 20, & 15, 16, 22, 21, & 16, 17, 23, 22, & 18, 19, 25, 24, & 19, 20, 26, 25, & 20, 21, 27, 26, & 21, 22, 28, 27 /), (/ element_order, element_num /) ) element_neighbor(1:4,1:element_num) = reshape ( (/ & -1, 2, 5, -1, & -1, 3, 6, 1, & -1, 4, 7, 2, & -1, -1, 8, 3, & 1, 6, 10, -1, & 2, 7, -2, 5, & 3, 8, 11, 6, & 4, 9, 12, 7, & -1, -1, 13, 8, & 5, -2, 14, -1, & 7, 12, 16, -2, & 8, 13, 17, 11, & 9, -1, -1, 12, & 10, 15, -1, -1, & -2, 16, -1, 14, & 11, 17, -1, 15, & 12, -1, -1, 16 /), (/ 4, element_num /) ) return end subroutine example1_q4_mesh_size ( node_num, element_num, hole_num ) !*****************************************************************************80 ! !! EXAMPLE1_Q4_MESH_SIZE sets sizes for example #1 Q4 mesh. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer NODE_NUM, the number of nodes. ! ! Output, integer ELEMENT_NUM, the number of elements. ! ! Output, integer HOLE_NUM, the number of holes. ! implicit none integer element_num integer hole_num integer node_num element_num = 17 hole_num = 1 node_num = 28 return end subroutine example2_q4_mesh ( node_num, element_num, node_xy, element_node, & element_neighbor ) !*****************************************************************************80 ! !! EXAMPLE2_Q4_MESH sets up example #2 Q4 mesh. ! ! Discussion: ! ! The region is a semicircle. This example includes degenerate elements ! (the first layer of elements is touching the origin, and so has a side ! of length zero). The elements are not parallelograms. And the elements ! vary in size. ! ! Because of the treatment of node 1, algorithms for counting boundary ! edges may become "confused". ! ! The appropriate values of NODE_NUM and ELEMENT_NUM can be found by ! calling EXAMPLE1_Q4_MESH_SIZE first. ! ! 29---30---31---32---33---34---35---36---37 ! | 25 | 26 | 27 | 28 | 29 | 30 | 31 | 32 | ! 20---21---22---23---24---25---26---27---28 ! | 17 | 18 | 19 | 20 | 21 | 22 | 23 | 24 | ! 11---12---13---14---15---16---17---18---19 ! | 9 | 10 | 11 | 12 | 13 | 14 | 15 | 16 | ! 2----3----4----5----6----7----8----9---10 ! | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | ! 1----1----1----1----1----1----1----1----1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 25 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Output, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of the ! nodes. ! ! Output, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes ! that make up the elements. ! ! Output, integer ELEMENT_NEIGHBOR(4,ELEMENT_NUM), the ! element neighbors on each side. Negative values indicate edges that ! lie on the exterior. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 integer element_num integer, parameter :: element_order = 4 integer node_num real ( kind = rk ) a integer col integer element integer k integer element_node(element_order,element_num) integer element_neighbor(4,element_num) real ( kind = rk ) node_xy(dim_num,node_num) real ( kind = rk ), parameter :: pi = 3.141592653589793D+00 real ( kind = rk ) r integer row k = 1 node_xy(1,k) = 0.0D+00 node_xy(2,k) = 0.0D+00 do row = 1, 4 r = real ( row, kind = rk ) do col = 0, 8 a = real ( 8 - col, kind = rk ) * pi / 8.0D+00 k = k + 1 node_xy(1,k) = r * cos ( a ) node_xy(2,k) = r * sin ( a ) end do end do element = 0 do row = 0, 3 do col = 0, 7 element = element + 1 if ( row == 0 ) then element_node(1,element) = 1 element_node(2,element) = 1 element_node(3,element) = col + 3 element_node(4,element) = col + 2 else element_node(1,element) = element_node(4,element-8) element_node(2,element) = element_node(3,element-8) element_node(3,element) = element_node(2,element) + 9 element_node(4,element) = element_node(1,element) + 9 end if end do end do element = 0 do row = 0, 3 do col = 0, 7 element = element + 1 if ( row == 0 ) then element_neighbor(1,element) = -1 else element_neighbor(1,element) = element - 8 end if if ( col == 7 ) then element_neighbor(2,element) = -1 else element_neighbor(2,element) = element + 1 end if if ( row == 3 ) then element_neighbor(3,element) = - 1 else element_neighbor(3,element) = element + 8 end if if ( col == 0 ) then element_neighbor(4,element) = - 1 else element_neighbor(4,element) = element - 1 end if end do end do return end subroutine example2_q4_mesh_size ( node_num, element_num, hole_num ) !*****************************************************************************80 ! !! EXAMPLE2_Q4_MESH_SIZE sets sizes for example #2 Q4 mesh. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 25 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer NODE_NUM, the number of nodes. ! ! Output, integer ELEMENT_NUM, the number of elements. ! ! Output, integer HOLE_NUM, the number of holes. ! implicit none integer element_num integer hole_num integer node_num element_num = 32 hole_num = 0 node_num = 37 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 column_num logical got_one character ( len = * ) input_file_name integer input_status integer input_unit character ( len = 255 ) 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_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 bad_num integer comment_num integer ierror character ( len = * ) input_file_name integer input_status integer input_unit character ( len = 255 ) 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 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 function i4_modp ( i, j ) !*****************************************************************************80 ! !! I4_MODP returns the nonnegative remainder of I4 division. ! ! Discussion: ! ! If ! NREM = I4_MODP ( I, J ) ! NMULT = ( I - NREM ) / J ! then ! I = J * NMULT + NREM ! where NREM is always nonnegative. ! ! 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. ! ! An I4 is an integer value. ! ! 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 i integer i4_modp integer j integer value if ( j == 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4_MODP - Fatal error!' write ( *, '(a,i8)' ) ' Illegal divisor J = ', j stop end if value = mod ( i, j ) if ( value < 0 ) then value = value + abs ( j ) end if i4_modp = value return end function i4_wrap ( ival, ilo, ihi ) !*****************************************************************************80 ! !! I4_WRAP forces an I4 to lie between given limits by wrapping. ! ! Discussion: ! ! An I4 is an integer value. ! ! Example: ! ! ILO = 4, IHI = 8 ! ! I Value ! ! -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, an integer value. ! ! Input, integer ILO, IHI, the desired bounds. ! ! Output, integer I4_WRAP, a "wrapped" version of IVAL. ! implicit none integer i4_modp integer i4_wrap integer ihi integer ilo integer ival integer jhi integer jlo integer value integer wide jlo = min ( ilo, ihi ) jhi = max ( ilo, ihi ) wide = jhi - jlo + 1 if ( wide == 1 ) then value = jlo else value = jlo + i4_modp ( ival - jlo, wide ) end if i4_wrap = value return end subroutine i4col_compare ( m, n, a, i, j, isgn ) !*****************************************************************************80 ! !! I4COL_COMPARE compares columns I and J of an I4COL. ! ! 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 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.' 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.' 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 I4COL. ! ! 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 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_sorted_unique_count ( m, n, a, unique_num ) !*****************************************************************************80 ! !! I4COL_SORTED_UNIQUE_COUNT counts unique elements in an I4COL. ! ! Discussion: ! ! The columns of the array may be ascending or descending sorted. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), a sorted array, containing ! N columns of data. ! ! Output, integer UNIQUE_NUM, the number of unique columns. ! implicit none integer m integer n integer a(m,n) integer j1 integer j2 integer unique_num if ( n <= 0 ) then unique_num = 0 return end if unique_num = 1 j1 = 1 do j2 = 2, n if ( any ( a(1:m,j1) /= a(1:m,j2) ) ) then unique_num = unique_num + 1 j1 = j2 end if end do return end subroutine i4col_swap ( m, n, a, i, j ) !*****************************************************************************80 ! !! I4COL_SWAP swaps columns I and J of an I4COL. ! ! 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 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: ! ! 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 table data. ! implicit none 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. ! ! 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 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 ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed. ! ! Discussion: ! ! An I4MAT is a rectangular array of I4 values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 December 2004 ! ! 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, character ( len = * ) TITLE, a title. ! implicit none integer m integer n integer a(m,n) character ( len = * ) title call i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) 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. ! ! Discussion: ! ! An I4MAT is a rectangular array of I4 values. ! ! 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, a title. ! implicit none integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 8 ) 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) 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), '(i8)' ) i end do write ( *, '('' Row '',10a8)' ) 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), '(i8)' ) a(i,j) end do write ( *, '(i5,1x,10a8)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine i4mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! I4MAT_WRITE writes an I4MAT file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 August 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Input, integer TABLE(M,N), the table data. ! implicit none integer m integer n integer j character ( len = * ) output_filename integer output_status integer output_unit character ( len = 30 ) string integer table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a4)' ) '(', m, 'i10)' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) return end subroutine i4vec_print ( n, a, title ) !*****************************************************************************80 ! !! I4VEC_PRINT prints an I4VEC. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 November 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, integer A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title to be printed first. ! TITLE may be blank. ! implicit none integer n integer 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,i12)' ) i, a(i) end do return end subroutine i4vec_heap_d ( n, a ) !*****************************************************************************80 ! !! i4vec_heap_d reorders an I4VEC into an descending heap. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! A descending heap is an array A with the property that, for every index J, ! A(J) >= A(2*J) and A(J) >= A(2*J+1), (as long as the indices ! 2*J and 2*J+1 are legal). ! ! A(1) ! / \ ! A(2) A(3) ! / \ / \ ! A(4) A(5) A(6) A(7) ! / \ / \ ! A(8) A(9) A(10) A(11) ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 April 1999 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms for Computers and Calculators, ! Academic Press, 1978, ! ISBN: 0-12-519260-6, ! LC: QA164.N54. ! ! Input: ! ! integer N, the size of the input array. ! ! integer A(N), an unsorted array. ! ! Output: ! ! integer A(N), the array has been reordered into a heap. ! implicit none integer n integer a(n) integer i integer ifree integer key integer m ! ! Only nodes N/2 down to 1 can be "parent" nodes. ! do i = n / 2, 1, -1 ! ! Copy the value out of the parent node. ! Position IFREE is now "open". ! key = a(i) ifree = i do ! ! Positions 2*IFREE and 2*IFREE + 1 are the descendants of position ! IFREE. (One or both may not exist because they exceed N.) ! m = 2 * ifree ! ! Does the first position exist? ! if ( n < m ) then exit end if ! ! Does the second position exist? ! if ( m + 1 <= n ) then ! ! If both positions exist, take the larger of the two values, ! and update M if necessary. ! if ( a(m) < a(m+1) ) then m = m + 1 end if end if ! ! If the large descendant is larger than KEY, move it up, ! and update IFREE, the location of the free position, and ! consider the descendants of THIS position. ! if ( a(m) <= key ) then exit end if a(ifree) = a(m) ifree = m end do ! ! Once there is no more shifting to do, KEY moves into the free spot IFREE. ! a(ifree) = key end do return end subroutine i4vec_sort_heap_a ( n, a ) !*****************************************************************************80 ! !! i4vec_sort_heap_a ascending sorts an I4VEC using heap sort. ! ! Discussion: ! ! An I4VEC is a vector of I4's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 30 September 2009 ! ! Author: ! ! John Burkardt ! ! Reference: ! ! Albert Nijenhuis, Herbert Wilf, ! Combinatorial Algorithms for Computers and Calculators, ! Academic Press, 1978, ! ISBN: 0-12-519260-6, ! LC: QA164.N54. ! ! Input: ! ! integer N, the number of entries in the array. ! ! integer A(N), the array to be sorted; ! ! Parameters: ! ! Iinteger A(N), the array has been sorted. ! implicit none integer n integer a(n) integer n1 integer t if ( n <= 1 ) then return end if ! ! 1: Put A into descending heap form. ! call i4vec_heap_d ( n, a ) ! ! 2: Sort A. ! ! The largest object in the heap is in A(1). ! Move it to position A(N). ! t = a(1) a(1) = a(n) a(n) = t ! ! Consider the diminished heap of size N1. ! do n1 = n - 1, 2, -1 ! ! Restore the heap structure of A(1) through A(N1). ! call i4vec_heap_d ( n1, a ) ! ! Take the largest object from A(1) and move it to A(N1). ! t = a(1) a(1) = a(n1) a(n1) = t end do return end subroutine mesh_base_one ( node_num, element_order, element_num, element_node ) !*****************************************************************************80 ! !! MESH_BASE_ONE ensures that the element definition is one-based. ! ! Discussion: ! ! The ELEMENT_NODE array contains nodes indices that form elements. ! The convention for node indexing might start at 0 or at 1. ! Since a FORTRAN90 program will naturally assume a 1-based indexing, it is ! necessary to check a given element definition and, if it is actually ! 0-based, to convert it. ! ! This function attempts to detect 9-based node indexing and correct it. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 October 2014 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, int NODE_NUM, the number of nodes. ! ! Input, int ELEMENT_ORDER, the order of the elements. ! ! Input, int ELEMENT_NUM, the number of elements. ! ! Input/output, int ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the element ! definitions. ! implicit none integer element_num integer element_order integer element_node(element_order,element_num) integer, parameter :: i4_huge = 2147483647 integer node_max integer node_min integer node_num node_min = + i4_huge node_max = - i4_huge node_min = minval ( element_node(1:element_order,1:element_num) ) node_max = maxval ( element_node(1:element_order,1:element_num) ) if ( node_min == 0 .and. node_max == node_num - 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' )'MESH_BASE_ONE:' write ( *, '(a)' )' The element indexing appears to be 0-based!' write ( *, '(a)' )' This will be converted to 1-based.' element_node(1:element_order,1:element_num) = & element_node(1:element_order,1:element_num) + 1 else if ( node_min == 1 .and. node_max == node_num ) then write ( *, '(a)' ) ' ' write ( *, '(a)' )'MESH_BASE_ONE:' write ( *, '(a)' )' The element indexing appears to be 1-based!' write ( *, '(a)' )' No conversion is necessary.' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MESH_BASE_ONE - Warning!' write ( *, '(a)' ) ' The element indexing is not of a recognized type.' write ( *, '(a,i8)' ) ' NODE_MIN = ', node_min write ( *, '(a,i8)' ) ' NODE_MAX = ', node_max write ( *, '(a,i8)' ) ' NODE_NUM = ', node_num end if return end subroutine neighbor_elements_q4_mesh ( element_num, element_node, & element_neighbor ) !*****************************************************************************80 ! !! NEIGHBOR_ELEMENTS_Q4_MESH determines element neighbors in a Q4 mesh. ! ! Discussion: ! ! A quadrilateral mesh of a set of nodes can be completely described by ! the coordinates of the nodes, and the list of nodes that make up ! each element. However, in some cases, it is necessary to know ! element adjacency information, that is, which element, if any, ! is adjacent to a given element on a particular side. ! ! This routine creates a data structure recording this information. ! ! The primary amount of work occurs in sorting a list of 4 * ELEMENT_NUM ! data items. ! ! Note that COL is a work array allocated dynamically inside this ! routine. It is possible, for very large values of ELEMENT_NUM, ! that the necessary amount of memory will not be accessible, and the ! routine will fail. This is a limitation of the implementation of ! dynamic arrays in FORTRAN90. One way to get around this would be ! to require the user to declare ROW in the calling routine ! as an allocatable array, get the necessary memory explicitly with ! an ALLOCATE statement, and then pass ROW into this routine. ! ! Of course, the point of dynamic arrays was to make it easy to ! hide these sorts of temporary work arrays from the poor user! ! ! This routine was revised to store the edge data in a column ! array rather than a row array. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes ! that make up each element. ! ! Output, integer ELEMENT_NEIGHBOR(4,ELEMENT_NUM), lists the ! neighboring element on each side of a given element, or -1 if there is ! no neighbor. ! implicit none integer element_num integer, parameter :: element_order = 4 integer, allocatable :: col(:,:) integer element_neighbor(4,element_num) integer element_node(element_order,element_num) integer i integer icol integer j integer k integer l integer q integer q1 integer q2 integer side1 integer side2 allocate ( col (4,4*element_num) ) ! ! Step 1. ! From the list of nodes for element Q, of the form: (I,J,K,L) ! construct the four neighbor relations: ! ! (I,J,1,Q) or (J,I,1,Q), ! (J,K,2,Q) or (K,J,2,Q), ! (K,L,3,Q) or (L,K,3,Q) ! (K,I,4,Q) or (I,K,4,Q) ! ! where we choose (I,J,1,Q) if I < J, or else (J,I,1,Q) ! do q = 1, element_num i = element_node(1,q) j = element_node(2,q) k = element_node(3,q) l = element_node(4,q) if ( i < j ) then col(1:4,4*(q-1)+1) = (/ i, j, 1, q /) else col(1:4,4*(q-1)+1) = (/ j, i, 1, q /) end if if ( j < k ) then col(1:4,4*(q-1)+2) = (/ j, k, 2, q /) else col(1:4,4*(q-1)+2) = (/ k, j, 2, q /) end if if ( k < l ) then col(1:4,4*(q-1)+3) = (/ k, l, 3, q /) else col(1:4,4*(q-1)+3) = (/ l, k, 3, q /) end if if ( l < i ) then col(1:4,4*(q-1)+4) = (/ l, i, 4, q /) else col(1:4,4*(q-1)+4) = (/ i, l, 4, q /) end if end do ! ! Step 2. Perform an ascending dictionary sort on the neighbor relations. ! We only intend to sort on rows 1 and 2; the routine we call here ! sorts on rows 1 through 4 but that won't hurt us. ! ! What we need is to find cases where two elements share an edge. ! Say they share an edge defined by the nodes I and J. Then there are ! two columns of COL that start out ( I, J, ?, ? ). By sorting COL, ! we make sure that these two columns occur consecutively. That will ! make it easy to notice that the elements are neighbors. ! call i4col_sort_a ( 4, 4*element_num, col ) ! ! Step 3. Neighboring elements show up as consecutive columns with ! identical first two entries. Whenever you spot this happening, ! make the appropriate entries in ELEMENT_NEIGHBOR. ! element_neighbor(1:4,1:element_num) = -1 icol = 1 do if ( 4 * element_num <= icol ) then exit end if if ( col(1,icol) /= col(1,icol+1) .or. col(2,icol) /= col(2,icol+1) ) then icol = icol + 1 cycle end if side1 = col(3,icol) q1 = col(4,icol) side2 = col(3,icol+1) q2 = col(4,icol+1) element_neighbor(side1,q1) = q2 element_neighbor(side2,q2) = q1 icol = icol + 2 end do deallocate ( col ) return end subroutine node_order_q4_mesh ( element_num, element_node, node_num, & node_order ) !*****************************************************************************80 ! !! NODE_ORDER_Q4_MESH determines the order of nodes in a Q4 mesh. ! ! Discussion: ! ! The order of a node is the number of elements that use that node ! as a vertex. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), ! the nodes that make up the elements. ! ! Input, integer NODE_NUM, the number of nodes. ! ! Output, integer NODE_ORDER(NODE_NUM), the order of each node. ! implicit none integer element_num integer, parameter :: element_order = 4 integer node_num integer element integer element_node(element_order,element_num) integer i integer node integer node_order(node_num) node_order(1:node_num) = 0 do element = 1, element_num do i = 1, element_order node = element_node(i,element) if ( node < 1 .or. node_num < node ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NODE_ORDER_Q4_MESH - Fatal error!' write ( *, '(a)' ) ' Illegal entry in ELEMENT_NODE.' stop else node_order(node) = node_order(node) + 1 end if end do end do return end subroutine plot_q4_mesh ( node_num, element_num, node_xy, element_node, & node_show, element_show, output_filename ) !*****************************************************************************80 ! !! PLOT_Q4_MESH plots a Q4 mesh. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes ! that form the elements. ! ! Input, integer NODE_SHOW, ! 0, do not show nodes; ! 1, show nodes; ! 2, show nodes and label them. ! ! Input, integer ELEMENT_SHOW, ! 0, do not show elements; ! 1, show elements; ! 2, show elements and label them. ! ! Input, character ( len = * ) OUTPUT_FILENAME, the name of the output file. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer node_num integer element_num integer, parameter :: element_order = 4 real ( kind = rk ) ave_x real ( kind = rk ) ave_y integer :: circle_size integer delta integer e integer element integer element_node(element_order,element_num) integer element_show integer i integer i4_wrap integer ios integer node integer node_show real ( kind = rk ) node_xy(2,node_num) character ( len = * ) output_filename integer output_unit 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 ( output_unit ) open ( unit = output_unit, file = output_filename, status = 'replace', & iostat = ios ) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'PLOT_Q4_MESH - Fatal error!' write ( *, '(a)' ) ' Can not open output file.' return end if write ( output_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( output_unit, '(a)' ) '%%Creator: plot_q4_mesh.f90' write ( output_unit, '(a)' ) '%%Title: ' // trim ( output_filename ) write ( output_unit, '(a)' ) '%%Pages: 1' write ( output_unit, '(a,i3,2x,i3,2x,i3,2x,i3)' ) '%%BoundingBox: ', & x_ps_min, y_ps_min, x_ps_max, y_ps_max write ( output_unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( output_unit, '(a)' ) '%%LanguageLevel: 1' write ( output_unit, '(a)' ) '%%EndComments' write ( output_unit, '(a)' ) '%%BeginProlog' write ( output_unit, '(a)' ) '/inch {72 mul} def' write ( output_unit, '(a)' ) '%%EndProlog' write ( output_unit, '(a)' ) '%%Page: 1 1' write ( output_unit, '(a)' ) 'save' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the RGB line color to very light gray.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '0.900 0.900 0.900 setrgbcolor' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Draw a gray border around the page.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) 'newpath' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_min, ' moveto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_max, y_ps_min, ' lineto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_max, y_ps_max, ' lineto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_max, ' lineto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', x_ps_min, y_ps_min, ' lineto' write ( output_unit, '(a)' ) 'stroke' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the RGB color to black.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '0.000 0.000 0.000 setrgbcolor' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the font and its size.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '/Times-Roman findfont' write ( output_unit, '(a)' ) '0.50 inch scalefont' write ( output_unit, '(a)' ) 'setfont' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Print a title.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% 210 702 moveto' write ( output_unit, '(a)' ) '% (Q4 Mesh) show' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Define a clipping polygon.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) 'newpath' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_min_clip, ' moveto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_max_clip, y_ps_min_clip, ' lineto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_max_clip, y_ps_max_clip, ' lineto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_max_clip, ' lineto' write ( output_unit, '(a,i3,2x,i3,2x,a)' ) ' ', & x_ps_min_clip, y_ps_min_clip, ' lineto' write ( output_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 ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Draw filled dots at the nodes.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the RGB color to blue.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '0.000 0.150 0.750 setrgbcolor' write ( output_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 ( output_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 ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Label the nodes:' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the RGB color to darker blue.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '0.000 0.250 0.850 setrgbcolor' write ( output_unit, '(a)' ) '/Times-Roman findfont' write ( output_unit, '(a)' ) '0.20 inch scalefont' write ( output_unit, '(a)' ) 'setfont' write ( output_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 ( string, '(i4)' ) node string = adjustl ( string ) write ( output_unit, '(i4,2x,i4,a)' ) x_ps, y_ps+5, & ' moveto (' // trim ( string ) // ') show' end do end if ! ! Draw the elements. ! if ( 1 <= element_show ) then write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the RGB color to red.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '0.900 0.200 0.100 setrgbcolor' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Draw the elements.' write ( output_unit, '(a)' ) '%' do element = 1, element_num write ( output_unit, '(a)' ) 'newpath' do i = 1, element_order + 1 e = i4_wrap ( i, 1, element_order ) node = element_node(e,element) 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 ) ) if ( i == 1 ) then write ( output_unit, '(i3,2x,i3,2x,a)' ) x_ps, y_ps, ' moveto' else write ( output_unit, '(i3,2x,i3,2x,a)' ) x_ps, y_ps, ' lineto' end if end do write ( output_unit, '(a)' ) 'stroke' end do end if ! ! Label the elements. ! if ( 2 <= element_show ) then write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Label the elements:' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% Set the RGB color to darker red.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '0.950 0.250 0.150 setrgbcolor' write ( output_unit, '(a)' ) '/Times-Roman findfont' write ( output_unit, '(a)' ) '0.20 inch scalefont' write ( output_unit, '(a)' ) 'setfont' write ( output_unit, '(a)' ) '%' do element = 1, element_num ave_x = 0.0D+00 ave_y = 0.0D+00 do i = 1, element_order node = element_node(i,element) ave_x = ave_x + node_xy(1,node) ave_y = ave_y + node_xy(2,node) end do ave_x = ave_x / real ( element_order, kind = rk ) ave_y = ave_y / real ( element_order, kind = rk ) 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)' ) element string = adjustl ( string ) write ( output_unit, '(i4,2x,i4,a)' ) x_ps, y_ps, ' moveto (' & // trim ( string ) // ') show' end do end if write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) 'restore showpage' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '% End of page.' write ( output_unit, '(a)' ) '%' write ( output_unit, '(a)' ) '%%Trailer' write ( output_unit, '(a)' ) '%%EOF' close ( unit = output_unit ) return end subroutine r8mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_DATA_READ reads data from an R8MAT file. ! ! 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 table 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. ! ! 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 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 ( m, n, a, title ) !*****************************************************************************80 ! !! R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! 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, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n real ( kind = rk ) a(m,n) character ( len = * ) title call r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) 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. ! ! Discussion: ! ! An R8MAT is an array of R8 values. ! ! 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, a 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 write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) 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), '(i8,6x)' ) 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 r8mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! R8MAT_WRITE writes an R8MAT file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Input, real ( kind = rk ) TABLE(M,N), the table data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer j character ( len = * ) output_filename integer output_status integer output_unit character ( len = 30 ) string real ( kind = rk ) table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'R8MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! ! For greater precision in the output file, try: ! ! '(', m, 'g', 24, '.', 16, ')' ! if ( 0 < m .and. 0 < n ) then write ( string, '(a1,i8,a1,i8,a1,i8,a1)' ) '(', m, 'g', 14, '.', 6, ')' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do end if ! ! Close the file. ! close ( unit = output_unit ) return end subroutine r8vec_bracket ( n, x, xval, left, right ) !*****************************************************************************80 ! !! R8VEC_BRACKET searches a sorted R8VEC for successive brackets of a value. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! If the values in the vector are thought of as defining intervals ! on the real line, then this routine searches for the interval ! nearest to or containing the given value. ! ! It is always true that RIGHT = LEFT+1. ! ! If XVAL < X(1), then LEFT = 1, RIGHT = 2, and ! XVAL < X(1) < X(2); ! If X(1) <= XVAL < X(N), then ! X(LEFT) <= XVAL < X(RIGHT); ! If X(N) <= XVAL, then LEFT = N-1, RIGHT = N, and ! X(LEFT) <= X(RIGHT) <= XVAL. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 06 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, length of input array. ! ! Input, real ( kind = rk ) X(N), an array that has been sorted into ! ascending order. ! ! Input, real ( kind = rk ) XVAL, a value to be bracketed. ! ! Output, integer LEFT, RIGHT, the results of the search. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer left integer right real ( kind = rk ) x(n) real ( kind = rk ) xval do i = 2, n - 1 if ( xval < x(i) ) then left = i - 1 right = i return end if end do left = n - 1 right = n return end subroutine r8vec_print ( n, a, title ) !*****************************************************************************80 ! !! R8VEC_PRINT prints an R8VEC. ! ! Discussion: ! ! An R8VEC is a vector of R8's. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 August 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer N, the number of components of the vector. ! ! Input, real ( kind = rk ) A(N), the vector to be printed. ! ! Input, character ( len = * ) TITLE, a title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) a(n) integer i character ( len = * ) title write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) write ( *, '(a)' ) ' ' do i = 1, n write ( *, '(2x,i8,2x,g16.8)' ) i, a(i) end do return end subroutine reference_to_physical_q4 ( q4, n, rs, xy ) !*****************************************************************************80 ! !! REFERENCE_TO_PHYSICAL_Q4 maps Q4 reference points to physical points. ! ! Discussion: ! ! XY(R,S) = XY(0,0) * (1-R) * (1-S) ! + XY(1,0) * R * (1-S) ! + XY(1,1) * R * S ! + XY(0,1) * (1-R) * S ! ! Reference Element Q4: ! ! | ! 1 4-----3 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1-----2 ! | ! +--0--R--1--> ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 25 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) Q4(2,4), the coordinates of the vertices. ! The vertices are assumed to be the images of the reference vertices ! (0,0), (1,0), (1,1) and (0,1) respectively. ! ! Input, integer N, the number of points to transform. ! ! Input, real ( kind = rk ) RS(2,N), (R,S) points in the reference element. ! ! Output, real ( kind = rk ) XY(2,N), (X,Y) points in the physical element. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) psi(4,n) real ( kind = rk ) q4(2,4) real ( kind = rk ) rs(2,n) real ( kind = rk ) xy(2,n) psi(1,1:n) = ( 1.0D+00 - rs(1,1:n) ) * ( 1.0D+00 - rs(2,1:n) ) psi(2,1:n) = rs(1,1:n) * ( 1.0D+00 - rs(2,1:n) ) psi(3,1:n) = rs(1,1:n) * rs(2,1:n) psi(4,1:n) = ( 1.0D+00 - rs(1,1:n) ) * rs(2,1:n) xy(1:2,1:n) = matmul ( q4(1:2,1:4), psi(1:4,1:n) ) return end subroutine s_blank_delete ( s ) !*****************************************************************************80 ! !! S_BLANK_DELETE removes blanks from a string, left justifying the remainder. ! ! Discussion: ! ! All TAB characters are also removed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 26 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character ( len = * ) S, the string to be transformed. ! implicit none character c integer get integer put integer nchar character ( len = * ) s character, parameter :: TAB = char ( 9 ) put = 0 nchar = len_trim ( s ) do get = 1, nchar c = s(get:get) if ( c /= ' ' .and. c /= TAB ) then put = put + 1 s(put:put) = c end if end do s(put+1:nchar) = ' ' 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 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 I4VEC 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 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 ) character c logical ch_eqi 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 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 sample_q4_mesh ( node_num, node_xy, element_num, element_node, & sample_num, sample_xy, sample_element ) !*****************************************************************************80 ! !! SAMPLE_Q4_MESH returns random points in a Q4 mesh. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 15 March 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NODE_NUM, the number of nodes. ! ! Input, real ( kind = rk ) NODE_XY(2,NODE_NUM), the coordinates of the nodes. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer ELEMENT_NODE(4,ELEMENT_NUM), the nodes ! that form the elements. ! ! Input, integer SAMPLE_NUM, the number of points to sample. ! ! Output, real ( kind = rk ) SAMPLE_XY(2,SAMPLE_NUM), the sample points. ! ! Output, integer SAMPLE_ELEMENT(SAMPLE_NUM), the elements from ! which each point was drawn. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer node_num integer sample_num real ( kind = rk ) area real ( kind = rk ) area_cum(0:element_num) real ( kind = rk ) area_total integer element integer element_node(4,element_num) integer i1 integer i2 integer i3 integer i4 integer left real ( kind = rk ) node_xy(2,node_num) real ( kind = rk ) quad_xy(2,4) real ( kind = rk ) r integer right integer sample integer sample_element(sample_num) real ( kind = rk ) sample_xy(2,sample_num) ! ! Compute the areas of the quadrilaterals. ! area_cum(0) = 0.0D+00 do element = 1, element_num i1 = element_node(1,element) i2 = element_node(2,element) i3 = element_node(3,element) i4 = element_node(4,element) quad_xy(1:2,1) = node_xy(1:2,i1) quad_xy(1:2,2) = node_xy(1:2,i2) quad_xy(1:2,3) = node_xy(1:2,i3) quad_xy(1:2,4) = node_xy(1:2,i4) call area_quad ( quad_xy, area ) area_cum(element) = area_cum(element-1) + area end do area_total = area_cum(element_num) area_cum(0:element_num) = area_cum(0:element_num) / area_total ! ! A random value R indicates the corresponding quadrilateral whose ! cumulative relative area first includes the number R. ! do sample = 1, sample_num call random_number ( harvest = r ) call r8vec_bracket ( element_num + 1, area_cum, r, left, right ) element = right - 1 i1 = element_node(1,element) i2 = element_node(2,element) i3 = element_node(3,element) i4 = element_node(4,element) quad_xy(1:2,1) = node_xy(1:2,i1) quad_xy(1:2,2) = node_xy(1:2,i2) quad_xy(1:2,3) = node_xy(1:2,i3) quad_xy(1:2,4) = node_xy(1:2,i4) call sample_quad ( quad_xy, 1, sample_xy(1:2,sample) ) sample_element(sample) = element end do return end subroutine sample_quad ( quad_xy, n, xy ) !*****************************************************************************80 ! !! SAMPLE_QUAD returns random points in a quadrilateral. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 February 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) QUAD_XY(2,4), the coordinates of the nodes. ! ! Input, integer N, the number of points to sample. ! ! Output, real ( kind = rk ) XY(2,N), the sample points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n real ( kind = rk ) area1 real ( kind = rk ) area2 real ( kind = rk ) area_total integer i real ( kind = rk ) quad_xy(2,4) real ( kind = rk ) r real ( kind = rk ) t1(2,3) real ( kind = rk ) t2(2,3) real ( kind = rk ) xy(2,n) t1(1:2,1) = quad_xy(1:2,1) t1(1:2,2) = quad_xy(1:2,2) t1(1:2,3) = quad_xy(1:2,3) call triangle_area ( t1, area1 ) t2(1:2,1) = quad_xy(1:2,3) t2(1:2,2) = quad_xy(1:2,4) t2(1:2,3) = quad_xy(1:2,1) call triangle_area ( t2, area2 ) if ( area1 < 0.0D+00 .or. area2 < 0.0D+00 ) then t1(1:2,1) = quad_xy(1:2,2) t1(1:2,2) = quad_xy(1:2,3) t1(1:2,3) = quad_xy(1:2,4) call triangle_area ( t1, area1 ) t2(1:2,1) = quad_xy(1:2,4) t2(1:2,2) = quad_xy(1:2,1) t2(1:2,3) = quad_xy(1:2,2) call triangle_area ( t2, area2 ) if ( area1 < 0.0D+00 .or. area2 < 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SAMPLE_QUAD - Fatal error!' write ( *, '(a)' ) ' The quadrilateral nodes seem to be listed in' write ( *, '(a)' ) ' the wrong order, or the quadrilateral is' write ( *, '(a)' ) ' degenerate.' stop end if end if area_total = area1 + area2 ! ! Choose a triangle at random, weighted by the areas. ! Then choose a point in that triangle. ! do i = 1, n call random_number ( harvest = r ) if ( r * area_total < area1 ) then call triangle_sample ( t1, 1, xy(1:2,i) ) else call triangle_sample ( t2, 1, xy(1:2,i) ) end if end do 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 Albert Nijenhuis, Herbert Wilf. ! FORTRAN90 version by John Burkardt. ! ! Reference: ! ! Albert Nijenhuis, Herbert 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 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: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! implicit none character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine triangle_area ( t, area ) !*****************************************************************************80 ! !! TRIANGLE_AREA computes the area of a triangle. ! ! Discussion: ! ! If the triangle's vertices are given in counter clockwise order, ! the area will be positive. If the triangle's vertices are given ! in clockwise order, the area will be negative! ! ! An earlier version of this routine always returned the absolute ! value of the computed area. I am convinced now that that is ! a less useful result! For instance, by returning the signed ! area of a triangle, it is possible to easily compute the area ! of a nonconvex polygon as the sum of the (possibly negative) ! areas of triangles formed by node 1 and successive pairs of vertices. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 17 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) T(2,3), the triangle vertices. ! ! Output, real ( kind = rk ) AREA, the area of the triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 real ( kind = rk ) area real ( kind = rk ) t(dim_num,3) area = 0.5D+00 * ( & 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 triangle_sample ( t, n, p ) !*****************************************************************************80 ! !! TRIANGLE_SAMPLE returns random points in a triangle. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 10 April 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = rk ) T(2,3), the triangle vertices. ! ! Input, integer N, the number of points to generate. ! ! Output, real ( kind = rk ) P(2,N), random points in the triangle. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: dim_num = 2 integer n real ( kind = rk ) alpha(n) integer dim real ( kind = rk ) p(dim_num,n) real ( kind = rk ) p12(dim_num,n) real ( kind = rk ) p13(dim_num,n) real ( kind = rk ) t(dim_num,3) call random_number ( harvest = alpha(1:n) ) ! ! Interpret R as a percentage of the triangle's area. ! ! Imagine a line L, parallel to side 1, so that the area between ! vertex 1 and line L is R percent of the full triangle's area. ! ! The line L will intersect sides 2 and 3 at a fraction ! ALPHA = SQRT ( R ) of the distance from vertex 1 to vertices 2 and 3. ! alpha(1:n) = sqrt ( alpha(1:n) ) ! ! Determine the coordinates of the points on sides 2 and 3 intersected ! by line L. ! do dim = 1, dim_num p12(dim,1:n) = ( 1.0D+00 - alpha(1:n) ) * t(dim,1) & + alpha(1:n) * t(dim,2) p13(dim,1:n) = ( 1.0D+00 - alpha(1:n) ) * t(dim,1) & + alpha(1:n) * t(dim,3) end do ! ! Now choose, uniformly at random, a point on the line L. ! call random_number ( harvest = alpha(1:n) ) do dim = 1, dim_num p(dim,1:n) = ( 1.0D+00 - alpha(1:n) ) * p12(dim,1:n) & + alpha(1:n) * p13(dim,1:n) end do return end