subroutine bandwidth_mesh ( element_order, element_num, element_node, & ml, mu, m ) !*****************************************************************************80 ! !! BANDWIDTH_MESH: bandwidth of 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. ! ! Modified: ! ! 02 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Output, integer ( kind = 4 ) ML, MU, the lower and upper bandwidths of the matrix. ! ! Output, integer ( kind = 4 ) M, the bandwidth of the matrix. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) element_order integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) global_i integer ( kind = 4 ) global_j integer ( kind = 4 ) local_i integer ( kind = 4 ) local_j integer ( kind = 4 ) m integer ( kind = 4 ) ml integer ( kind = 4 ) 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 bandwidth_var ( element_order, element_num, element_node, & node_num, var_node, var_num, var, ml, mu, m ) !*****************************************************************************80 ! !! BANDWIDTH_VAR determines the bandwidth for finite element variables. ! ! Discussion: ! ! We assume that, attached to each node in the finite element mesh ! there are a (possibly zero) number of finite element variables. ! We wish to determine the bandwidth necessary to store the stiffness ! matrix associated with these variables. ! ! An entry K(I,J) of the stiffness matrix must be zero unless the ! variables I and J correspond to nodes N(I) and N(J) which are ! common to some element. ! ! In order to determine the bandwidth of the stiffness matrix, we ! essentially seek a nonzero entry K(I,J) for which abs ( I - J ) ! is maximized. ! ! 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. ! ! We assume the finite element variable adjacency relationship is ! symmetric, so we are guaranteed that ML = MU. ! ! Note that the user is free to number the variables in any way ! whatsoever, and to associate variables to nodes in any way, ! so that some nodes have no variables, some have one, and some ! have several. ! ! The storage of the indices of the variables is fairly simple. ! In VAR, simply list all the variables associated with node 1, ! then all those associated with node 2, and so on. Then set up ! the pointer array VAR_NODE so that we can jump to the section of ! VAR where the list begins for any particular node. ! ! The routine does not check that each variable is only associated ! with a single node. This would normally be the case in a finite ! element setting. ! ! Modified: ! ! 03 September 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM); ! ELEMENT_NODE(I,J) is the global index of local node I in element J. ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, integer ( kind = 4 ) VAR_NODE(NODE_NUM+1), used to find the variables ! associated with a given node, which are in VAR in locations ! VAR_NODE(NODE) to VAR_NODE(NODE+1)-1. Note that the last entry of ! this array points to the location just after the last location in VAR. ! ! Input, integer ( kind = 4 ) VAR_NUM, the number of variables. ! ! Input, integer ( kind = 4 ) VAR(VAR_NUM), the indexes of the variables, which ! are presumably (but not necessarily) 1, 2, 3, ..., VAR_NUM. ! ! Output, integer ( kind = 4 ) ML, MU, the lower and upper bandwidths of the matrix. ! ! Output, integer ( kind = 4 ) M, the bandwidth of the matrix. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) element_order integer ( kind = 4 ) node_num integer ( kind = 4 ) var_num integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) m integer ( kind = 4 ) ml integer ( kind = 4 ) mu integer ( kind = 4 ) node_global_i integer ( kind = 4 ) node_global_j integer ( kind = 4 ) node_local_i integer ( kind = 4 ) node_local_j integer ( kind = 4 ) var(var_num) integer ( kind = 4 ) var_global_i integer ( kind = 4 ) var_global_j integer ( kind = 4 ) var_local_i integer ( kind = 4 ) var_local_j integer ( kind = 4 ) var_node(node_num+1) ml = 0 mu = 0 do element = 1, element_num do node_local_i = 1, element_order node_global_i = element_node(node_local_i,element) do var_local_i = var_node(node_global_i), var_node(node_global_i+1)-1 var_global_i = var(var_local_i) do node_local_j = 1, element_order node_global_j = element_node(node_local_j,element) do var_local_j = var_node(node_global_j), var_node(node_global_j+1)-1 var_global_j = var(var_local_j) mu = max ( mu, var_global_j - var_global_i ) ml = max ( ml, var_global_i - var_global_j ) end do end do end do end do end do m = ml + 1 + mu return end subroutine basis_11_t3 ( t, i, p, qi, dqidx, dqidy ) !*****************************************************************************80 ! !! BASIS_11_T3: one basis at one point for the T3 element. ! ! Discussion: ! ! The routine is given the coordinates of the nodes of a triangle. ! ! 3 ! / \ ! / \ ! / \ ! 1-------2 ! ! It evaluates the linear basis function Q(I)(X,Y) associated with ! node I, which has the property that it is a linear function ! which is 1 at node I and zero at the other two nodes. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,3), the coordinates of the nodes. ! ! Input, integer ( kind = 4 ) I, the index of the desired basis function. ! I should be between 1 and 3. ! ! Input, real ( kind = 8 ) P(2), the coordinates of a point at which the ! basis function is to be evaluated. ! ! Output, real ( kind = 8 ) QI, DQIDX, DQIDY, the values of the basis ! function and its X and Y derivatives. ! implicit none real ( kind = 8 ) area real ( kind = 8 ) dqidx real ( kind = 8 ) dqidy integer ( kind = 4 ) i integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) ip1 integer ( kind = 4 ) ip2 real ( kind = 8 ) p(2) real ( kind = 8 ) qi real ( kind = 8 ) t(2,3) area = abs ( t(1,1) * ( t(2,2) - t(2,3) ) & + t(1,2) * ( t(2,3) - t(2,1) ) & + t(1,3) * ( t(2,1) - t(2,2) ) ) if ( area == 0.0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_11_T3 - Fatal error!' write ( *, '(a)' ) ' Element has zero area.' stop end if if ( i < 1 .or. 3 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_11_T3 - Fatal error!' write ( *, '(a)' ) ' Basis index I is not between 1 and 3.' write ( *, '(a,i8)' ) ' I = ', i stop end if ip1 = i4_wrap ( i + 1, 1, 3 ) ip2 = i4_wrap ( i + 2, 1, 3 ) qi = ( ( t(1,ip2) - t(1,ip1) ) * ( p(2) - t(2,ip1) ) & - ( t(2,ip2) - t(2,ip1) ) * ( p(1) - t(1,ip1) ) ) / area dqidx = - ( t(2,ip2) - t(2,ip1) ) / area dqidy = ( t(1,ip2) - t(1,ip1) ) / area return end subroutine basis_11_t3_test ( ) !*****************************************************************************80 ! !! BASIS_11_T3_TEST verifies BASIS_11_T3. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer ( kind = 4 ), parameter :: node_num = 3 real ( kind = 8 ) dphidx(node_num,node_num) real ( kind = 8 ) dphidy(node_num,node_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) phi(node_num,node_num) real ( kind = 8 ) sum_x real ( kind = 8 ) sum_y real ( kind = 8 ), dimension ( 2, node_num ) :: t = reshape ( (/ & 2.0D+00, 0.0D+00, & 4.0D+00, 3.0D+00, & 0.0D+00, 4.0D+00 /), (/ 2, node_num /) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_11_T3_TEST:' write ( *, '(a)' ) ' Verify basis functions for element T3.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Physical Nodes:' write ( *, '(a)' ) ' ' do j = 1, node_num write ( *, '(2x,i8,2x,f7.3,2x,f7.3)' ) j, t(1:2,j) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' do i = 1, node_num do j = 1, node_num call basis_11_t3 ( t, i, t(1:2,j), phi(i,j), dphidx(i,j), dphidy(i,j) ) end do end do do i = 1, node_num write ( *, '(2x,10f7.3)' ) phi(i,1:node_num) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The X and Y derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dPhidX sum dPhidY sum' write ( *, '(a)' ) ' ' do j = 1, node_num sum_x = sum ( dphidx(1:node_num,j) ) sum_y = sum ( dphidy(1:node_num,j) ) write ( *, '(2x,f14.8,f14.8)' ) sum_x, sum_y end do return end subroutine basis_11_t6 ( t, i, p, bi, dbidx, dbidy ) !*****************************************************************************80 ! !! BASIS_11_T6: one basis at one point for the T6 element. ! ! Discussion: ! ! The routine is given the coordinates of the nodes of a triangle. ! ! 3 ! / \ ! 6 5 ! / \ ! 1---4---2 ! ! It evaluates the quadratic basis function B(I)(X,Y) associated with ! node I, which has the property that it is a quadratic function ! which is 1 at node I and zero at the other five nodes. ! ! This routine assumes that the sides of the triangle are straight, ! so that the midside nodes fall on the line between two vertices. ! ! This routine relies on the fact that each basis function can be ! written as the product of two linear factors, which are easily ! computed and normalized. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,6), the coordinates of the nodes. ! ! Input, integer ( kind = 4 ) I, the index of the desired basis function. ! I should be between 1 and 6. ! ! Input, real ( kind = 8 ) P(2), the coordinates of a point at which the ! basis function is to be evaluated. ! ! Output, real ( kind = 8 ) BI, DBIDX, DBIDY, the values of the basis ! function and its X and Y derivatives. ! implicit none real ( kind = 8 ) bi real ( kind = 8 ) dbidx real ( kind = 8 ) dbidy real ( kind = 8 ) gf real ( kind = 8 ) gn real ( kind = 8 ) hf real ( kind = 8 ) hn integer ( kind = 4 ) i integer ( kind = 4 ) i4_wrap integer ( kind = 4 ) j1 integer ( kind = 4 ) j2 integer ( kind = 4 ) k1 integer ( kind = 4 ) k2 real ( kind = 8 ) p(2) real ( kind = 8 ) t(2,6) if ( i < 1 .or. 6 < i ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_11_T6 - Fatal error!' write ( *, '(a)' ) ' Basis index I is not between 1 and 6.' write ( *, '(a,i8)' ) ' I = ', i stop end if ! ! Determine the pairs of nodes. ! if ( i <= 3 ) then j1 = i4_wrap ( i + 1, 1, 3 ) j2 = i4_wrap ( i + 2, 1, 3 ) k1 = i + 3 k2 = i4_wrap ( i + 5, 4, 6 ) else j1 = i - 3 j2 = i4_wrap ( i - 3 + 2, 1, 3 ) k1 = i4_wrap ( i - 3 + 1, 1, 3 ) k2 = i4_wrap ( i - 3 + 2, 1, 3 ) end if ! ! Evaluate the two linear factors GF and HF, ! and their normalizers GN and HN. ! gf = ( p(1) - t(1,j1) ) * ( t(2,j2) - t(2,j1) ) & - ( t(1,j2) - t(1,j1) ) * ( p(2) - t(2,j1) ) gn = ( t(1,i) - t(1,j1) ) * ( t(2,j2) - t(2,j1) ) & - ( t(1,j2) - t(1,j1) ) * ( t(2,i) - t(2,j1) ) hf = ( p(1) - t(1,k1) ) * ( t(2,k2) - t(2,k1) ) & - ( t(1,k2) - t(1,k1) ) * ( p(2) - t(2,k1) ) hn = ( t(1,i) - t(1,k1) ) * ( t(2,k2) - t(2,k1) ) & - ( t(1,k2) - t(1,k1) ) * ( t(2,i) - t(2,k1) ) ! ! Construct the basis function and its derivatives. ! bi = ( gf / gn ) * ( hf / hn ) dbidx = ( ( t(2,j2) - t(2,j1) ) / gn ) * ( hf / hn ) & + ( gf / gn ) * ( ( t(2,k2) - t(2,k1) ) / hn ) dbidy = - ( ( t(1,j2) - t(1,j1) ) / gn ) * ( hf / hn ) & - ( gf / gn ) * ( ( t(1,k2) - t(1,k1) ) / hn ) return end subroutine basis_11_t6_test ( ) !*****************************************************************************80 ! !! BASIS_11_T6_TEST verifies BASIS_11_T6. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer ( kind = 4 ), parameter :: node_num = 6 real ( kind = 8 ) dphidx(node_num,node_num) real ( kind = 8 ) dphidy(node_num,node_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) phi(node_num,node_num) real ( kind = 8 ) sum_x real ( kind = 8 ) sum_y real ( kind = 8 ), dimension ( 2, node_num ) :: t = reshape ( (/ & 2.0D+00, 0.0D+00, & 4.0D+00, 3.0D+00, & 0.0D+00, 4.0D+00, & 3.0D+00, 1.5D+00, & 2.0D+00, 3.5D+00, & 1.0D+00, 2.0D+00 /), (/ 2, node_num /) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_11_T6_TEST:' write ( *, '(a)' ) ' Verify basis functions for element T6.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Physical Nodes:' write ( *, '(a)' ) ' ' do j = 1, node_num write ( *, '(2x,i8,2x,f7.3,2x,f7.3)' ) j, t(1:2,j) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' do i = 1, node_num do j = 1, node_num call basis_11_t6 ( t, i, t(1:2,j), phi(i,j), dphidx(i,j), dphidy(i,j) ) end do end do do i = 1, node_num write ( *, '(2x,10f7.3)' ) phi(i,1:node_num) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The X and Y derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dPhidX sum dPhidY sum' write ( *, '(a)' ) ' ' do j = 1, node_num sum_x = sum ( dphidx(1:node_num,j) ) sum_y = sum ( dphidy(1:node_num,j) ) write ( *, '(2x,f14.8,2x,f14.8)' ) sum_x, sum_y end do return end subroutine basis_mn_q4 ( q, n, p, phi, dphidx, dphidy ) !*****************************************************************************80 ! !! BASIS_MN_Q4: all bases at N points for a Q4 element. ! ! Discussion: ! ! The routine is given the coordinates of the vertices of a quadrilateral. ! It works directly with these coordinates, and does not refer to a ! reference element. ! ! The sides of the element are presumed to lie along coordinate axes. ! ! The routine evaluates the basis functions, and their X and Y derivatives. ! ! Physical Element Q4: ! ! | ! | 4-----3 ! | | | ! | | | ! Y | | ! | | | ! | | | ! | 1-----2 ! | ! +-----X------> ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) Q(2,4), the coordinates of the vertices. ! It is common to list these points in counter clockwise order. ! ! Input, integer ( kind = 4 ) N, the number of evaluation points. ! ! Input, real ( kind = 8 ) P(2,N), the evaluation points. ! ! Output, real ( kind = 8 ) PHI(4,N), the basis functions ! at the evaluation points. ! ! Output, real ( kind = 8 ) DPHIDX(4,N), DPHIDY(4,N), the basis ! derivatives at the evaluation points. ! ! Local Parameter: ! ! Local, real ( kind = 8 ) AREA, the area of the rectangle. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) area real ( kind = 8 ) dphidx(4,n) real ( kind = 8 ) dphidy(4,n) real ( kind = 8 ) p(2,n) real ( kind = 8 ) phi(4,n) real ( kind = 8 ) q(2,4) area = ( q(1,3) - q(1,1) ) & * ( q(2,3) - q(2,1) ) phi(1,1:n) = ( q(1,3) - p(1,1:n) ) & * ( q(2,3) - p(2,1:n) ) phi(2,1:n) = ( p(1,1:n) - q(1,1) ) & * ( q(2,3) - p(2,1:n) ) phi(3,1:n) = ( p(1,1:n) - q(1,1) ) & * ( p(2,1:n) - q(2,1) ) phi(4,1:n) = ( q(1,3) - p(1,1:n) ) & * ( p(2,1:n) - q(2,1) ) dphidx(1,1:n) = - ( q(2,3) - p(2,1:n) ) dphidx(2,1:n) = ( q(2,3) - p(2,1:n) ) dphidx(3,1:n) = ( p(2,1:n) - q(2,1) ) dphidx(4,1:n) = - ( p(2,1:n) - q(2,1) ) dphidy(1,1:n) = - ( q(1,3) - p(1,1:n) ) dphidy(2,1:n) = - ( p(1,1:n) - q(1,1) ) dphidy(3,1:n) = ( p(1,1:n) - q(1,1) ) dphidy(4,1:n) = ( q(1,3) - p(1,1:n) ) ! ! Normalize. ! phi(1:4,1:n) = phi(1:4,1:n) / area dphidx(1:4,1:n) = dphidx(1:4,1:n) / area dphidy(1:4,1:n) = dphidy(1:4,1:n) / area return end subroutine basis_mn_q4_test ( ) !*****************************************************************************80 ! !! BASIS_MN_Q4_TEST verifies BASIS_MN_Q4. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer ( kind = 4 ), parameter :: node_num = 4 real ( kind = 8 ) dphidx(node_num,node_num) real ( kind = 8 ) dphidy(node_num,node_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) phi(node_num,node_num) real ( kind = 8 ), dimension ( 2, node_num ) :: q = reshape ( (/ & 3.0D+00, 1.0D+00, & 5.0D+00, 1.0D+00, & 5.0D+00, 4.0D+00, & 3.0D+00, 4.0D+00 /), (/ 2, node_num /) ) real ( kind = 8 ) sum_x real ( kind = 8 ) sum_y write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' BASIS_MN_Q4_TEST: ' write ( *, '(a)' ) ' Verify basis functions for element Q4.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Physical Nodes:' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' I X Y' write ( *, '(a)' ) ' ' do j = 1, node_num write ( *, '(2x,i8,2x,f7.3,2x,f7.3)' ) j, q(1:2,j) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' call basis_mn_q4 ( q, node_num, q, phi, dphidx, dphidy ) do i = 1, node_num write ( *, '(2x,10f7.3)' ) phi(i,1:node_num) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The X and Y derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dPhidX sum dPhidY sum' write ( *, '(a)' ) ' ' do j = 1, node_num sum_x = sum ( dphidx(1:node_num,j) ) sum_y = sum ( dphidy(1:node_num,j) ) write ( *, '(2x,f14.8,2x,f14.8)' ) sum_x, sum_y end do return end subroutine basis_mn_t3 ( t, n, p, phi, dphidx, dphidy ) !*****************************************************************************80 ! !! BASIS_MN_T3: all bases at N points for a T3 element. ! ! Discussion: ! ! The routine is given the coordinates of the vertices of a triangle. ! It works directly with these coordinates, and does not refer to a ! reference element. ! ! The sides of the triangle DO NOT have to lie along a coordinate ! axis. ! ! The routine evaluates the basis functions associated with each vertex, ! and their derivatives with respect to X and Y. ! ! Physical Element T3: ! ! 3 ! / \ ! / \ ! / \ ! / \ ! 1---------2 ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,3), the coordinates of the vertices ! of the triangle. It is common to list these points in counter clockwise ! order. ! ! Input, integer ( kind = 4 ) N, the number of evaluation points. ! ! Input, real ( kind = 8 ) P(2,N), the points where the basis functions ! are to be evaluated. ! ! Output, real ( kind = 8 ) PHI(3,N), the value of the basis functions ! at the evaluation points. ! ! Output, real ( kind = 8 ) DPHIDX(3,N), DPHIDY(3,N), the value of the ! derivatives at the evaluation points. ! ! Local parameters: ! ! Local, real ( kind = 8 ) AREA, is (twice) the area of the triangle. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) area real ( kind = 8 ) dphidx(3,n) real ( kind = 8 ) dphidy(3,n) real ( kind = 8 ) p(2,n) real ( kind = 8 ) phi(3,n) real ( kind = 8 ) t(2,3) area = 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) ) if ( area == 0.0D+00 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_MN_T3 - Fatal error!' write ( *, '(a)' ) ' Element has zero area.' stop end if phi(1,1:n) = ( ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) ) & - ( t(2,3) - t(2,2) ) * ( p(1,1:n) - t(1,2) ) ) dphidx(1,1:n) = - ( t(2,3) - t(2,2) ) dphidy(1,1:n) = ( t(1,3) - t(1,2) ) phi(2,1:n) = ( ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) ) & - ( t(2,1) - t(2,3) ) * ( p(1,1:n) - t(1,3) ) ) dphidx(2,1:n) = - ( t(2,1) - t(2,3) ) dphidy(2,1:n) = ( t(1,1) - t(1,3) ) phi(3,1:n) = ( ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) ) & - ( t(2,2) - t(2,1) ) * ( p(1,1:n) - t(1,1) ) ) dphidx(3,1:n) = - ( t(2,2) - t(2,1) ) dphidy(3,1:n) = ( t(1,2) - t(1,1) ) ! ! Normalize. ! phi(1:3,1:n) = phi(1:3,1:n) / area dphidx(1:3,1:n) = dphidx(1:3,1:n) / area dphidy(1:3,1:n) = dphidy(1:3,1:n) / area return end subroutine basis_mn_t3_test ( ) !*****************************************************************************80 ! !! BASIS_MN_T3_TEST verifies BASIS_MN_T3. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer ( kind = 4 ), parameter :: node_num = 3 real ( kind = 8 ) dphidx(node_num,node_num) real ( kind = 8 ) dphidy(node_num,node_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) phi(node_num,node_num) real ( kind = 8 ) sum_x real ( kind = 8 ) sum_y real ( kind = 8 ), dimension ( 2, node_num ) :: t = reshape ( (/ & 2.0D+00, 0.0D+00, & 4.0D+00, 3.0D+00, & 0.0D+00, 4.0D+00 /), (/ 2, node_num /) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_MN_T3_TEST:' write ( *, '(a)' ) ' Verify basis functions for element T3.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Physical Nodes:' write ( *, '(a)' ) ' ' do j = 1, node_num write ( *, '(2x,i8,2x,f7.3,2x,f7.3)' ) j, t(1:2,j) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' call basis_mn_t3 ( t, node_num, t, phi, dphidx, dphidy ) do i = 1, node_num write ( *, '(2x,10f7.3)' ) phi(i,1:node_num) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The X and Y derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dPhidX sum dPhidY sum' write ( *, '(a)' ) ' ' do j = 1, node_num sum_x = sum ( dphidx(1:node_num,j) ) sum_y = sum ( dphidy(1:node_num,j) ) write ( *, '(2x,f14.8,2x,f14.8)' ) sum_x, sum_y end do return end subroutine basis_mn_t4 ( t, n, p, phi, dphidx, dphidy ) !*****************************************************************************80 ! !! BASIS_MN_T4: all bases at N points for a T4 element. ! ! Discussion: ! ! The T4 element is the cubic bubble triangle. ! ! The routine is given the coordinates of the vertices of a triangle. ! It works directly with these coordinates, and does not refer to a ! reference element. ! ! The sides of the triangle DO NOT have to lie along a coordinate ! axis. ! ! The routine evaluates the basis functions associated with each vertex, ! and their derivatives with respect to X and Y. ! ! Physical Element T4: ! ! 3 ! / \ ! / \ ! / 4 \ ! / \ ! 1---------2 ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,4), the coordinates of the vertices ! of the triangle, and the coordinates of the centroid. ! It is common to list the first three points in counter clockwise ! order. ! ! Input, integer ( kind = 4 ) N, the number of evaluation points. ! ! Input, real ( kind = 8 ) P(2,N), the points where the basis functions ! are to be evaluated. ! ! Output, real ( kind = 8 ) PHI(4,N), the value of the basis functions ! at the evaluation points. ! ! Output, real ( kind = 8 ) DPHIDX(4,N), DPHIDY(4,N), the value of the ! derivatives at the evaluation points. ! ! Local parameters: ! ! Local, real ( kind = 8 ) AREA, is (twice) the area of the triangle. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) area real ( kind = 8 ) dphidx(4,n) real ( kind = 8 ) dphidy(4,n) integer ( kind = 4 ) i real ( kind = 8 ) p(2,n) real ( kind = 8 ) phi(4,n) real ( kind = 8 ) t(2,4) area = 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) ) phi(1,1:n) = ( ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) ) & - ( t(2,3) - t(2,2) ) * ( p(1,1:n) - t(1,2) ) ) dphidx(1,1:n) = - ( t(2,3) - t(2,2) ) dphidy(1,1:n) = ( t(1,3) - t(1,2) ) phi(2,1:n) = ( ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) ) & - ( t(2,1) - t(2,3) ) * ( p(1,1:n) - t(1,3) ) ) dphidx(2,1:n) = - ( t(2,1) - t(2,3) ) dphidy(2,1:n) = ( t(1,1) - t(1,3) ) phi(3,1:n) = ( ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) ) & - ( t(2,2) - t(2,1) ) * ( p(1,1:n) - t(1,1) ) ) dphidx(3,1:n) = - ( t(2,2) - t(2,1) ) dphidy(3,1:n) = ( t(1,2) - t(1,1) ) ! ! Normalize the first three functions. ! phi(1:3,1:n) = phi(1:3,1:n) / area dphidx(1:3,1:n) = dphidx(1:3,1:n) / area dphidy(1:3,1:n) = dphidy(1:3,1:n) / area ! ! Compute the cubic bubble function. ! phi(4,1:n) = 27.0D+00 * phi(1,1:n) * phi(2,1:n) * phi(3,1:n) dphidx(4,1:n) = 27.0D+00 * ( & dphidx(1,1:n) * phi(2,1:n) * phi(3,1:n) & + phi(1,1:n) * dphidx(2,1:n) * phi(3,1:n) & + phi(1,1:n) * phi(2,1:n) * dphidx(3,1:n) ) dphidy(4,1:n) = 27.0D+00 * ( & dphidy(1,1:n) * phi(2,1:n) * phi(3,1:n) & + phi(1,1:n) * dphidy(2,1:n) * phi(3,1:n) & + phi(1,1:n) * phi(2,1:n) * dphidy(3,1:n) ) ! ! Subtract 1/3 of the cubic bubble function from each of the three linears. ! do i = 1, 3 phi(i,1:n) = phi(i,1:n) - phi(4,1:n) / 3.0D+00 dphidx(i,1:n) = dphidx(i,1:n) - dphidx(4,1:n) / 3.0D+00 dphidy(i,1:n) = dphidy(i,1:n) - dphidy(4,1:n) / 3.0D+00 end do return end subroutine basis_mn_t4_test ( ) !*****************************************************************************80 ! !! BASIS_MN_T4_TEST verifies BASIS_MN_T4. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer ( kind = 4 ), parameter :: node_num = 4 real ( kind = 8 ) dphidx(node_num,node_num) real ( kind = 8 ) dphidy(node_num,node_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) phi(node_num,node_num) real ( kind = 8 ) sum_x real ( kind = 8 ) sum_y real ( kind = 8 ), dimension ( 2, node_num ) :: t = reshape ( (/ & 2.0D+00, 0.0D+00, & 4.0D+00, 2.0D+00, & 0.0D+00, 4.0D+00, & 2.0D+00, 2.0D+00 /), (/ 2, node_num /) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_MN_T4_TEST:' write ( *, '(a)' ) ' Verify basis functions for element T4.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Physical Nodes:' write ( *, '(a)' ) ' ' do j = 1, node_num write ( *, '(2x,i8,2x,f7.3,2x,f7.3)' ) j, t(1:2,j) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' call basis_mn_t4 ( t, node_num, t, phi, dphidx, dphidy ) do i = 1, node_num write ( *, '(2x,10f7.3)' ) phi(i,1:node_num) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The X and Y derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dPhidX sum dPhidY sum' write ( *, '(a)' ) ' ' do j = 1, node_num sum_x = sum ( dphidx(1:node_num,j) ) sum_y = sum ( dphidy(1:node_num,j) ) write ( *, '(2x,f14.8,2x,f14.8)' ) sum_x, sum_y end do return end subroutine basis_mn_t6 ( t, n, p, phi, dphidx, dphidy ) !*****************************************************************************80 ! !! BASIS_MN_T6: all bases at N points for a T6 element. ! ! Discussion: ! ! The routine is given the coordinates of the vertices and midside ! nodes of a triangle. It works directly with these coordinates, and does ! not refer to a reference element. ! ! This routine requires that the midside nodes be "in line" ! with the vertices, that is, that the sides of the triangle be ! straight. However, the midside nodes do not actually have to ! be halfway along the side of the triangle. ! ! Physical element T6: ! ! This picture indicates the assumed ordering of the six nodes ! of the triangle. ! ! 3 ! / \ ! / \ ! 6 5 ! / \ ! / \ ! 1-----4-----2 ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) T(2,6), the nodal oordinates of the element. ! It is common to list these points in counter clockwise order. ! ! Input, integer ( kind = 4 ) N, the number of evaluation points. ! ! Input, real ( kind = 8 ) P(2,N), the coordinates of the point where ! the basis functions are to be evaluated. ! ! Output, real ( kind = 8 ) PHI(6,N), the basis functions at the ! evaluation points. ! ! Output, real ( kind = 8 ) DPHIDX(6,N), DPHIDY(6,N), the derivatives ! of the basis functions at the evaluation points. ! ! Local Parameters: ! ! Local, real ( kind = 8 ) AREA, is (twice) the area of the triangle. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) dphidx(6,n) real ( kind = 8 ) dphidy(6,n) real ( kind = 8 ) gn(n) real ( kind = 8 ) gx(n) real ( kind = 8 ) hn(n) real ( kind = 8 ) hx(n) real ( kind = 8 ) p(2,n) real ( kind = 8 ) phi(6,n) real ( kind = 8 ) t(2,6) ! ! Basis function 1: PHI(X,Y) = G(3,2) * H(6,4) / normalization. ! gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3) - t(2,2) ) & - ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) ) gn(1:n) = ( t(1,1) - t(1,2) ) * ( t(2,3) - t(2,2) ) & - ( t(1,3) - t(1,2) ) * ( t(2,1) - t(2,2) ) hx(1:n) = ( p(1,1:n) - t(1,4) ) * ( t(2,6) - t(2,4) ) & - ( t(1,6) - t(1,4) ) * ( p(2,1:n) - t(2,4) ) hn(1:n) = ( t(1,1) - t(1,4) ) * ( t(2,6) - t(2,4) ) & - ( t(1,6) - t(1,4) ) * ( t(2,1) - t(2,4) ) phi(1,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) ) dphidx(1,1:n) = ( ( t(2,3) - t(2,2) ) * hx(1:n) & + gx(1:n) * ( t(2,6) - t(2,4) ) ) / ( gn(1:n) * hn(1:n) ) dphidy(1,1:n) = -( ( t(1,3) - t(1,2) ) * hx(1:n) & + gx(1:n) * ( t(1,6) - t(1,4) ) ) / ( gn(1:n) * hn(1:n) ) ! ! Basis function 2: PHI(X,Y) = G(3,1) * H(4,5) / normalization. ! gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3) - t(2,1) ) & - ( t(1,3) - t(1,1) ) * ( p(2,1:n) - t(2,1) ) gn(1:n) = ( t(1,2) - t(1,1) ) * ( t(2,3) - t(2,1) ) & - ( t(1,3) - t(1,1) ) * ( t(2,2) - t(2,1) ) hx(1:n) = ( p(1,1:n) - t(1,5) ) * ( t(2,4) - t(2,5) ) & - ( t(1,4) - t(1,5) ) * ( p(2,1:n) - t(2,5) ) hn(1:n) = ( t(1,2) - t(1,5) ) * ( t(2,4) - t(2,5) ) & - ( t(1,4) - t(1,5) ) * ( t(2,2) - t(2,5) ) phi(2,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) ) dphidx(2,1:n) = ( ( t(2,3) - t(2,1) ) * hx(1:n) & + gx(1:n) * ( t(2,4) - t(2,5) ) ) / ( gn(1:n) * hn(1:n) ) dphidy(2,1:n) = -( ( t(1,3) - t(1,1) ) * hx(1:n) & + gx(1:n) * ( t(1,4) - t(1,5) ) ) / ( gn(1:n) * hn(1:n) ) ! ! Basis function 3: PHI(X,Y) = G(1,2) * H(5,6) / normalization. ! gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1) - t(2,2) ) & - ( t(1,1) - t(1,2) ) * ( p(2,1:n) - t(2,2) ) gn(1:n) = ( t(1,3) - t(1,2) ) * ( t(2,1) - t(2,2) ) & - ( t(1,1) - t(1,2) ) * ( t(2,3) - t(2,2) ) hx(1:n) = ( p(1,1:n) - t(1,6) ) * ( t(2,5) - t(2,6) ) & - ( t(1,5) - t(1,6) ) * ( p(2,1:n) - t(2,6) ) hn(1:n) = ( t(1,3) - t(1,6) ) * ( t(2,5) - t(2,6) ) & - ( t(1,5) - t(1,6) ) * ( t(2,3) - t(2,6) ) phi(3,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) ) dphidx(3,1:n) = ( ( t(2,1) - t(2,2) ) * hx(1:n) & + gx(1:n) * ( t(2,5) - t(2,6) ) ) / ( gn(1:n) * hn(1:n) ) dphidy(3,1:n) = -( ( t(1,1) - t(1,2) ) * hx(1:n) & + gx(1:n) * ( t(1,5) - t(1,6) ) ) / ( gn(1:n) * hn(1:n) ) ! ! Basis function 4: PHI(X,Y) = G(1,3) * H(2,3) / normalization. ! gx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,1) - t(2,3) ) & - ( t(1,1) - t(1,3) ) * ( p(2,1:n) - t(2,3) ) gn(1:n) = ( t(1,4) - t(1,3) ) * ( t(2,1) - t(2,3) ) & - ( t(1,1) - t(1,3) ) * ( t(2,4) - t(2,3) ) hx(1:n) = ( p(1,1:n) - t(1,3) ) * ( t(2,2) - t(2,3) ) & - ( t(1,2) - t(1,3) ) * ( p(2,1:n) - t(2,3) ) hn(1:n) = ( t(1,4) - t(1,3) ) * ( t(2,2) - t(2,3) ) & - ( t(1,2) - t(1,3) ) * ( t(2,4) - t(2,3) ) phi(4,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) ) dphidx(4,1:n) = ( ( t(2,1) - t(2,3) ) * hx(1:n) & + gx(1:n) * ( t(2,2) - t(2,3) ) ) / ( gn(1:n) * hn(1:n) ) dphidy(4,1:n) = -( ( t(1,1) - t(1,3) ) * hx(1:n) & + gx(1:n) * ( t(1,2) - t(1,3) ) ) / ( gn(1:n) * hn(1:n) ) ! ! Basis function 5: PHI(X,Y) = G(2,1) * H(3,1) / normalization. ! gx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,2) - t(2,1) ) & - ( t(1,2) - t(1,1) ) * ( p(2,1:n) - t(2,1) ) gn(1:n) = ( t(1,5) - t(1,1) ) * ( t(2,2) - t(2,1) ) & - ( t(1,2) - t(1,1) ) * ( t(2,5) - t(2,1) ) hx(1:n) = ( p(1,1:n) - t(1,1) ) * ( t(2,3) - t(2,1) ) & - ( t(1,3) - t(1,1) ) * ( p(2,1:n) - t(2,1) ) hn(1:n) = ( t(1,5) - t(1,1) ) * ( t(2,3) - t(2,1) ) & - ( t(1,3) - t(1,1) ) * ( t(2,5) - t(2,1) ) phi(5,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) ) dphidx(5,1:n) = ( ( t(2,2) - t(2,1) ) * hx(1:n) & + gx(1:n) * ( t(2,3) - t(2,1) ) ) / ( gn(1:n) * hn(1:n) ) dphidy(5,1:n) = -( ( t(1,2) - t(1,1) ) * hx(1:n) & + gx(1:n) * ( t(1,3) - t(1,1) ) ) / ( gn(1:n) * hn(1:n) ) ! ! Basis function 6: PHI(X,Y) = G(1,2) * H(3,2) / normalization. ! gx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,1) - t(2,2) ) & - ( t(1,1) - t(1,2) ) * ( p(2,1:n) - t(2,2) ) gn(1:n) = ( t(1,6) - t(1,2) ) * ( t(2,1) - t(2,2) ) & - ( t(1,1) - t(1,2) ) * ( t(2,6) - t(2,2) ) hx(1:n) = ( p(1,1:n) - t(1,2) ) * ( t(2,3) - t(2,2) ) & - ( t(1,3) - t(1,2) ) * ( p(2,1:n) - t(2,2) ) hn(1:n) = ( t(1,6) - t(1,2) ) * ( t(2,3) - t(2,2) ) & - ( t(1,3) - t(1,2) ) * ( t(2,6) - t(2,2) ) phi(6,1:n) = ( gx(1:n) * hx(1:n) ) / ( gn(1:n) * hn(1:n) ) dphidx(6,1:n) = ( ( t(2,1) - t(2,2) ) * hx(1:n) & + gx(1:n) * ( t(2,3) - t(2,2) ) ) / ( gn(1:n) * hn(1:n) ) dphidy(6,1:n) = -( ( t(1,1) - t(1,2) ) * hx(1:n) & + gx(1:n) * ( t(1,3) - t(1,2) ) ) / ( gn(1:n) * hn(1:n) ) return end subroutine basis_mn_t6_test ( ) !*****************************************************************************80 ! !! BASIS_MN_T6_TEST verifies BASIS_MN_T6. ! ! Modified: ! ! 08 February 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer ( kind = 4 ), parameter :: node_num = 6 real ( kind = 8 ) dphidx(node_num,node_num) real ( kind = 8 ) dphidy(node_num,node_num) integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) phi(node_num,node_num) real ( kind = 8 ) sum_x real ( kind = 8 ) sum_y real ( kind = 8 ), dimension ( 2, node_num ) :: t = reshape ( (/ & 2.0D+00, 0.0D+00, & 4.0D+00, 3.0D+00, & 0.0D+00, 4.0D+00, & 3.0D+00, 1.5D+00, & 2.0D+00, 3.5D+00, & 1.0D+00, 2.0D+00 /), (/ 2, node_num /) ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'BASIS_MN_T6_TEST:' write ( *, '(a)' ) ' Verify basis functions for element T6.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of nodes = ', node_num write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Physical Nodes:' write ( *, '(a)' ) ' ' do j = 1, node_num write ( *, '(2x,i8,2x,f7.3,2x,f7.3)' ) j, t(1:2,j) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' call basis_mn_t6 ( t, node_num, t, phi, dphidx, dphidy ) do i = 1, node_num write ( *, '(2x,10f7.3)' ) phi(i,1:node_num) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' The X and Y derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dPhidX sum dPhidY sum' write ( *, '(a)' ) ' ' do j = 1, node_num sum_x = sum ( dphidx(1:node_num,j) ) sum_y = sum ( dphidy(1:node_num,j) ) write ( *, '(2x,2f14.8)' ) sum_x, sum_y end do return end subroutine ch_cap ( ch ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character CH, the character to capitalize. ! implicit none character ch integer ( kind = 4 ) itemp itemp = ichar ( ch ) if ( 97 <= itemp .and. itemp <= 122 ) then ch = char ( itemp - 32 ) end if return end function degrees_to_radians ( angle ) !*****************************************************************************80 ! !! DEGREES_TO_RADIANS converts an angle from degrees to radians. ! ! Modified: ! ! 10 July 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real ( kind = 8 ) ANGLE, an angle in degrees. ! ! Output, real ( kind = 8 ) DEGREES_TO_RADIANS, the equivalent angle ! in radians. ! implicit none real ( kind = 8 ) angle real ( kind = 8 ) degrees_to_radians real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 degrees_to_radians = ( angle / 180.0D+00 ) * pi return end subroutine derivative_average_t3 ( node_num, node_xy, element_num, & element_node, c, dcdx, dcdy ) !*****************************************************************************80 ! !! DERIVATIVE_AVERAGE_T3 averages derivatives at the nodes of a T3 mesh. ! ! Discussion: ! ! This routine can be used to compute an averaged nodal value of any ! quantity associated with the finite element function. At a node ! that is shared by several elements, the fundamental function ! U will be continuous, but its spatial derivatives, for instance, ! will generally be discontinuous. This routine computes the ! value of the spatial derivatives in each element, and averages ! them, to make a reasonable assignment of a nodal value. ! ! In this version of the routine, the average is not weighted. ! ! Modified: ! ! 09 June 2006 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of ! the nodes. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(3,ELEMENT_NUM), ! the element->node data. ! ! Input, real ( kind = 8 ) C(NODE_NUM), the finite element coefficient ! vector. ! ! Output, real ( kind = 8 ) DCDX(NODE_NUM), DCDY(NODE_NUM), the averaged ! values of dCdX and dCdY at the nodes. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ), parameter :: element_order = 3 integer ( kind = 4 ) node_num real ( kind = 8 ) c(node_num) real ( kind = 8 ) dcdx(node_num) real ( kind = 8 ) dcdy(node_num) real ( kind = 8 ) dphidx(element_order,element_order) real ( kind = 8 ) dphidy(element_order,element_order) integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) node_count(node_num) integer ( kind = 4 ) node_global1 integer ( kind = 4 ) node_global2 integer ( kind = 4 ) node_local1 integer ( kind = 4 ) node_local2 real ( kind = 8 ) node_xy(2,node_num) real ( kind = 8 ) phi(element_order,element_order) real ( kind = 8 ) t(2,element_order) node_count(1:node_num) = 0 dcdx(1:node_num) = 0.0D+00 dcdy(1:node_num) = 0.0D+00 ! ! Consider every element. ! do element = 1, element_num ! ! Get the coordinates of the nodes of the element. ! t(1:2,1:element_order) = node_xy(1:2,element_node(1:element_order,element)) ! ! Evaluate the X and Y derivatives of the basis functions at the nodes. ! call basis_mn_t3 ( t, element_order, t, phi, dphidx, dphidy ) ! ! Evaluate dCdX and dCdY at each node in the element, and add ! them to the running totals. ! do node_local1 = 1, element_order node_global1 = element_node(node_local1,element) do node_local2 = 1, element_order node_global2 = element_node(node_local2,element) dcdx(node_global1) = dcdx(node_global1) & + c(node_global2) * dphidx(node_local2,node_local1) dcdy(node_global1) = dcdy(node_global1) & + c(node_global2) * dphidy(node_local2,node_local1) end do node_count(node_global1) = node_count(node_global1) + 1 end do end do ! ! Average the running totals. ! dcdx(1:node_num) = dcdx(1:node_num) & / real ( node_count(1:node_num), kind = 8 ) dcdy(1:node_num) = dcdy(1:node_num) & / real ( node_count(1:node_num), kind = 8 ) return end subroutine derivative_average_t6 ( node_num, node_xy, element_num, & element_node, c, dcdx, dcdy ) !*****************************************************************************80 ! !! DERIVATIVE_AVERAGE_T6 averages derivatives at the nodes of a T6 mesh. ! ! Discussion: ! ! This routine can be used to compute an averaged nodal value of any ! quantity associated with the finite element function. At a node ! that is shared by several elements, the fundamental function ! U will be continuous, but its spatial derivatives, for instance, ! will generally be discontinuous. This routine computes the ! value of the spatial derivatives in each element, and averages ! them, to make a reasonable assignment of a nodal value. ! ! In this version of the routine, the average is not weighted. ! ! Modified: ! ! 22 April 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the coordinates of ! the nodes. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(6,ELEMENT_NUM), ! the element->node data. ! ! Input, real ( kind = 8 ) C(NODE_NUM), the finite element coefficient ! vector. ! ! Output, real ( kind = 8 ) DCDX(NODE_NUM), DCDY(NODE_NUM), the averaged ! values of dCdX and dCdY at the nodes. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ), parameter :: element_order = 6 integer ( kind = 4 ) node_num real ( kind = 8 ) c(node_num) real ( kind = 8 ) dcdx(node_num) real ( kind = 8 ) dcdy(node_num) real ( kind = 8 ) dphidx(element_order,element_order) real ( kind = 8 ) dphidy(element_order,element_order) integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) node_count(node_num) integer ( kind = 4 ) node_global1 integer ( kind = 4 ) node_global2 integer ( kind = 4 ) node_local1 integer ( kind = 4 ) node_local2 real ( kind = 8 ) node_xy(2,node_num) real ( kind = 8 ) phi(element_order,element_order) real ( kind = 8 ) t(2,element_order) node_count(1:node_num) = 0 dcdx(1:node_num) = 0.0D+00 dcdy(1:node_num) = 0.0D+00 ! ! Consider every element. ! do element = 1, element_num ! ! Get the coordinates of the nodes of the element. ! t(1:2,1:element_order) = node_xy(1:2,element_node(1:element_order,element)) ! ! Evaluate the X and Y derivatives of the basis functions at the nodes. ! call basis_mn_t6 ( t, element_order, t, phi, dphidx, dphidy ) ! ! Evaluate dCdX and dCdY at each node in the element, and add ! them to the running totals. ! do node_local1 = 1, element_order node_global1 = element_node(node_local1,element) do node_local2 = 1, element_order node_global2 = element_node(node_local2,element) dcdx(node_global1) = dcdx(node_global1) & + c(node_global2) * dphidx(node_local2,node_local1) dcdy(node_global1) = dcdy(node_global1) & + c(node_global2) * dphidy(node_local2,node_local1) end do node_count(node_global1) = node_count(node_global1) + 1 end do end do ! ! Average the running totals. ! dcdx(1:node_num) = dcdx(1:node_num) & / real ( node_count(1:node_num), kind = 8 ) dcdy(1:node_num) = dcdy(1:node_num) & / real ( node_count(1:node_num), kind = 8 ) return end subroutine div_q4 ( m, n, u, v, q, div, vort ) !*****************************************************************************80 ! !! DIV_Q4 estimates the divergence and vorticity of a discrete field. ! ! Discussion: ! ! The routine is given the values of a vector field ( U(X,Y), V(X,Y) ) at ! an array of points ( X(1:M), Y(1:N) ). ! ! The routine models the vector field over the interior of this region using ! a bilinear interpolant. It then uses the interpolant to estimate the ! value of the divergence: ! ! DIV(X,Y) = dU/dX + dV/dY ! ! and the vorticity: ! ! VORT(X,Y) = dV/dX - dU/dY ! ! at the center point of each of the bilinear elements. ! ! | | | ! (3,1)---(3,2)---(3,3)--- ! | | | ! | [2,1] | [2,2] | ! | | | ! (2,1)---(2,2)---(2,3)--- ! | | | ! | [1,1] | [1,2] | ! | | | ! (1,1)---(1,2)---(1,3)--- ! ! Here, the nodes labeled with parentheses represent the points at ! which the original (U,V) data is given, while the nodes labeled ! with square brackets represent the centers of the bilinear ! elements, where the approximations to the divergence and vorticity ! are made. ! ! The reason for evaluating the divergence and vorticity in this way ! is that the bilinear interpolant to the (U,V) data is not ! differentiable at the boundaries of the elements, nor especially at ! the nodes, but is an (infinitely differentiable) bilinear function ! in the interior of each element. If a value at the original nodes ! is strongly desired, then the average at the four surrounding ! central nodes may be taken. ! ! Reference Element Q4: ! ! | ! 1 4-----3 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 22 June 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) M, the number of data rows. M must be at least 2. ! ! Input, integer ( kind = 4 ) N, the number of data columns. N must be at least 2. ! ! Input, real ( kind = 8 ) U(M,N), V(M,N), the value of the components ! of a vector quantity whose divergence and vorticity are desired. ! A common example would be that U and V are the horizontal and ! vertical velocity component of a flow field. ! ! Input, real ( kind = 8 ) Q(2,4), the coordinates of the nodes of ! the quadrilateral, in counterclockwise order. ! ! Output, real ( kind = 8 ) DIV(M-1,N-1), an estimate for ! the divergence in the bilinear element that lies between ! data rows I and I+1, and data columns J and J+1. ! ! Output, real ( kind = 8 ) VORT(M-1,N-1), an estimate for ! the vorticity in the bilinear element that lies between ! data rows I and I+1, and data columns J and J+1. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) div(m-1,n-1) real ( kind = 8 ) dphidx(4) real ( kind = 8 ) dphidy(4) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ), parameter :: node_num = 1 real ( kind = 8 ) p(2) real ( kind = 8 ) phi(4) real ( kind = 8 ) q(2,4) real ( kind = 8 ) q2(2,4) real ( kind = 8 ) u(m,n) real ( kind = 8 ) v(m,n) real ( kind = 8 ) vort(m-1,n-1) if ( m <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_Q4 - Fatal error!' write ( *, '(a)' ) ' M must be at least 2,' write ( *, '(a,i8)' ) ' but the input value of M is ', m stop end if if ( n <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_Q4 - Fatal error!' write ( *, '(a)' ) ' N must be at least 2,' write ( *, '(a,i8)' ) ' but the input value of N is ', n stop end if do i = 1, m-1 q2(2,1) = ( real ( 2 * m - 2 * i, kind = 8 ) * q(2,1) & + real ( 2 * i - 2, kind = 8 ) * q(2,3) ) & / real ( 2 * m - 2, kind = 8 ) p(2) = ( real ( 2 * m - 2 * i - 1, kind = 8 ) * q(2,1) & + real ( 2 * i - 1, kind = 8 ) * q(2,3) ) & / real ( 2 * m - 2, kind = 8 ) q2(2,3) = ( real ( 2 * m - 2 * i - 2, kind = 8 ) * q(2,1) & + real ( 2 * i, kind = 8 ) * q(2,3) ) & / real ( 2 * m - 2, kind = 8 ) q2(2,2) = q2(2,1) q2(2,4) = q2(2,3) do j = 1, n-1 q2(1,1) = ( real ( 2 * n - 2 * j, kind = 8 ) * q(1,1) & + real ( 2 * j - 2, kind = 8 ) * q(1,3) ) & / real ( 2 * n - 2, kind = 8 ) p(1) = ( real ( 2 * n - 2 * j - 1, kind = 8 ) * q(1,1) & + real ( 2 * j - 1, kind = 8 ) * q(1,3) ) & / real ( 2 * n - 2, kind = 8 ) q2(1,3) = ( real ( 2 * n - 2 * j - 2, kind = 8 ) * q(1,1) & + real ( 2 * j, kind = 8 ) * q(1,3) ) & / real ( 2 * n - 2, kind = 8 ) q2(1,2) = q2(1,3) q2(1,4) = q2(1,1) call basis_mn_q4 ( q2, node_num, p, phi, dphidx, dphidy ) ! ! Note the following formula for the value of U and V at the same ! point that the divergence and vorticity are being evaluated. ! ! umid = u(i ,j ) * phi(1) & ! + u(i ,j+1) * phi(2) & ! + u(i+1,j+1) * phi(3) & ! + u(i+1,j ) * phi(4) ! ! vmid = v(i ,j ) * phi(1) & ! + v(i ,j+1) * phi(2) & ! + v(i+1,j+1) * phi(3) & ! + v(i+1,j ) * phi(4) ! div(i,j) = u(i ,j ) * dphidx(1) + v(i ,j ) * dphidy(1) & + u(i ,j+1) * dphidx(2) + v(i ,j+1) * dphidy(2) & + u(i+1,j+1) * dphidx(3) + v(i+1,j+1) * dphidy(3) & + u(i+1,j ) * dphidx(4) + v(i+1,j ) * dphidy(4) vort(i,j) = v(i ,j ) * dphidx(1) - u(i ,j ) * dphidy(1) & + v(i ,j+1) * dphidx(2) - u(i ,j+1) * dphidy(2) & + v(i+1,j+1) * dphidx(3) - u(i+1,j+1) * dphidy(3) & + v(i+1,j ) * dphidx(4) - u(i+1,j ) * dphidy(4) end do end do return end subroutine div_t3 ( m, n, u, v, q, div, vor ) !*****************************************************************************80 ! !! DIV_T3 estimates the divergence and vorticity of a discrete field. ! ! Discussion: ! ! The routine is given the values of a vector field ( U(X,Y), V(X,Y) ) at ! a regularly spaced grid of points ( X(1:M), Y(1:N) ). This grid is ! described implicitly by giving the values M, N, and the coordinates ! Q(2,4) of the bounding quadrilateral. (Note that Q need not be a ! rectangle.) ! ! The quadrilateral is suggested by the following diagram: ! ! ^ Q(1:2,4)-----Q(1:2,3) ! | | | ! N | | ! | | | ! V Q(1:2,1)-----Q(1:2,2) ! ! <--(M)---> ! ! The routine models the vector field over the interior of this region using ! a linear interpolant over 2*(M-1)*(N-1) triangles. It then uses the ! interpolant to estimate the value of the divergence: ! ! DIV(X,Y) = dU/dX + dV/dY ! ! and the vorticity: ! ! VOR(X,Y) = dV/dX - dU/dY ! ! at the centroid of each of the triangular elements. ! ! The grid is (somewhat arbitrarily) subdivided into triangular elements ! as suggested here: ! ! (3,1)---(3,2)---(3,3) ! | \ | \ | ! | \ | \ | ! | \ | \ | ! | \ | \ | ! (2,1)---(2,2)---(2,3) ! | \ | \ | ! | \ | \ | ! | \ | \ | ! | \ | \ | ! (1,1)---(1,2)---(1,3) ! ! In each triangular element, linear functions are used to interpolate ! the U and V data. The divergence and vorticity functions are then ! evaluated at the centroid of each element. ! ! This means that, given a grid of M X coordinates and N Y coordinates, ! we will construct 2 * ( M - 1 ) * ( N - 1 ) triangular elements. ! ! Modified: ! ! 16 November 2007 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, the number of data rows. M must be at least 2. ! ! Input, integer N, the number of data columns. N must be at least 2. ! ! Input, real ( kind = 8 ) U(M,N), V(M,N), the value of the components ! of a vector quantity whose divergence and vorticity are desired. ! A common example would be that U and V are the horizontal and ! vertical velocity component of a flow field. ! ! Input, real ( kind = 8 ) Q(2,4), the coordinates of the nodes of ! the quadrilateral, in counterclockwise order. ! ! Output, real ( kind = 8 ) DIV(2,M-1,N-1), an estimate for ! the divergence in the two linear elements that lie between ! data rows I and I+1, and data columns J and J+1. ! ! Output, real ( kind = 8 ) VOR(2,M-1,N-1), an estimate for ! the vorticity in the two linear elements that lie between ! data rows I and I+1, and data columns J and J+1. ! implicit none integer ( kind = 4 ) m integer ( kind = 4 ) n real ( kind = 8 ) div(2,m-1,n-1) real ( kind = 8 ) dphidx(3) real ( kind = 8 ) dphidy(3) real ( kind = 8 ) dudx real ( kind = 8 ) dudy real ( kind = 8 ) dvdx real ( kind = 8 ) dvdy integer ( kind = 4 ) i integer ( kind = 4 ) j real ( kind = 8 ) p(2) real ( kind = 8 ) phi(3) real ( kind = 8 ) q(2,4) real ( kind = 8 ) t(2,3) real ( kind = 8 ) u(m,n) real ( kind = 8 ) v(m,n) real ( kind = 8 ) vor(2,m-1,n-1) real ( kind = 8 ) xlb real ( kind = 8 ) xlt real ( kind = 8 ) xrb real ( kind = 8 ) xrt real ( kind = 8 ) xxlb real ( kind = 8 ) xxlt real ( kind = 8 ) xxrb real ( kind = 8 ) xxrt real ( kind = 8 ) ylb real ( kind = 8 ) ylt real ( kind = 8 ) yrb real ( kind = 8 ) yrt real ( kind = 8 ) yylb real ( kind = 8 ) yylt real ( kind = 8 ) yyrb real ( kind = 8 ) yyrt if ( m <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_T3 - Fatal error!' write ( *, '(a)' ) ' M must be at least 2,' write ( *, '(a,i8)' ) ' but the input value of M is ', m stop end if if ( n <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_T3 - Fatal error!' write ( *, '(a)' ) ' N must be at least 2,' write ( *, '(a,i8)' ) ' but the input value of N is ', n stop end if ! ! Consider the data between logical rows I and I + 1. ! do i = 1, m - 1 ! ! Consider the data between logical columns J and J + 1. ! do j = 1, n - 1 xlb = q(1,1) ylb = q(2,1) xrb = q(1,2) yrb = q(2,2) xrt = q(1,3) yrt = q(2,3) xlt = q(1,4) ylt = q(2,4) yylb = ( & real ( n - j, kind = 8 ) * ( real ( m - i, kind = 8 ) * ylb & + real ( i - 1, kind = 8 ) * yrb ) & / real ( m - 1, kind = 8 ) & + real ( j - 1, kind = 8 ) * ( real ( m - i, kind = 8 ) * ylt & + real ( i - 1, kind = 8 ) * yrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) yyrb = ( & real ( n - j, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * ylb & + real ( i, kind = 8 ) * yrb ) & / real ( m - 1, kind = 8 ) & + real ( j - 1, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * ylt & + real ( i, kind = 8 ) * yrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) yylt = ( & real ( n - j - 1, kind = 8 ) * ( real ( m - i, kind = 8 ) * ylb & + real ( i - 1, kind = 8 ) * yrb ) & / real ( m - 1, kind = 8 ) & + real ( j, kind = 8 ) * ( real ( m - i, kind = 8 ) * ylt & + real ( i - 1, kind = 8 ) * yrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) yyrt = ( & real ( n - j - 1, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * ylb & + real ( i, kind = 8 ) * yrb ) & / real ( m - 1, kind = 8 ) & + real ( j, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * ylt & + real ( i, kind = 8 ) * yrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) xxlb = ( & real ( n - j, kind = 8 ) * ( real ( m - i, kind = 8 ) * xlb & + real ( i - 1, kind = 8 ) * xrb ) & / real ( m - 1, kind = 8 ) & + real ( j - 1, kind = 8 ) * ( real ( m - i, kind = 8 ) * xlt & + real ( i - 1, kind = 8 ) * xrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) xxlt = ( & real ( n - j - 1, kind = 8 ) * ( real ( m - i, kind = 8 ) * xlb & + real ( i - 1, kind = 8 ) * xrb ) & / real ( m - 1, kind = 8 ) & + real ( j, kind = 8 ) * ( real ( m - i, kind = 8 ) * xlt & + real ( i - 1, kind = 8 ) * xrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) xxrb = ( & real ( n - j, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * xlb & + real ( i, kind = 8 ) * xrb ) & / real ( m - 1, kind = 8 ) & + real ( j - 1, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * xlt & + real ( i, kind = 8 ) * xrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) xxrt = ( & real ( n - j - 1, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * xlb & + real ( i, kind = 8 ) * xrb ) & / real ( m - 1, kind = 8 ) & + real ( j, kind = 8 ) * ( real ( m - i - 1, kind = 8 ) * xlt & + real ( i, kind = 8 ) * xrt ) & / real ( m - 1, kind = 8 ) & ) / real ( n - 1, kind = 8 ) ! write(*,'(i4,i4,8f8.3)')i,j,xxlb, yylb, xxrb, yyrb, xxrt,yyrt, xxlt,yylt ! ! (I,J+1) = LT-----RT = (I+1,J+1) ! |\ | ! | \ T2 | ! | \ | ! | \ | ! | T1 \ | ! | \ | ! (I,J) = LB-----RB = (I+1,J) ! t(1:2,1:3) = reshape ( (/ xxlb, yylb, xxrb, yyrb, xxrt, yyrt /), (/ 2, 3 /) ) p(1:2) = (/ xxlb + xxrb + xxrt, yylb + yyrb + yyrt /) / 3.0D+00 call basis_mn_t3 ( t, 1, p, phi, dphidx, dphidy ) dudx = u(i,j) * dphidx(1) + u(i+1,j) * dphidx(2) + u(i+1,j+1) * dphidx(3) dudy = u(i,j) * dphidy(1) + u(i+1,j) * dphidy(2) + u(i+1,j+1) * dphidy(3) dvdx = v(i,j) * dphidx(1) + v(i+1,j) * dphidx(2) + v(i+1,j+1) * dphidx(3) dvdy = v(i,j) * dphidy(1) + v(i+1,j) * dphidy(2) + v(i+1,j+1) * dphidy(3) div(1,i,j) = dudx + dvdy vor(1,i,j) = dvdx - dudy ! write ( *, '(4g14.6)' ) p(1), p(2), div(1,i,j), vor(1,i,j) t = reshape ( (/ xxrt, yyrt, xxlt, yylt, xxlb, yylb /), (/ 2, 3 /) ) p(1:2) = (/ xxrt + xxlt + xxlb, yyrt + yylt + yyrb /) / 3.0D+00 call basis_mn_t3 ( t, 1, p, phi, dphidx, dphidy ) dudx = u(i+1,j+1) * dphidx(1) + u(i,j+1) * dphidx(2) + u(i,j) * dphidx(3) dudy = u(i+1,j+1) * dphidy(1) + u(i,j+1) * dphidy(2) + u(i,j) * dphidy(3) dvdx = v(i+1,j+1) * dphidx(1) + v(i,j+1) * dphidx(2) + v(i,j) * dphidx(3) dvdy = v(i+1,j+1) * dphidy(1) + v(i,j+1) * dphidy(2) + v(i,j) * dphidy(3) div(2,i,j) = dudx + dvdy vor(2,i,j) = dvdx - dudy ! write ( *, '(4g14.6)' ) p(1), p(2), div(2,i,j), vor(2,i,j) end do end do return end function element_code ( i ) !*****************************************************************************80 ! !! ELEMENT_CODE returns the code for each element. ! ! List: ! ! I ELEMENT_CODE Definition ! - ------------ ---------- ! 1 Q4 4 node linear Lagrange/serendipity quadrilateral; ! 2 Q8 8 node quadratic serendipity quadrilateral; ! 3 Q9 9 node quadratic Lagrange quadrilateral; ! 4 Q12 12 node cubic serendipity quadrilateral; ! 5 Q16 16 node cubic Lagrange quadrilateral; ! 6 QL 6 node linear/quadratic quadrilateral; ! 7 T3 3 node linear triangle; ! 8 T4 4 node cubic bubble triangle ! 9 T6 6 node quadratic triangle; ! 10 T10 10 node cubic triangle. ! ! Modified: ! ! 03 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) I, the index of the element type. ! ! Output, character ( len = 3 ) ELEMENT_CODE, the code for the element type. ! implicit none character ( len = 3 ) element_code integer ( kind = 4 ) i if ( i == 1 ) then element_code = 'Q4' else if ( i == 2 ) then element_code = 'Q8' else if ( i == 3 ) then element_code = 'Q9' else if ( i == 4 ) then element_code = 'Q12' else if ( i == 5 ) then element_code = 'Q16' else if ( i == 6 ) then element_code = 'QL' else if ( i == 7 ) then element_code = 'T3' else if ( i == 8 ) then element_code = 'T4' else if ( i == 9 ) then element_code = 'T6' else if ( i == 10 ) then element_code = 'T10' else element_code = '???' end if return end subroutine elements_eps ( file_name, node_num, node_xy, code, & element_order, element_num, element_mask, element_node, title ) !*****************************************************************************80 ! !! ELEMENTS_EPS creates an EPS file image of the elements of a grid. ! ! Modified: ! ! 10 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) FILE_NAME, the name of the file to create. ! ! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes. ! ! Input, real ( kind = 8 ) NODE_XY(2,NODE_NUM), the ! coordinates of the nodes. ! ! Input, character ( len = * ) CODE, the code for the element. ! ! Input, integer ( kind = 4 ) ELEMENT_ORDER, the element order. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, logical ELEMENT_MASK(ELEMENT_NUM), a mask for the elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), the ! element->node data. ! ! Input, character ( len = * ) TITLE, a title for the plot. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) element_order integer ( kind = 4 ) node_num real ( kind = 8 ) ave_x real ( kind = 8 ) ave_y integer ( kind = 4 ), parameter :: circle_size = 3 character ( len = * ) code real ( kind = 8 ) dif integer ( kind = 4 ) element logical element_mask(element_num) integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) eps_unit integer ( kind = 4 ) eps_x integer ( kind = 4 ) eps_y character ( len = * ) file_name integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) j integer ( kind = 4 ) local integer ( kind = 4 ) next_boundary_node integer ( kind = 4 ) node logical node_mask(node_num) real ( kind = 8 ) node_x_max real ( kind = 8 ) node_x_min real ( kind = 8 ) node_xy(2,node_num) real ( kind = 8 ) node_y_max real ( kind = 8 ) node_y_min real ( kind = 8 ) r8_huge real ( kind = 8 ) scale character ( len = 40 ) string character ( len = * ) title ! ! Determine the range of the unmasked elements. ! node_x_min = r8_huge ( node_x_min ) node_x_max = -r8_huge ( node_x_max ) node_y_min = r8_huge ( node_y_min ) node_y_max = -r8_huge ( node_y_max ) node_mask(1:node_num) = .false. do element = 1, element_num if ( element_mask(element) ) then do j = 1, element_order node = element_node(j,element) node_mask(node) = .true. node_x_min = min ( node_x_min, node_xy(1,node) ) node_x_max = max ( node_x_max, node_xy(1,node) ) node_y_min = min ( node_y_min, node_xy(2,node) ) node_y_max = max ( node_y_max, node_xy(2,node) ) end do end if end do if ( node_y_max - node_y_min < node_x_max - node_x_min ) then scale = node_x_max - node_x_min dif = ( node_x_max - node_x_min ) - ( node_y_max - node_y_min ) node_y_max = node_y_max + 0.5D+00 * dif node_y_min = node_y_min - 0.5D+00 * dif else scale = node_y_max - node_y_min dif = ( node_y_max - node_y_min ) - ( node_x_max - node_x_min ) node_x_max = node_x_max + 0.5D+00 * dif node_x_min = node_x_min - 0.5D+00 * dif end if call get_unit ( eps_unit ) open ( unit = eps_unit, file = file_name, status = 'replace', iostat = ios ) if ( ios /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'ELEMENTS_EPS - Fatal error!' write ( *, '(a)' ) ' Could not open the output EPS file.' stop end if call timestring ( string ) write ( eps_unit, '(a)' ) '%!PS-Adobe-3.0 EPSF-3.0' write ( eps_unit, '(a)' ) '%%Creator: elements_eps(fempack.f90)' write ( eps_unit, '(a)' ) '%%Title: ' // trim ( file_name ) write ( eps_unit, '(a)' ) '%%CreationDate: ' // trim ( string ) write ( eps_unit, '(a)' ) '%%Pages: 1' write ( eps_unit, '(a)' ) '%%BoundingBox: 36 36 576 756' write ( eps_unit, '(a)' ) '%%Document-Fonts: Times-Roman' write ( eps_unit, '(a)' ) '%%LanguageLevel: 1' write ( eps_unit, '(a)' ) '%%EndComments' write ( eps_unit, '(a)' ) '%%BeginProlog' write ( eps_unit, '(a)' ) '/inch {72 mul} def' write ( eps_unit, '(a)' ) '%%EndProlog' write ( eps_unit, '(a)' ) '%%Page: 1 1' write ( eps_unit, '(a)' ) 'save' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Set RGB line color.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.9000 0.9000 0.9000 setrgbcolor' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw a gray border around the page.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) 'newpath' write ( eps_unit, '(a)' ) ' 36 126 moveto' write ( eps_unit, '(a)' ) ' 576 126 lineto' write ( eps_unit, '(a)' ) ' 576 666 lineto' write ( eps_unit, '(a)' ) ' 36 666 lineto' write ( eps_unit, '(a)' ) ' 36 126 lineto' write ( eps_unit, '(a)' ) 'stroke' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Set RGB line color.' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the plot:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) '/Times-Roman findfont 0.50 inch scalefont setfont' write ( eps_unit, '(a)' ) ' 36 666 moveto' write ( eps_unit, '(a)' ) '(' // trim ( title ) // ') show' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Define a clipping polygon' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 36 126 moveto' write ( eps_unit, '(a)' ) ' 576 126 lineto' write ( eps_unit, '(a)' ) ' 576 666 lineto' write ( eps_unit, '(a)' ) ' 36 666 lineto' write ( eps_unit, '(a)' ) ' 36 126 lineto' write ( eps_unit, '(a)' ) 'clip newpath' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw filled dots at each node:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 0.9000 setrgbcolor' do node = 1, node_num if ( node_mask(node) ) then eps_x = int & ( ( node_x_max - node_xy(1,node) ) * 61.0D+00 & + ( + node_xy(1,node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_xy(2,node) ) * 151.0D+00 & + ( node_xy(2,node) - node_y_min ) * 641.0D+00 ) & / scale write ( eps_unit, '(a,i4,2x,i4,2x,i4,a)' ) & 'newpath ', eps_x, eps_y, circle_size, ' 0 360 arc closepath fill' end if end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the nodes:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.0000 0.0000 1.0000 setrgbcolor' write ( eps_unit, '(a)' ) '/Times-Roman findfont 0.20 inch scalefont setfont' do node = 1, node_num if ( node_mask(node) ) then eps_x = int & ( ( node_x_max - node_xy(1,node) ) * 61.0D+00 & + ( + node_xy(1,node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_xy(2,node) ) * 151.0D+00 & + ( node_xy(2,node) - node_y_min ) * 641.0D+00 ) & / scale write ( string, '(i4)' ) node string = adjustl ( string ) write ( eps_unit, '(i4,2x,i4,a)' ) eps_x, eps_y+5, & ' moveto (' // trim ( string ) // ') show' end if end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Draw the element sides:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 0.9000 0.0000 0.0000 setrgbcolor' do element = 1, element_num if ( .not. element_mask(element) ) then cycle end if local = 1 node = element_node(local,element) eps_x = int & ( ( node_x_max - node_xy(1,node) ) * 61.0D+00 & + ( + node_xy(1,node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_xy(2,node) ) * 151.0D+00 & + ( node_xy(2,node) - node_y_min ) * 641.0D+00 ) & / scale write ( eps_unit, '(a,i4,2x,i4,a)' ) 'newpath ', eps_x, eps_y, ' moveto' do local = next_boundary_node ( local, code ) node = element_node(local,element) eps_x = int & ( ( node_x_max - node_xy(1,node) ) * 61.0D+00 & + ( + node_xy(1,node) - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - node_xy(2,node) ) * 151.0D+00 & + ( node_xy(2,node) - node_y_min ) * 641.0D+00 ) & / scale write ( eps_unit, '(i4,2x,i4,a)' ) eps_x, eps_y, ' lineto' if ( local == 1 ) then exit end if end do write ( eps_unit, '(a)' ) 'stroke' end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% Label the elements:' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) ' 1.0000 0.0000 0.0000 setrgbcolor' write ( eps_unit, '(a)' ) '/Times-Roman findfont 0.30 inch scalefont setfont' do element = 1, element_num if ( .not. element_mask(element) ) then cycle end if 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 = 8 ) ave_y = ave_y / real ( element_order, kind = 8 ) eps_x = int & ( ( node_x_max - ave_x ) * 61.0D+00 & + ( + ave_x - node_x_min ) * 551.0D+00 ) & / scale eps_y = int & ( ( node_y_max - ave_y ) * 151.0D+00 & + ( ave_y - node_y_min ) * 641.0D+00 ) & / scale write ( string, '(i4)' ) element string = adjustl ( string ) write ( eps_unit, '(i4,2x,i4,a)' ) eps_x, eps_y, ' moveto (' & // trim ( string ) // ') show' end do write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) 'restore showpage' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '% End of page' write ( eps_unit, '(a)' ) '%' write ( eps_unit, '(a)' ) '%%Trailer' write ( eps_unit, '(a)' ) '%%EOF' close ( unit = eps_unit ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' ELEMENTS_EPS: An encapsulated PostScript file' write ( *, '(a)' ) ' was created containing an image of the nodes and' write ( *, '(a)' ) ' elements. The file is named "' & // trim ( file_name ) // '".' 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. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer ( kind = 4 ) IUNIT, the free unit number. ! implicit none integer ( kind = 4 ) i integer ( kind = 4 ) ios integer ( kind = 4 ) iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine grid_element ( code, element_order, nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_ELEMENT returns the element grid associated with any available element. ! ! Modified: ! ! 03 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element desired. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3', ! 'T4', 'T6' and 'T10'. ! ! Input, integer ( kind = 4 ) ELEMENT_ORDER, the number of nodes ! per element. ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of quadrilaterals ! along the X and Y directions. The number of elements generated will be ! NELEMX * NELEMY for quadrilaterals, or 2 * NELEMX * NELEMY for ! triangles. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), ! the nodes that form each element. ! implicit none integer ( kind = 4 ) element_order character ( len = * ) code integer ( kind = 4 ) element_node(element_order,*) integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy logical s_eqi if ( s_eqi ( code, 'Q4' ) ) then call grid_q4_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'Q8' ) ) then call grid_q8_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'Q9' ) ) then call grid_q9_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'Q12' ) ) then call grid_q12_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'Q16' ) ) then call grid_q16_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'QL' ) ) then call grid_ql_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'T3' ) ) then call grid_t3_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'T4' ) ) then call grid_t4_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'T6' ) ) then call grid_t6_element ( nelemx, nelemy, element_node ) else if ( s_eqi ( code, 'T10' ) ) then call grid_t10_element ( nelemx, nelemy, element_node ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_ELEMENT - Fatal error!' write ( *, '(a)' ) ' Illegal value of CODE = "' // trim ( code ) // '".' stop end if return end subroutine grid_element_num ( code, nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_ELEMENT_NUM returns the number of elements in a grid. ! ! Discussion: ! ! The number of elements generated will be NELEMX * NELEMY for ! quadrilaterals, or 2 * NELEMX * NELEMY for triangles. ! ! Modified: ! ! 03 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element desired. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3', ! 'T4', 'T6' and 'T10'. ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none character ( len = * ) code integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy logical s_eqi if ( s_eqi ( code, 'Q4' ) ) then call grid_q4_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'Q8' ) ) then call grid_q8_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'Q9' ) ) then call grid_q9_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'Q12' ) ) then call grid_q12_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'Q16' ) ) then call grid_q16_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'QL' ) ) then call grid_ql_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'T3' ) ) then call grid_t3_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'T4' ) ) then call grid_t4_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'T6' ) ) then call grid_t6_element_num ( nelemx, nelemy, element_num ) else if ( s_eqi ( code, 'T10' ) ) then call grid_t10_element_num ( nelemx, nelemy, element_num ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_ELEMENT_NUM - Fatal error!' write ( *, '(a)' ) ' Illegal value of CODE = "' // trim ( code ) // '".' element_num = -1 stop end if return end subroutine grid_node_num ( code, nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_NODE_NUM returns the number of nodes in a grid. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element desired. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3', ! 'T4', 'T6' and 'T10'. ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of elements in the grid. ! implicit none character ( len = * ) code integer ( kind = 4 ) node_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy logical s_eqi if ( s_eqi ( code, 'Q4' ) ) then call grid_q4_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'Q8' ) ) then call grid_q8_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'Q9' ) ) then call grid_q9_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'Q12' ) ) then call grid_q12_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'Q16' ) ) then call grid_q16_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'QL' ) ) then call grid_ql_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'T3' ) ) then call grid_t3_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'T4' ) ) then call grid_t4_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'T6' ) ) then call grid_t6_node_num ( nelemx, nelemy, node_num ) else if ( s_eqi ( code, 'T10' ) ) then call grid_t10_node_num ( nelemx, nelemy, node_num ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_NODE_NUM - Fatal error!' write ( *, '(a)' ) ' Illegal value of CODE = "' // trim ( code ) // '".' node_num = -1 stop end if return end subroutine grid_nodes_01 ( x_num, y_num, node_xy ) !*****************************************************************************80 ! !! GRID_NODES_01 returns an equally spaced rectangular grid of nodes in the unit square. ! ! Example: ! ! X_NUM = 5 ! Y_NUM = 3 ! ! NODE_XY = ! ( 0, 0.25, 0.5, 0.75, 1, 0, 0.25, 0.5, 0.75, 1, 0, 0.25, 0.5, 0.75, 1; ! 0, 0, 0, 0, 0, 0.5, 0.5, 0.5, 0.5, 0.5, 1, 1.0, 1.0, 1.0, 1 ) ! ! Modified: ! ! 14 May 2008 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) X_NUM, Y_NUM, the number of nodes in the ! X and Y directions. ! ! Output, real ( kind = 8 ) NODE_XY(2,X_NUM*Y_NUM), the coordinates of ! the nodes. ! implicit none integer ( kind = 4 ) x_num integer ( kind = 4 ) y_num integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) node_num real ( kind = 8 ) node_xy(2,x_num*y_num) node_num = x_num * y_num node_xy(1:2,1:node_num) = 0.0D+00 if ( x_num == 1 ) then node_xy(1,1:node_num) = 0.5D+00 else do i = 1, x_num node_xy(1,i:i+(y_num-1)*x_num:x_num) = real ( i - 1, kind = 8 ) & / real ( x_num - 1, kind = 8 ) end do end if if ( y_num == 1 ) then node_xy(2,1:node_num) = 0.5D+00 else do j = 1, y_num node_xy(2,1+(j-1)*x_num:j*x_num) = real ( j - 1, kind = 8 ) & / real ( y_num - 1, kind = 8 ) end do end if return end subroutine grid_print ( element_order, element_num, element_node ) !*****************************************************************************80 ! !! GRID_PRINT prints the elements that form a grid. ! ! Modified: ! ! 24 March 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) ELEMENT_ORDER, the number of nodes ! per element. ! ! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements. ! ! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM), ! the nodes that form each element. ! implicit none integer ( kind = 4 ) element_order integer ( kind = 4 ) element_num integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,element_num) integer ( kind = 4 ) i write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' GRID_PRINT: Element -> Node table.' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Number of elements = ', element_num write ( *, '(a,i8)' ) ' Element order = ', element_order write ( *, '(a)' ) ' ' write ( *, '(2x,a,3x,20i3)' ) ' #', ( i, i = 1, element_order ) write ( *, '(a)' ) ' ' do element = 1, element_num write ( *, '(2x,i3,3x,20i3)' ) & element, element_node(1:element_order,element) end do return end subroutine grid_q4_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_Q4_ELEMENT produces an element grid of 4 node quadrilaterals. ! ! Discussion: ! ! For each element, the nodes are listed in counter-clockwise order. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 6, 5; ! 2, 3, 7, 6; ! 3, 4, 8, 7; ! 5, 6, 10, 9; ! 6, 7, 11, 10; ! 7, 8, 12, 11. ! ! Grid: ! ! 9---10---11---12 ! | | | | ! | | | | ! | 4 | 5 | 6 | ! | | | | ! 5----6----7----8 ! | | | | ! | | | | ! | 1 | 2 | 3 | ! | | | | ! 1----2----3----4 ! ! Reference Element Q4: ! ! | ! 1 4-----3 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 07 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(4,NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 4 integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) ne integer ( kind = 4 ) nw integer ( kind = 4 ) se integer ( kind = 4 ) sw ! ! Node labeling: ! ! NW---NE ! | | ! SW---SE ! element = 0 do j = 1, nelemy do i = 1, nelemx sw = i + ( j - 1 ) * ( nelemx + 1 ) se = i + 1 + ( j - 1 ) * ( nelemx + 1 ) nw = i + j * ( nelemx + 1 ) ne = i + 1 + j * ( nelemx + 1 ) element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = ne element_node(4,element) = nw end do end do return end subroutine grid_q4_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_Q4_ELEMENT_NUM counts the elements in a grid of 4 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = NELEMX * NELEMY = 6 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = nelemx * nelemy return end subroutine grid_q4_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_Q4_NODE_NUM counts the nodes in a grid of 4 node quadrilaterals. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = ( nelemx + 1 ) * ( nelemy + 1 ) return end subroutine grid_q8_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_Q8_ELEMENT produces an element grid of 8 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 3, 14, 12, 2, 9, 13, 8; ! 3, 5, 16, 14, 4, 10, 15, 9; ! 5, 7, 18, 16, 6, 11, 17, 10; ! 12, 14, 25, 23, 13, 20, 24, 19; ! 14, 16, 27, 25, 15, 21, 26, 20; ! 16, 18, 29, 27, 17, 22, 28, 21. ! ! Diagram: ! ! 23---24---25---26---27---28---29 ! | | | | ! | | | | ! 19 20 21 22 ! | | | | ! | 4 | 5 | 6 | ! 12---13---14---15---16---17---18 ! | | | | ! | | | | ! 8 9 10 11 ! | | | | ! | 1 | 2 | 3 | ! 1----2----3----4----5----6----7 ! ! Reference Element Q8: ! ! | ! 1 4--7--3 ! | | | ! | | | ! S 8 6 ! | | | ! | | | ! 0 1--5--2 ! | ! +--0--R--1--> ! ! Modified: ! ! 08 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(8,NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 8 integer ( kind = 4 ) e integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) n integer ( kind = 4 ) ne integer ( kind = 4 ) nw integer ( kind = 4 ) s integer ( kind = 4 ) se integer ( kind = 4 ) sw integer ( kind = 4 ) w ! ! Node labeling: ! ! NW----N----NE ! | | ! W (C) E ! | | ! SW----S----SE ! element = 0 do j = 1, nelemy do i = 1, nelemx sw = ( j - 1 ) * ( 3 * nelemx + 2 ) + 2 * i - 1 w = sw + 2 * nelemx + 2 - i nw = sw + 3 * nelemx + 2 s = sw + 1 n = sw + ( 3 * nelemx + 2 ) + 1 se = sw + 2 e = sw + 2 * nelemx + 2 - i + 1 ne = sw + ( 3 * nelemx + 2 ) + 2 element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = ne element_node(4,element) = nw element_node(5,element) = s element_node(6,element) = e element_node(7,element) = n element_node(8,element) = w end do end do return end subroutine grid_q8_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_Q8_ELEMENT_NUM counts the elements in a grid of 8 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = NELEMX * NELEMY = 6 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = nelemx * nelemy return end subroutine grid_q8_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_Q8_NODE_NUM counts the nodes in a grid of 8 node quadrilaterals. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of node in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = 3 * nelemx * nelemy + 2 * nelemx + 2 * nelemy + 1 return end subroutine grid_q9_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_Q9_ELEMENT produces an element grid of 9 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 3, 17, 15, 2, 10, 16, 8, 9; ! 3, 5, 19, 17, 4, 12, 18, 10, 11; ! 5, 7, 21, 19, 6, 14, 20, 12, 13; ! 15, 17, 31, 29, 16, 24, 30, 22, 23; ! 17, 19, 33, 31, 18, 26, 32, 24, 25; ! 19, 21, 35, 33, 20, 28, 34, 26, 27. ! ! Grid: ! ! 29---30---31---32---33---34---35 ! | . | . | . | ! | . | . | . | ! 22 . 23 . 24 . 25 . 26 . 27 . 28 ! | . | . | . | ! | 4 . | 5 . | 6 . | ! 15---16---17---18---19---20---21 ! | . | . | . | ! | . | . | . | ! 8 . 9 . 10 . 11 . 12 . 13 . 14 ! | . | . | . | ! | 1 . | 2 . | 3 . | ! 1----2----3----4----5----6----7 ! ! Reference Element Q9: ! ! | ! 1 4--7--3 ! | | | ! | | | ! S 8 9 6 ! | | | ! | | | ! 0 1--5--2 ! | ! +--0--R--1--> ! ! Modified: ! ! 09 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(9,NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 9 integer ( kind = 4 ) c integer ( kind = 4 ) e integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) n integer ( kind = 4 ) ne integer ( kind = 4 ) nw integer ( kind = 4 ) s integer ( kind = 4 ) se integer ( kind = 4 ) sw integer ( kind = 4 ) w ! ! Node labeling: ! ! NW----N----NE ! | | ! W C E ! | | ! SW----S----SE ! element = 0 do j = 1, nelemy do i = 1, nelemx sw = 2 * ( j - 1 ) * ( 2 * nelemx + 1 ) + 2 * ( i - 1 ) + 1 w = sw + 2 * nelemx + 1 nw = sw + 2 * ( 2 * nelemx + 1 ) s = sw + 1 c = sw + 1 + 2 * nelemx + 1 n = sw + 1 + 2 * ( 2 * nelemx + 1 ) se = sw + 2 e = sw + 2 + 2 * nelemx + 1 ne = sw + 2 + 2 * ( 2 * nelemx + 1 ) element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = ne element_node(4,element) = nw element_node(5,element) = s element_node(6,element) = e element_node(7,element) = n element_node(8,element) = w element_node(9,element) = c end do end do return end subroutine grid_q9_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_Q9_ELEMENT_NUM counts the elements in a grid of 9 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = NELEMX * NELEMY = 6 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = nelemx * nelemy return end subroutine grid_q9_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_Q9_NODE_NUM counts the nodes in a grid of 9 node quadrilaterals. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = ( 2 * nelemx + 1 ) * ( 2 * nelemy + 1 ) return end subroutine grid_q12_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_Q12_ELEMENT produces an element grid of 12 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 3, 4, 11, 12, 15, 16, 19, 20, 21, 22; ! 4, 5, 6, 7, 12, 13, 16, 17, 22, 23, 24, 25; ! 7, 8, 9, 10, 13, 14, 17, 18, 25, 26, 27, 28; ! 19, 20, 21, 22, 29, 30, 33, 34, 37, 38, 39, 40; ! 22, 23, 24, 25, 30, 31, 34, 35, 40, 41, 42, 43; ! 25, 26, 27, 28, 31, 32, 35, 36, 43, 44, 45, 46. ! ! Grid: ! ! 37-38-39-40-41-42-43-44-45-46 ! | | | | ! 33 34 35 36 ! | | | | ! 29 30 31 32 ! | 4 | 5 | 6 | ! 19-20-21-22-23-24-25-26-27-28 ! | | | | ! 15 16 17 18 ! | | | | ! 11 12 13 14 ! | 1 | 2 | 3 | ! 1--2--3--4--5--6--7--8--9-10 ! ! Reference Element Q12: ! ! | ! 1 9-10-11-12 ! | | | ! | 7 8 ! S | | ! | 5 6 ! | | | ! 0 1--2--3--4 ! | ! +--0---R---1--> ! ! Modified: ! ! 07 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(12,NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 12 integer ( kind = 4 ) base integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j element = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * ( 5 * nelemx + 3 ) + 1 element = element + 1 element_node( 1,element) = base + ( i - 1 ) * 3 element_node( 2,element) = base + ( i - 1 ) * 3 + 1 element_node( 3,element) = base + ( i - 1 ) * 3 + 2 element_node( 4,element) = base + ( i - 1 ) * 3 + 3 element_node( 5,element) = base + 3 * nelemx + i element_node( 6,element) = base + 3 * nelemx + i + 1 element_node( 7,element) = base + 4 * nelemx + i + 1 element_node( 8,element) = base + 4 * nelemx + i + 2 element_node( 9,element) = base + 5 * nelemx + 3 * i element_node(10,element) = base + 5 * nelemx + 3 * i + 1 element_node(11,element) = base + 5 * nelemx + 3 * i + 2 element_node(12,element) = base + 5 * nelemx + 3 * i + 3 end do end do return end subroutine grid_q12_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_Q12_ELEMENT_NUM counts the elements in a grid of 12 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = NELEMX * NELEMY = 6 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = nelemx * nelemy return end subroutine grid_q12_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_Q12_NODE_NUM counts the nodes in a grid of 12 node quadrilaterals. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of node in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = 5 * nelemx * nelemy + 3 * nelemx + 3 * nelemy + 1 return end subroutine grid_q16_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_Q16_ELEMENT produces an element grid of 16 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 2, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 3, 4, 8, 9, 10, 11, 15, 16, 17, 18, 22, 23, 24, 25; ! 4, 5, 6, 7, 11, 12, 13, 14, 18, 19, 20, 21, 25, 26, 27, 28; ! 22, 23, 24, 25, 29, 30, 31, 32, 36, 37, 38, 39, 43, 44, 45, 46; ! 25, 26, 27, 28, 32, 33, 34, 35, 39, 40, 41, 42, 46, 47, 48, 49. ! ! Grid: ! ! 43-44-45-46-47-48-49 ! | | | ! | | | ! 36 37 38 39 40 41 42 ! | | | ! | | | ! 29 30 31 32 33 34 35 ! | | | ! | 3 | 4 | ! 22-23-24-25-26-27-28 ! | | | ! | | | ! 15 16 17 18 19 20 21 ! | | | ! | | | ! 8 9 10 11 12 13 14 ! | | | ! | 1 | 2 | ! 1--2--3--4--5--6--7 ! ! Reference Element Q16: ! ! | ! 1 13--14--15--16 ! | | : : | ! | | : : | ! | 9..10..11..12 ! S | : : | ! | | : : | ! | 5...6...7...8 ! | | : : | ! | | : : | ! 0 1---2---3---4 ! | ! +--0-----R-----1--> ! ! Modified: ! ! 08 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(16,NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 16 integer ( kind = 4 ) base integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j element = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2 element = element + 1 element_node( 1,element) = base element_node( 2,element) = base + 1 element_node( 3,element) = base + 2 element_node( 4,element) = base + 3 element_node( 5,element) = base + ( 3 * nelemx + 1 ) element_node( 6,element) = base + ( 3 * nelemx + 1 ) + 1 element_node( 7,element) = base + ( 3 * nelemx + 1 ) + 2 element_node( 8,element) = base + ( 3 * nelemx + 1 ) + 3 element_node( 9,element) = base + 2 * ( 3 * nelemx + 1 ) element_node(10,element) = base + 2 * ( 3 * nelemx + 1 ) + 1 element_node(11,element) = base + 2 * ( 3 * nelemx + 1 ) + 2 element_node(12,element) = base + 2 * ( 3 * nelemx + 1 ) + 3 element_node(13,element) = base + 3 * ( 3 * nelemx + 1 ) element_node(14,element) = base + 3 * ( 3 * nelemx + 1 ) + 1 element_node(15,element) = base + 3 * ( 3 * nelemx + 1 ) + 2 element_node(16,element) = base + 3 * ( 3 * nelemx + 1 ) + 3 end do end do return end subroutine grid_q16_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_Q16_ELEMENT_NUM counts the elements in a grid of 16 node quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = NELEMX * NELEMY = 6 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = nelemx * nelemy return end subroutine grid_q16_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_Q16_NODE_NUM counts the nodes in a grid of 16 node quadrilaterals. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = ( 3 * nelemx + 1 ) * ( 3 * nelemy + 1 ) return end subroutine grid_ql_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_QL_ELEMENT produces an element grid of 6 node quadratics/linears. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 3, 8, 9, 10; ! 3, 4, 5, 10, 11, 12; ! 5, 6, 7, 12, 13, 14; ! 8, 9, 10, 15, 16, 17; ! 10, 11, 12, 17, 18, 19; ! 12, 13, 14, 19, 20, 21. ! ! Grid: ! ! 15---16---17---18---19---20---21 ! | | | | ! | | | | ! | 4 | 5 | 6 | ! | | | | ! | | | | ! 8----9---10---11---12---13---14 ! | | | | ! | | | | ! | 1 | 2 | 3 | ! | | | | ! | | | | ! 1----2----3----4----5----6----7 ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. X will the the "quadratic direction", and ! Y will be the "linear direction". ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(6,NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 6 integer ( kind = 4 ) base integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j element = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * ( 2 * nelemx + 1 ) + 2 * i - 1 element = element + 1 element_node(1,element) = base element_node(2,element) = base + 1 element_node(3,element) = base + 2 element_node(4,element) = base + ( 2 * nelemx + 1 ) element_node(5,element) = base + ( 2 * nelemx + 1 ) + 1 element_node(6,element) = base + ( 2 * nelemx + 1 ) + 2 end do end do return end subroutine grid_ql_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_QL_ELEMENT_NUM counts the elements in a grid of QL quadrilaterals. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = NELEMX * NELEMY = 6 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = nelemx * nelemy return end subroutine grid_ql_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_QL_NODE_NUM counts the nodes in a grid of QL quadrilaterals. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = 2 * nelemx * nelemy + 2 * nelemx + nelemy + 1 return end subroutine grid_shape_2d ( n, a, n1, n2 ) !*****************************************************************************80 ! !! GRID_SHAPE_2D guesses the shape N1 by N2 of a vector of data. ! ! Discussion: ! ! The data vector A is assumed to contain N1 * N2 values, with ! where each of N2 values is repeated N1 times. ! ! Example: ! ! Input: ! ! A = ( 2, 2, 2, 7, 7, 7 ) ! ! Output: ! ! N1 = 3, N2 = 2 ! ! Modified: ! ! 19 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) N, the number of data values. ! ! Input, real ( kind = 8 ) A(N), the data, which should have the properties ! described above. ! ! Output, integer ( kind = 4 ) N1, N2, the "shape" of the data in the array. ! implicit none integer ( kind = 4 ) n real ( kind = 8 ) a(n) integer ( kind = 4 ) i integer ( kind = 4 ) n1 integer ( kind = 4 ) n2 ! ! Make a guess for N1. ! i = 1 n1 = 1 do i = 2, n if ( a(i) /= a(1) ) then exit end if n1 = n1 + 1 end do ! ! Guess that N2 = N / N1. ! n2 = n / n1 return end subroutine grid_t3_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_T3_ELEMENT produces an element grid of pairs of 3 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 5; ! 6, 5, 2; ! 2, 3, 6; ! 7, 6, 3; ! 3, 4, 7; ! 8, 7, 4; ! 5, 6, 9; ! 10, 9, 6; ! 6, 7, 10; ! 11, 10, 7; ! 7, 8, 11; ! 12, 11, 8. ! ! Grid: ! ! 9---10---11---12 ! |\ 8 |\10 |\12 | ! | \ | \ | \ | ! | \ | \ | \ | ! | 7\| 9\| 11\| ! 5----6----7----8 ! |\ 2 |\ 4 |\ 6 | ! | \ | \ | \ | ! | \ | \ | \ | ! | 1\| 3\| 5\| ! 1----2----3----4 ! ! Reference Element T3: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 07 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(3,2*NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 3 integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,2*nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) ne integer ( kind = 4 ) nw integer ( kind = 4 ) se integer ( kind = 4 ) sw ! ! Node labeling: ! ! NW--NE ! |\ | ! | \| ! SW--SE ! element = 0 do j = 1, nelemy do i = 1, nelemx sw = i + ( j - 1 ) * ( nelemx + 1 ) se = i + 1 + ( j - 1 ) * ( nelemx + 1 ) nw = i + j * ( nelemx + 1 ) ne = i + 1 + j * ( nelemx + 1 ) element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = nw element = element + 1 element_node(1,element) = ne element_node(2,element) = nw element_node(3,element) = se end do end do return end subroutine grid_t3_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_T3_ELEMENT_NUM counts the elements in a grid of 3 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = 2 * NELEMX * NELEMY = 12 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = 2 * nelemx * nelemy return end subroutine grid_t3_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_T3_NODE_NUM counts the nodes in a grid of 3 node triangles. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = ( nelemx + 1 ) * ( nelemy + 1 ) return end subroutine grid_t4_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_T4_ELEMENT produces an element grid of pairs of 4 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 11, 5; ! 12, 11, 2, 8; ! 2, 3, 12, 6; ! 13, 12, 3, 9; ! 3 4 13, 7; ! 14, 13, 4, 10; ! 11, 12, 21, 15; ! 22, 21, 12, 18; ! 12, 13, 22, 16; ! 23, 22, 13, 19; ! 13 14 23, 17; ! 24, 23, 14, 20; ! ! Grid: ! ! 21---22---23---24 ! |\18 |\19 |\20 | ! | \ | \ | \ | ! | \ | \ | \ | ! | 15\| 16\| 17\| ! 11---12---13---14 ! |\ 8 |\ 9 |\10 | ! | \ | \ | \ | ! | \ | \ | \ | ! | 5 \| 6\| 7\| ! 1----2----3----4 ! ! Reference Element T4: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | 4 \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(4,2*NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 4 integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,2*nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) nc integer ( kind = 4 ) ne integer ( kind = 4 ) nw integer ( kind = 4 ) sc integer ( kind = 4 ) se integer ( kind = 4 ) sw ! ! Node labeling: ! ! NW----NE ! |\ | ! | \NC| ! |SC\ | ! | \| ! SW---SE ! element = 0 do j = 1, nelemy do i = 1, nelemx sw = i + ( j - 1 ) * ( 3 * nelemx + 1 ) se = sw + 1 sc = sw + nelemx + 1 nc = sw + 2 * nelemx + 1 nw = sw + 3 * nelemx + 1 ne = sw + 3 * nelemx + 2 element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = nw element_node(4,element) = sc element = element + 1 element_node(1,element) = ne element_node(2,element) = nw element_node(3,element) = se element_node(4,element) = nc end do end do return end subroutine grid_t4_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_T4_ELEMENT_NUM counts the elements in a grid of 4 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = 2 * NELEMX * NELEMY = 12 ! ! Modified: ! ! 03 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = 2 * nelemx * nelemy return end subroutine grid_t4_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_T4_NODE_NUM counts the nodes in a grid of 4 node triangles. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = ( nelemx + 1 ) * ( nelemy + 1 ) + 2 * nelemx * nelemy return end subroutine grid_t6_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_T6_ELEMENT produces an element grid of pairs of 6 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 3, 15, 2, 9, 8; ! 17, 15, 3, 16, 9, 10; ! 3, 5, 17, 4, 11, 10; ! 19, 17, 5, 18, 11, 12; ! 5, 7, 19, 6, 13, 12; ! 21, 19, 7, 20, 13, 14; ! 15, 17, 29, 16, 23, 22; ! 31, 29, 17, 30, 23, 24; ! 17, 19, 31, 18, 25, 24; ! 33, 31, 19, 32, 25, 26; ! 19, 21, 33, 20, 27, 26; ! 35, 33, 21, 34, 27, 28. ! ! Grid: ! ! 29-30-31-32-33-34-35 ! |\ 8 |\10 |\12 | ! | \ | \ | \ | ! 22 23 24 25 26 27 28 ! | \ | \ | \ | ! | 7 \| 9 \| 11 \| ! 15-16-17-18-19-20-21 ! |\ 2 |\ 4 |\ 6 | ! | \ | \ | \ | ! 8 9 10 11 12 13 14 ! | \ | \ | \ | ! | 1 \| 3 \| 5 \| ! 1--2--3--4--5--6--7 ! ! Reference Element T6: ! ! | ! 1 3 ! | |\ ! | | \ ! S 6 5 ! | | \ ! | | \ ! 0 1--4--2 ! | ! +--0--R--1--> ! ! Modified: ! ! 10 April 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(6,2*NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 6 integer ( kind = 4 ) c integer ( kind = 4 ) e integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,2*nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j integer ( kind = 4 ) n integer ( kind = 4 ) ne integer ( kind = 4 ) nw integer ( kind = 4 ) s integer ( kind = 4 ) se integer ( kind = 4 ) sw integer ( kind = 4 ) w ! ! Node labeling: ! ! NW---N--NE ! | \ | ! W C E ! | \ | ! SW---S--SE ! element = 0 do j = 1, nelemy do i = 1, nelemx sw = 2 * ( j - 1 ) * ( 2 * nelemx + 1 ) + 2 * ( i - 1 ) + 1 w = sw + 2 * nelemx + 1 nw = sw + 2 * ( 2 * nelemx + 1 ) s = sw + 1 c = sw + 1 + 2 * nelemx + 1 n = sw + 1 + 2 * ( 2 * nelemx + 1 ) se = sw + 2 e = sw + 2 + 2 * nelemx + 1 ne = sw + 2 + 2 * ( 2 * nelemx + 1 ) element = element + 1 element_node(1,element) = sw element_node(2,element) = se element_node(3,element) = nw element_node(4,element) = s element_node(5,element) = c element_node(6,element) = w element = element + 1 element_node(1,element) = ne element_node(2,element) = nw element_node(3,element) = se element_node(4,element) = n element_node(5,element) = c element_node(6,element) = e end do end do return end subroutine grid_t6_element_num ( nelemx, nelemy, element_num ) !*****************************************************************************80 ! !! GRID_T6_ELEMENT_NUM counts the elements in a grid of 6 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! ELEMENT_NUM = 2 * NELEMX * NELEMY = 12 ! ! Modified: ! ! 14 April 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NUM, the number of elements in the grid. ! implicit none integer ( kind = 4 ) element_num integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy element_num = 2 * nelemx * nelemy return end subroutine grid_t6_node_num ( nelemx, nelemy, node_num ) !*****************************************************************************80 ! !! GRID_T6_NODE_NUM counts the nodes in a grid of 6 node triangles. ! ! Modified: ! ! 05 October 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) NODE_NUM, the number of nodes in the grid. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ) node_num node_num = ( 2 * nelemx + 1 ) * ( 2 * nelemy + 1 ) return end subroutine grid_t10_element ( nelemx, nelemy, element_node ) !*****************************************************************************80 ! !! GRID_T10_ELEMENT produces an element grid of pairs of 10 node triangles. ! ! Example: ! ! Input: ! ! NELEMX = 2, NELEMY = 2 ! ! Output: ! ! ELEMENT_NODE = ! 1, 2, 3, 4, 10, 16, 22, 15, 8, 9; ! 25, 24, 23, 22, 16, 10, 4, 11, 18, 17; ! 4, 5, 6, 7, 13, 19, 25, 18, 11, 12; ! 28, 27, 26, 25, 19, 13, 7, 14, 21, 20; ! 22, 23, 24, 25, 31, 37, 43, 36, 29, 30; ! 46, 45, 44, 43, 37, 31, 25, 32, 39, 38; ! 25, 26, 27, 28, 34, 40, 46, 39, 31, 33; ! 49, 48, 47, 46, 40, 34, 28, 35, 42, 41. ! ! Grid: ! ! 43-44-45-46-47-48-49 ! |\ 6 |\ 8 | ! | \ | \ | ! 36 37 38 39 40 41 42 ! | \ | \ | ! | \ | \ | ! 29 30 31 32 33 34 35 ! | \ | \ | ! | 5 \| 7 \| ! 22-23-24-25-26-27-28 ! |\ 2 |\ 4 | ! | \ | \ | ! 15 16 17 18 19 20 21 ! | \ | \ | ! | \ | \ | ! 8 9 10 11 12 13 14 ! | \ | \ | ! | 1 \| 3 \| ! 1--2--3--4--5--6--7 ! ! Reference Element T10: ! ! | ! 1 10 ! | |\ ! | | \ ! | 8 9 ! | | \ ! S | \ ! | 5 6 7 ! | | \ ! | | \ ! 0 1--2--3--4 ! | ! +--0----R---1--> ! ! Modified: ! ! 09 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer ( kind = 4 ) NELEMX, NELEMY, the number of elements along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer ( kind = 4 ) ELEMENT_NODE(10,2*NELEMX*NELEMY), the nodes that form ! each element. ! implicit none integer ( kind = 4 ) nelemx integer ( kind = 4 ) nelemy integer ( kind = 4 ), parameter :: element_order = 10 integer ( kind = 4 ) base integer ( kind = 4 ) element integer ( kind = 4 ) element_node(element_order,2*nelemx*nelemy) integer ( kind = 4 ) i integer ( kind = 4 ) j element = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2 element = element +