subroutine basis_q4 ( dphidx, dphidy, phi, xl, xm, xr, yb, ym, yt ) ! !******************************************************************************* ! !! BASIS_Q4 evalutes basis functions for a rectangular bilinear element. ! ! ! Discussion: ! ! The routine is given the corners of a rectangular bilinear element. ! It then evaluates the basis functions associated with each corner, ! and their derivatives with respect to X and Y. ! ! The "local node" numbering is as follows: ! ! 3-------4 <-- Y = YT ! | | ! | | ! 1-------2 <-- Y = YB ! ! ^ ^ ! | | ! x = XL x = XR ! ! Modified: ! ! 28 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real DPHIDX(4), the value of the X-derivatives of the basis ! functions associated with local nodes 1, 2, 3 and 4 at (XM,YM). ! ! Output, real DPHIDY(4), the value of the Y-derivatives of the basis ! functions associated with local nodes 1, 2, 3 and 4 at (XM,YM). ! ! Output, real PHI(4), the value of the basis functions associated with ! local nodes 1, 2, 3 and 4 at (XM,YM). ! ! Input, real XL, the X coordinate of the left boundary of the element. ! ! Input, real XM. ! XM is the X coordinate of the point at which the basis functions ! and their derivatives should be evaluated. It need not be the ! center point, but generally should lie within the element. ! ! Input, real XR, the X coordinate of the right boundary of the element. ! ! Input, real YB, the Y coordinate of the bottom boundary of the element. ! ! Input, real YM. ! YM is the Y coordinate of the point at which the basis functions ! and their derivatives should be evaluated. It need not be the ! center point, but generally should lie within the element. ! ! Input, real YT, the Y coordinate of the top boundary of the element. ! implicit none ! real dphidx(4) real dphidy(4) real phi(4) real xl real xm real xr real yb real ym real yt ! phi(1) = ( xm - xr ) * ( ym - yt ) / ( ( xl - xr ) * ( yb - yt ) ) phi(2) = ( xm - xl ) * ( ym - yt ) / ( ( xr - xl ) * ( yb - yt ) ) phi(3) = ( xm - xr ) * ( ym - yb ) / ( ( xl - xr ) * ( yt - yb ) ) phi(4) = ( xm - xl ) * ( ym - yb ) / ( ( xr - xl ) * ( yt - yb ) ) dphidx(1) = ( ym - yt ) / ( ( xl - xr ) * ( yb - yt ) ) dphidx(2) = ( ym - yt ) / ( ( xr - xl ) * ( yb - yt ) ) dphidx(3) = ( ym - yb ) / ( ( xl - xr ) * ( yt - yb ) ) dphidx(4) = ( ym - yb ) / ( ( xr - xl ) * ( yt - yb ) ) dphidy(1) = ( xm - xr ) / ( ( xl - xr ) * ( yb - yt ) ) dphidy(2) = ( xm - xl ) / ( ( xr - xl ) * ( yb - yt ) ) dphidy(3) = ( xm - xr ) / ( ( xl - xr ) * ( yt - yb ) ) dphidy(4) = ( xm - xl ) / ( ( xr - xl ) * ( yt - yb ) ) return end subroutine div_q4 ( div, m, n, u, v, vort, xhi, xlo, yhi, ylo ) ! !******************************************************************************* ! !! 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(I), Y(J) ) for I = 1,M and J = 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. ! ! Modified: ! ! 28 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, real 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. ! ! 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 U(M,N), the value of the ! first component of a vector quantity whose divergence and ! vorticity are desired. A common example would be that U ! represents the horizontal velocity component of a flow field. ! ! Input, real V(M,N), the value of the ! second component of a vector quantity whose divergence and ! vorticity are desired. A common example would be that V ! represents the vertical velocity component of a flow field. ! ! Output, real 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. ! ! Input, real XHI, the X coordinate of the rightmost data column. ! ! Input, real XLO, the X coordinate of the leftmost data column. ! XHI and XLO must be distinct. ! ! Input, real YHI, the Y coordinate of the uppermost data row. ! ! Input, real YLO, the Y coordinate of the lowermost data row. ! YHI and YLO must be distinct. ! implicit none ! integer m integer n ! real div(m-1,n-1) real dphidx(4) real dphidy(4) integer i integer j real phi(4) real u(m,n) real v(m,n) real vort(m-1,n-1) real xhi real xl real xlo real xm real xr real yhi real yb real ylo real ym real yt ! if ( m <= 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_Q4 - Fatal error!' write ( *, '(a)' ) ' M must be at least 2,' write ( *, '(a,i6)' ) ' 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,i6)' ) ' but the input value of N is ', n stop end if if ( xhi == xlo ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_Q4 - Fatal error!' write ( *, '(a)' ) ' XHI and XLO must be distinct,' write ( *, '(a,g14.6)' ) ' but the input value of XLO is ', xlo write ( *, '(a,g14.6)' ) ' and the input value of XHI is ', xhi stop end if if ( yhi == ylo ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'DIV_Q4 - Fatal error!' write ( *, '(a)' ) ' YHI and YLO must be distinct,' write ( *, '(a,g14.6)' ) ' but the input value of YLO is ', ylo write ( *, '(a,g14.6)' ) ' and the input value of YHI is ', yhi stop end if do i = 1, m-1 yb = ( real (2*m-2*i ) * ylo + real (2*i-2) * yhi ) / real ( 2*m-2 ) ym = ( real (2*m-2*i-1) * ylo + real (2*i-1) * yhi ) / real ( 2*m-2 ) yt = ( real (2*m-2*i-2) * ylo + real (2*i ) * yhi ) / real ( 2*m-2 ) do j = 1, n-1 xl = ( real (2*n-2*j ) * xlo + real (2*j-2) * xhi ) / real ( 2*n-2 ) xm = ( real (2*n-2*j-1) * xlo + real (2*j-1) * xhi ) / real ( 2*n-2 ) xr = ( real (2*n-2*j-2) * xlo + real (2*j ) * xhi ) / real ( 2*n-2 ) call basis_q4 ( dphidx, dphidy, phi, xl, xm, xr, yb, ym, yt ) ! ! 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 ) * phi(3) + u(i+1,j+1) * phi(4) ! ! vmid = v(i ,j ) * phi(1) + v(i ,j+1) * phi(2) & ! + v(i+1,j ) * phi(3) + v(i+1,j+1) * phi(4) ! div(i,j) = u(i ,j ) * dphidx(1) + u(i ,j+1) * dphidx(2) & + u(i+1,j ) * dphidx(3) + u(i+1,j+1) * dphidx(4) & + v(i ,j ) * dphidy(1) + v(i ,j+1) * dphidy(2) & + v(i+1,j ) * dphidy(3) + v(i+1,j+1) * dphidy(4) vort(i,j) = v(i ,j ) * dphidx(1) + v(i ,j+1) * dphidx(2) & + v(i+1,j ) * dphidx(3) + v(i+1,j+1) * dphidx(4) & - u(i ,j ) * dphidy(1) - u(i ,j+1) * dphidy(2) & - u(i+1,j ) * dphidy(3) - u(i+1,j+1) * dphidy(4) end do end do return end function element_code ( i ) ! !******************************************************************************* ! !! 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 T6 6 node quadratic triangle; ! 9 T10 10 node cubic triangle. ! ! Modified: ! ! 04 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer I, the number of the element. ! ! Output, character ( len = 4 ) ELEMENT_CODE, the code for the element. ! implicit none ! character ( len = 4 ) element_code integer 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 = 'T6' else if ( i == 9 ) then element_code = 'T10' else element_code = '????' end if return end subroutine grid ( code, maxelem, n, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID returns the grid associated with any available element. ! ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element desired. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3', ! 'T6' and 'T10'. ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer N, the order of the element. ! ! Input, integer NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer NODES(N,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem integer n ! character ( len = * ) code integer nelemx integer nelemy integer nodes(n,maxelem) ! if ( code == 'Q4' ) then call grid_q4 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'Q8' ) then call grid_q8 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'Q9' ) then call grid_q9 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'Q12' ) then call grid_q12 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'Q16' ) then call grid_q16 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'QL' ) then call grid_ql ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'T3' ) then call grid_t3 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'T6' ) then call grid_t6 ( maxelem, nelemx, nelemy, nodes ) else if ( code == 'T10' ) then call grid_t10 ( maxelem, nelemx, nelemy, nodes ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID - Fatal error!' write ( *, '(a,a)' ) ' Illegal value of CODE = ', trim ( code ) stop end if return end subroutine grid_shape_2d ( n, a, n1, n2 ) ! !******************************************************************************* ! !! 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 N, the number of data values. ! ! Input, real A(N), the data, which should have the properties ! described above. ! ! Output, integer N1, N2, the "shape" of the data in the array. ! implicit none ! integer n ! real a(n) integer i integer n1 integer 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_test ( code ) ! !******************************************************************************* ! !! GRID_TEST tests the grid routines. ! ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, the code for the element. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! implicit none ! integer, parameter :: maxelem = 12 ! character ( len = * ) code integer i integer ielem integer n integer nelem integer nelemx integer nelemy integer nodes(16*maxelem) integer order_code integer width ! ! NODES is defined as a vector rather than a two dimensional array, ! so that we can handle the various cases using a single array. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_TEST' write ( *, '(a,a)' ) ' Test the grid routine for element ', trim ( code ) nelemx = 3 nelemy = 2 nelem = nelemx * nelemy n = order_code ( code ) call grid ( code, maxelem, n, nelemx, nelemy, nodes ) do ielem = 1, nelem write ( *, '(i3,3x,20i3)' ) ielem, ( nodes((ielem-1)*n+i), i = 1, n ) end do call grid_width ( maxelem, n, nelem, nodes, width ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'Grid width is ', width return end subroutine grid_width ( maxelm, n, nelem, node, width ) ! !******************************************************************************* ! !! GRID_WIDTH computes the width of a given grid. ! ! ! Definition: ! ! The grid width is defined to be 1 plus the maximum absolute ! difference of global indices of nodes in the same element. ! ! Example: ! ! For the following simple grid, the grid width is 14. ! ! 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 ! ! Modified: ! ! 12 October 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELM, the maximum number of elements. ! ! Input, integer N, the order of the elements. ! ! Input, integer NELEM, the number of elements. ! ! Input, integer NODE(N,MAXELM), the nodes that make up each element. ! ! Output, integer WIDTH, the grid width. ! implicit none ! integer maxelm integer n ! integer ielem integer inode1 integer inode2 integer ip1 integer ip2 integer nelem integer node(n,maxelm) integer width ! width = 0 do ielem = 1, nelem do inode1 = 1, n ip1 = node(inode1,ielem) do inode2 = 1, n ip2 = node(inode2,ielem) width = max ( width, abs ( ip1 - ip2 ) ) end do end do end do return end subroutine grid_q4 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_Q4 produces a grid of 4 node quadrilaterals. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 1, 2, 5, 6; ! 2, 3, 6, 7; ! 3, 4, 7, 8; ! 5, 6, 9, 10; ! 6, 7, 10, 11; ! 7, 8, 11, 12. ! ! Diagram: ! ! 9---10---11---12 ! | | | | ! | | | | ! | 4 | 5 | 6 | ! | | | | ! 5----6----7----8 ! | | | | ! | | | | ! | 1 | 2 | 3 | ! | | | | ! 1----2----3----4 ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer NODES(4,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(4,maxelem) ! nelem = nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_Q4 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx ielem = ielem + 1 nodes(1,ielem) = ( j - 1 ) * ( nelemx + 1 ) + i nodes(2,ielem) = ( j - 1 ) * ( nelemx + 1 ) + i + 1 nodes(3,ielem) = j * ( nelemx + 1 ) + i nodes(4,ielem) = j * ( nelemx + 1 ) + i + 1 end do end do return end subroutine grid_q8 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_Q8 produces a grid of 8 node quadrilaterals. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 3, 14, 12, 1, 2, 9, 13, 8; ! 5, 16, 14, 3, 4, 10, 15, 9; ! 7, 18, 16, 5, 6, 11, 17, 10; ! 14, 25, 23, 12, 13, 20, 24, 19; ! 16, 27, 25, 14, 15, 21, 26, 20; ! 18, 29, 27, 16, 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 ! ! Modified: ! ! 12 October 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer NODES(8,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(8,maxelem) ! nelem = nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_Q8 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx ielem = ielem + 1 base = ( j - 1 ) * ( 3 * nelemx + 2 ) + 2 * i - 1 nodes(1,ielem) = base + 2 nodes(2,ielem) = base + ( 3 * nelemx + 2 ) + 2 nodes(3,ielem) = base + ( 3 * nelemx + 2 ) nodes(4,ielem) = base nodes(5,ielem) = base + 1 nodes(6,ielem) = base + 2 * nelemx + 2 - i + 1 nodes(7,ielem) = base + ( 3 * nelemx + 2 ) + 1 nodes(8,ielem) = base + 2 * nelemx + 2 - i end do end do return end subroutine grid_q9 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_Q9 produces a grid of 9 node quadrilaterals. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 3, 17, 15, 1, 2, 10, 16, 8, 9; ! 5, 19, 17, 3, 4, 12, 18, 10, 11; ! 7, 21, 19, 5, 6, 14, 20, 12, 13; ! 17, 31, 29, 15, 16, 24, 30, 22, 23; ! 19, 33, 31, 17, 18, 26, 32, 24, 25; ! 21, 35, 33, 19, 20, 28, 34, 26, 27. ! ! Diagram: ! ! 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 ! ! Modified: ! ! 11 October 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer NODES(9,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(9,maxelem) ! nelem = nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_Q9 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx ielem = ielem + 1 base = ( 2 * j - 2 ) * ( 2 * nelemx + 1 ) + 2 * i - 1 nodes(1,ielem) = base + 2 nodes(2,ielem) = base + 2 * ( 2 * nelemx + 1 ) + 2 nodes(3,ielem) = base + 2 * ( 2 * nelemx + 1 ) nodes(4,ielem) = base nodes(5,ielem) = base + 1 nodes(6,ielem) = base + ( 2 * nelemx + 1 ) + 2 nodes(7,ielem) = base + 2 * ( 2 * nelemx + 1 ) + 1 nodes(8,ielem) = base + ( 2 * nelemx + 1 ) nodes(9,ielem) = base + ( 2 * nelemx + 1 ) + 1 end do end do return end subroutine grid_q12 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_Q12 produces a grid of 12 node quadrilaterals. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 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. ! ! Diagram: ! ! 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 ! ! Modified: ! ! 07 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of quadrilaterals along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer NODES(12,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(12,maxelem) ! nelem = nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_Q12 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx ielem = ielem + 1 base = ( j - 1 ) * ( 5 * nelemx + 3 ) + 1 nodes(1,ielem) = base + ( i - 1 ) * 3 nodes(2,ielem) = base + ( i - 1 ) * 3 + 1 nodes(3,ielem) = base + ( i - 1 ) * 3 + 2 nodes(4,ielem) = base + ( i - 1 ) * 3 + 3 nodes(5,ielem) = base + 3 * nelemx + i nodes(6,ielem) = base + 3 * nelemx + i + 1 nodes(7,ielem) = base + 4 * nelemx + i + 1 nodes(8,ielem) = base + 4 * nelemx + i + 2 nodes(9,ielem) = base + 5 * nelemx + 3 * i nodes(10,ielem) = base + 5 * nelemx + 3 * i + 1 nodes(11,ielem) = base + 5 * nelemx + 3 * i + 2 nodes(12,ielem) = base + 5 * nelemx + 3 * i + 3 end do end do return end subroutine grid_q16 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_Q16 produces a grid of 16 node quadrilaterals. ! ! ! Example: ! ! Input: ! ! NELEMX = 2, NELEMY = 2 ! ! Output: ! ! NODES = ! 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. ! ! ! Diagram: ! ! 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 ! ! Modified: ! ! 08 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of triangles along the ! X and Y directions. The number of elements generated will be ! NELEMX * NELEMY. ! ! Output, integer NODES(16,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(16,maxelem) ! nelem = nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_Q16 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2 ielem = ielem + 1 nodes( 1,ielem) = base nodes( 2,ielem) = base + 1 nodes( 3,ielem) = base + 2 nodes( 4,ielem) = base + 3 nodes( 5,ielem) = base + ( 3 * nelemx + 1 ) nodes( 6,ielem) = base + ( 3 * nelemx + 1 ) + 1 nodes( 7,ielem) = base + ( 3 * nelemx + 1 ) + 2 nodes( 8,ielem) = base + ( 3 * nelemx + 1 ) + 3 nodes( 9,ielem) = base + 2 * ( 3 * nelemx + 1 ) nodes(10,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 1 nodes(11,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 2 nodes(12,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 3 nodes(13,ielem) = base + 3 * ( 3 * nelemx + 1 ) nodes(14,ielem) = base + 3 * ( 3 * nelemx + 1 ) + 1 nodes(15,ielem) = base + 3 * ( 3 * nelemx + 1 ) + 2 nodes(16,ielem) = base + 3 * ( 3 * nelemx + 1 ) + 3 end do end do return end subroutine grid_ql ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_QL produces a grid of 6 node quadratics/linears. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 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. ! ! Diagram: ! ! 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 MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of quadrilaterals 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 NODES(6,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(6,maxelem) ! nelem = nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_QL - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx ielem = ielem + 1 base = ( j - 1 ) * ( 2 * nelemx + 1 ) + 2 * i - 1 nodes(1,ielem) = base nodes(2,ielem) = base + 1 nodes(3,ielem) = base + 2 nodes(4,ielem) = base + ( 2 * nelemx + 1 ) nodes(5,ielem) = base + ( 2 * nelemx + 1 ) + 1 nodes(6,ielem) = base + ( 2 * nelemx + 1 ) + 2 end do end do return end subroutine grid_t3 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_T3 produces a grid of pairs of 3 node triangles. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 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. ! ! Diagram: ! ! 9---10---11---12 ! |\ 8 |\10 |\12 | ! | \ | \ | \ | ! | \ | \ | \ | ! | 7\| 9\| 11\| ! 5----6----7----8 ! |\ 2 |\ 4 |\ 6 | ! | \ | \ | \ | ! | \ | \ | \ | ! | 1\| 3\| 5\| ! 1----2----3----4 ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least 2 * NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of triangles along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer NODES(3,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(3,maxelem) ! nelem = 2 * nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_T3 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx ielem = ielem + 1 nodes(1,ielem) = ( j - 1 ) * ( nelemx + 1 ) + i nodes(2,ielem) = ( j - 1 ) * ( nelemx + 1 ) + i + 1 nodes(3,ielem) = j * ( nelemx + 1 ) + i ielem = ielem + 1 nodes(1,ielem) = j * ( nelemx + 1 ) + i + 1 nodes(2,ielem) = j * ( nelemx + 1 ) + i nodes(3,ielem) = ( j - 1 ) * ( nelemx + 1 ) + i + 1 end do end do return end subroutine grid_t6 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_T6 produces a grid of pairs of 6 node triangles. ! ! ! Example: ! ! Input: ! ! NELEMX = 3, NELEMY = 2 ! ! Output: ! ! NODES = ! 3, 15, 1, 2, 9, 8; ! 15, 3, 17, 16, 9, 10; ! 5, 17, 3 4, 11, 10; ! 17, 5, 19, 18, 11, 12; ! 7, 19, 5, 6, 13, 12; ! 19, 7, 21, 20, 13, 14; ! 17, 29, 15, 16, 23, 22; ! 29, 17, 31, 30, 23, 24; ! 19, 31, 17, 18, 25, 24; ! 31, 19, 33, 32, 25, 26; ! 21, 33, 19, 20, 27, 26; ! 33, 21, 35, 34, 27, 28. ! ! Diagram: ! ! 29-30-31-32-33-34-35 ! |\ 8 |\10 |\12 | ! | \ | \ | \ | ! 22 23 24 25 26 27 28 ! | \ | \ | \ | ! | 7 \| 9 \| 11 \| ! 15-16-17-18-19-20-21 ! |\ 2 |\ 4 |\ 6 | ! | \ | \ | \ | ! 8 9 10 11 12 13 14 ! | \ | \ | \ | ! | 1 \| 3 \| 5 \| ! 1--2--3--4--5--6--7 ! ! Modified: ! ! 12 October 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least 2 * NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of triangles along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer NODES(6,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(6,maxelem) ! nelem = 2 * nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_T6 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * 2 * ( 2 * nelemx + 1 ) + 2 * i - 1 ielem = ielem + 1 nodes(1,ielem) = base + 2 nodes(2,ielem) = base + 2 * ( 2 * nelemx + 1 ) nodes(3,ielem) = base nodes(4,ielem) = base + 1 nodes(5,ielem) = base + ( 2 * nelemx + 1 ) + 1 nodes(6,ielem) = base + ( 2 * nelemx + 1 ) ielem = ielem + 1 nodes(1,ielem) = base + 2 * ( 2 * nelemx + 1 ) nodes(2,ielem) = base + 2 nodes(3,ielem) = base + 2 * ( 2 * nelemx + 1 ) + 2 nodes(4,ielem) = base + 2 * ( 2 * nelemx + 1 ) + 1 nodes(5,ielem) = base + ( 2 * nelemx + 1 ) + 1 nodes(6,ielem) = base + ( 2 * nelemx + 1 ) + 2 end do end do return end subroutine grid_t10 ( maxelem, nelemx, nelemy, nodes ) ! !******************************************************************************* ! !! GRID_T10 produces a grid of pairs of 10 node triangles. ! ! ! Example: ! ! Input: ! ! NELEMX = 2, NELEMY = 2 ! ! Output: ! ! NODES = ! 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. ! ! ! Diagram: ! ! 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 ! ! Modified: ! ! 09 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer MAXELEM, the maximum number of elements for which ! storage is provided, which must be at least 2 * NELEMX * NELEMY. ! ! Input, integer NELEMX, NELEMY, the number of triangles along the ! X and Y directions. The number of elements generated will be ! 2 * NELEMX * NELEMY. ! ! Output, integer NODES(10,MAXELEM); NODES(I,J) contains the index ! of the I-th node of the J-th element. ! implicit none ! integer maxelem ! integer base integer ielem integer i integer j integer nelem integer nelemx integer nelemy integer nodes(10,maxelem) ! nelem = 2 * nelemx * nelemy if ( nelem > maxelem ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'GRID_T10 - Fatal error!' write ( *, '(a)' ) ' Not enough storage for NODE array.' write ( *, '(a,i6)' ) ' Increase MAXELEM to ', nelem stop end if ielem = 0 do j = 1, nelemy do i = 1, nelemx base = ( j - 1 ) * 3 * ( 3 * nelemx + 1 ) + 3 * i - 2 ielem = ielem + 1 nodes( 1,ielem) = base nodes( 2,ielem) = base + 1 nodes( 3,ielem) = base + 2 nodes( 4,ielem) = base + 3 nodes( 5,ielem) = base + ( 3 * nelemx + 1 ) + 2 nodes( 6,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 1 nodes( 7,ielem) = base + 3 * ( 3 * nelemx + 1 ) nodes( 8,ielem) = base + 2 * ( 3 * nelemx + 1 ) nodes( 9,ielem) = base + ( 2 * nelemx + 1 ) + 2 nodes(10,ielem) = base + ( 2 * nelemx + 1 ) + 3 ielem = ielem + 1 nodes( 1,ielem) = base + 3 * ( 3 * nelemx + 1 ) + 3 nodes( 2,ielem) = base + 3 * ( 3 * nelemx + 1 ) + 2 nodes( 3,ielem) = base + 3 * ( 3 * nelemx + 1 ) + 1 nodes( 4,ielem) = base + 3 * ( 3 * nelemx + 1 ) nodes( 5,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 1 nodes( 6,ielem) = base + ( 3 * nelemx + 1 ) + 2 nodes( 7,ielem) = base + 3 nodes( 8,ielem) = base + ( 3 * nelemx + 1 ) + 3 nodes( 9,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 3 nodes(10,ielem) = base + 2 * ( 3 * nelemx + 1 ) + 2 end do end do return end subroutine interp ( code, dtdr, dtds, maxn, n, r, s, t, ubase, u, dudr, duds ) ! !******************************************************************************* ! !! INTERP interpolates a quantity in an element from basis node values. ! ! ! Modified: ! ! 27 October 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! ! Output, real DTDR(MAXN), DTDS(MAXN), the derivatives of the ! basis functions at (R,S). ! ! Input, integer MAXN, the maximum value of N. ! ! Output, integer N, the order of the element. ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(MAXN), the value of the basis functions at (R,S). ! ! Input, real UBASE(MAXN), the value of the quantity at the basis nodes. ! ! Output, real U, DUDR, DUDS, the interpolated value of the ! quantity and its derivatives at the point (R,S). ! implicit none ! integer maxn ! character ( len = * ) code real dtdr(maxn) real dtds(maxn) real dudr real duds integer n real r real s real t(maxn) real u real ubase(maxn) ! call shape ( code, n, r, s, t, dtdr, dtds ) u = dot_product ( ubase(1:n), t(1:n) ) dudr = dot_product ( ubase(1:n), dtdr(1:n) ) duds = dot_product ( ubase(1:n), dtds(1:n) ) return end subroutine legendre_com ( norder, xtab, weight ) ! !******************************************************************************* ! !! LEGENDRE_COM computes abscissas and weights for Gauss-Legendre quadrature. ! ! ! Integration interval: ! ! [ -1, 1 ] ! ! Weight function: ! ! 1. ! ! Integral to approximate: ! ! Integral ( -1 <= X <= 1 ) F(X) dX. ! ! Approximate integral: ! ! Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( XTAB(I) ). ! ! Modified: ! ! 04 November 1998 ! ! Parameters: ! ! Input, integer NORDER, the order of the rule. ! NORDER must be greater than 0. ! ! Output, real XTAB(NORDER), the abscissas of the rule. ! ! Output, real WEIGHT(NORDER), the weights of the rule. ! The weights are positive, symmetric, and should sum to 2. ! implicit none ! integer norder ! real d1 real d2pn real d3pn real d4pn real dp real dpn real e1 real fx real h integer i integer iback integer k integer m integer mp1mi integer ncopy integer nmove real p real, parameter :: pi = & 3.14159265358979323846264338327950288419716939937510E+00 real pk real pkm1 real pkp1 real t real u real v real x0 real xtab(norder) real xtemp real weight(norder) ! if ( norder < 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEGENDRE_COM - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of NORDER = ', norder stop end if e1 = real ( norder * ( norder + 1 ) ) m = ( norder + 1 ) / 2 do i = 1, ( norder + 1 ) / 2 mp1mi = m + 1 - i t = real ( 4 * i - 1 ) * pi / real ( 4 * norder + 2 ) x0 = cos(t) * ( 1.0E+00 - ( 1.0E+00 - 1.0E+00 / & real ( norder ) ) / real ( 8 * norder**2) ) pkm1 = 1.0E+00 pk = x0 do k = 2, norder pkp1 = 2.0E+00 * x0 * pk - pkm1 - ( x0 * pk - pkm1 ) / real ( k ) pkm1 = pk pk = pkp1 end do d1 = real ( norder ) * ( pkm1 - x0 * pk ) dpn = d1 / ( 1.0E+00 - x0**2 ) d2pn = ( 2.0E+00 * x0 * dpn - e1 * pk ) / ( 1.0E+00 - x0**2 ) d3pn = ( 4.0E+00 * x0 * d2pn + ( 2.0E+00 - e1 ) * dpn ) / & ( 1.0E+00 - x0**2 ) d4pn = ( 6.0E+00 * x0 * d3pn + ( 6.0E+00 - e1 ) * d2pn ) & / ( 1.0E+00 - x0**2 ) u = pk / dpn v = d2pn / dpn ! ! Initial approximation H: ! h = - u * ( 1.0E+00 + 0.5E+00 * u * ( v + u * ( v**2 - d3pn & / ( 3.0E+00 * dpn ) ) ) ) ! ! Refine H using one step of Newton's method: ! p = pk + h * ( dpn + 0.5E+00 * h * ( d2pn + h / 3.0E+00 & * ( d3pn + 0.25E+00 * h * d4pn ) ) ) dp = dpn + h * ( d2pn + 0.5E+00 * h * ( d3pn + h * d4pn / 3.0E+00 ) ) h = h - p / dp xtemp = x0 + h xtab(mp1mi) = xtemp fx = d1 - h * e1 * ( pk + 0.5E+00 * h * ( dpn + h / 3.0E+00 & * ( d2pn + 0.25E+00 * h * ( d3pn + 0.2E+00 * h * d4pn ) ) ) ) weight(mp1mi) = 2.0E+00 * ( 1.0E+00 - xtemp**2 ) / fx**2 end do if ( mod ( norder, 2 ) == 1 ) then xtab(1) = 0.0E+00 end if ! ! Shift the data up. ! nmove = ( norder + 1 ) / 2 ncopy = norder - nmove do i = 1, nmove iback = norder + 1 - i xtab(iback) = xtab(iback-ncopy) weight(iback) = weight(iback-ncopy) end do ! ! Reflect values for the negative abscissas. ! do i = 1, norder - nmove xtab(i) = - xtab(norder+1-i) weight(i) = weight(norder+1-i) end do return end subroutine legendre_set ( norder, xtab, weight ) ! !******************************************************************************* ! !! LEGENDRE_SET sets abscissas and weights for Gauss-Legendre quadrature. ! ! ! Integration interval: ! ! [ -1, 1 ] ! ! Weight function: ! ! 1.0D+00 ! ! Integral to approximate: ! ! Integral ( -1 <= X <= 1 ) F(X) dX ! ! Approximate integral: ! ! Sum ( 1 <= I <= NORDER ) WEIGHT(I) * F ( XTAB(I) ) ! ! Precision: ! ! The quadrature rule will integrate exactly all polynomials up to ! X**(2*NORDER-1). ! ! Note: ! ! The abscissas of the rule are the zeroes of the Legendre polynomial ! P(NORDER)(X). ! ! The integral produced by a Gauss-Legendre rule is equal to the ! integral of the unique polynomial of degree NORDER-1 which ! agrees with the function at the NORDER abscissas of the rule. ! ! Reference: ! ! Abramowitz and Stegun, ! Handbook of Mathematical Functions, ! National Bureau of Standards, 1964. ! ! Vladimir Krylov, ! Approximate Calculation of Integrals, ! MacMillan, 1962. ! ! Arthur Stroud and Don Secrest, ! Gaussian Quadrature Formulas, ! Prentice Hall, 1966. ! ! Daniel Zwillinger, editor, ! Standard Mathematical Tables and Formulae, ! 30th Edition, ! CRC Press, 1996. ! ! Modified: ! ! 18 December 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer NORDER, the order of the rule. ! NORDER must be between 1 and 20, 32 or 64. ! ! Output, double precision XTAB(NORDER), the abscissas of the rule. ! ! Output, double precision WEIGHT(NORDER), the weights of the rule. ! The weights are positive, symmetric and should sum to 2. ! implicit none ! integer norder ! double precision xtab(norder) double precision weight(norder) ! if ( norder == 1 ) then xtab(1) = 0.0D+00 weight(1) = 2.0D+00 else if ( norder == 2 ) then xtab(1) = - 0.577350269189625764509148780502D+00 xtab(2) = 0.577350269189625764509148780502D+00 weight(1) = 1.0D+00 weight(2) = 1.0D+00 else if ( norder == 3 ) then xtab(1) = - 0.774596669241483377035853079956D+00 xtab(2) = 0.0D+00 xtab(3) = 0.774596669241483377035853079956D+00 weight(1) = 5.0D+00 / 9.0D+00 weight(2) = 8.0D+00 / 9.0D+00 weight(3) = 5.0D+00 / 9.0D+00 else if ( norder == 4 ) then xtab(1) = - 0.861136311594052575223946488893D+00 xtab(2) = - 0.339981043584856264802665759103D+00 xtab(3) = 0.339981043584856264802665759103D+00 xtab(4) = 0.861136311594052575223946488893D+00 weight(1) = 0.347854845137453857373063949222D+00 weight(2) = 0.652145154862546142626936050778D+00 weight(3) = 0.652145154862546142626936050778D+00 weight(4) = 0.347854845137453857373063949222D+00 else if ( norder == 5 ) then xtab(1) = - 0.906179845938663992797626878299D+00 xtab(2) = - 0.538469310105683091036314420700D+00 xtab(3) = 0.0D+00 xtab(4) = 0.538469310105683091036314420700D+00 xtab(5) = 0.906179845938663992797626878299D+00 weight(1) = 0.236926885056189087514264040720D+00 weight(2) = 0.478628670499366468041291514836D+00 weight(3) = 0.568888888888888888888888888889D+00 weight(4) = 0.478628670499366468041291514836D+00 weight(5) = 0.236926885056189087514264040720D+00 else if ( norder == 6 ) then xtab(1) = - 0.932469514203152027812301554494D+00 xtab(2) = - 0.661209386466264513661399595020D+00 xtab(3) = - 0.238619186083196908630501721681D+00 xtab(4) = 0.238619186083196908630501721681D+00 xtab(5) = 0.661209386466264513661399595020D+00 xtab(6) = 0.932469514203152027812301554494D+00 weight(1) = 0.171324492379170345040296142173D+00 weight(2) = 0.360761573048138607569833513838D+00 weight(3) = 0.467913934572691047389870343990D+00 weight(4) = 0.467913934572691047389870343990D+00 weight(5) = 0.360761573048138607569833513838D+00 weight(6) = 0.171324492379170345040296142173D+00 else if ( norder == 7 ) then xtab(1) = - 0.949107912342758524526189684048D+00 xtab(2) = - 0.741531185599394439863864773281D+00 xtab(3) = - 0.405845151377397166906606412077D+00 xtab(4) = 0.0D+00 xtab(5) = 0.405845151377397166906606412077D+00 xtab(6) = 0.741531185599394439863864773281D+00 xtab(7) = 0.949107912342758524526189684048D+00 weight(1) = 0.129484966168869693270611432679D+00 weight(2) = 0.279705391489276667901467771424D+00 weight(3) = 0.381830050505118944950369775489D+00 weight(4) = 0.417959183673469387755102040816D+00 weight(5) = 0.381830050505118944950369775489D+00 weight(6) = 0.279705391489276667901467771424D+00 weight(7) = 0.129484966168869693270611432679D+00 else if ( norder == 8 ) then xtab(1) = - 0.960289856497536231683560868569D+00 xtab(2) = - 0.796666477413626739591553936476D+00 xtab(3) = - 0.525532409916328985817739049189D+00 xtab(4) = - 0.183434642495649804939476142360D+00 xtab(5) = 0.183434642495649804939476142360D+00 xtab(6) = 0.525532409916328985817739049189D+00 xtab(7) = 0.796666477413626739591553936476D+00 xtab(8) = 0.960289856497536231683560868569D+00 weight(1) = 0.101228536290376259152531354310D+00 weight(2) = 0.222381034453374470544355994426D+00 weight(3) = 0.313706645877887287337962201987D+00 weight(4) = 0.362683783378361982965150449277D+00 weight(5) = 0.362683783378361982965150449277D+00 weight(6) = 0.313706645877887287337962201987D+00 weight(7) = 0.222381034453374470544355994426D+00 weight(8) = 0.101228536290376259152531354310D+00 else if ( norder == 9 ) then xtab(1) = - 0.968160239507626089835576202904D+00 xtab(2) = - 0.836031107326635794299429788070D+00 xtab(3) = - 0.613371432700590397308702039341D+00 xtab(4) = - 0.324253423403808929038538014643D+00 xtab(5) = 0.0D+00 xtab(6) = 0.324253423403808929038538014643D+00 xtab(7) = 0.613371432700590397308702039341D+00 xtab(8) = 0.836031107326635794299429788070D+00 xtab(9) = 0.968160239507626089835576202904D+00 weight(1) = 0.812743883615744119718921581105D-01 weight(2) = 0.180648160694857404058472031243D+00 weight(3) = 0.260610696402935462318742869419D+00 weight(4) = 0.312347077040002840068630406584D+00 weight(5) = 0.330239355001259763164525069287D+00 weight(6) = 0.312347077040002840068630406584D+00 weight(7) = 0.260610696402935462318742869419D+00 weight(8) = 0.180648160694857404058472031243D+00 weight(9) = 0.812743883615744119718921581105D-01 else if ( norder == 10 ) then xtab(1) = - 0.973906528517171720077964012084D+00 xtab(2) = - 0.865063366688984510732096688423D+00 xtab(3) = - 0.679409568299024406234327365115D+00 xtab(4) = - 0.433395394129247190799265943166D+00 xtab(5) = - 0.148874338981631210884826001130D+00 xtab(6) = 0.148874338981631210884826001130D+00 xtab(7) = 0.433395394129247190799265943166D+00 xtab(8) = 0.679409568299024406234327365115D+00 xtab(9) = 0.865063366688984510732096688423D+00 xtab(10) = 0.973906528517171720077964012084D+00 weight(1) = 0.666713443086881375935688098933D-01 weight(2) = 0.149451349150580593145776339658D+00 weight(3) = 0.219086362515982043995534934228D+00 weight(4) = 0.269266719309996355091226921569D+00 weight(5) = 0.295524224714752870173892994651D+00 weight(6) = 0.295524224714752870173892994651D+00 weight(7) = 0.269266719309996355091226921569D+00 weight(8) = 0.219086362515982043995534934228D+00 weight(9) = 0.149451349150580593145776339658D+00 weight(10) = 0.666713443086881375935688098933D-01 else if ( norder == 11 ) then xtab(1) = - 0.978228658146056992803938001123D+00 xtab(2) = - 0.887062599768095299075157769304D+00 xtab(3) = - 0.730152005574049324093416252031D+00 xtab(4) = - 0.519096129206811815925725669459D+00 xtab(5) = - 0.269543155952344972331531985401D+00 xtab(6) = 0.0D+00 xtab(7) = 0.269543155952344972331531985401D+00 xtab(8) = 0.519096129206811815925725669459D+00 xtab(9) = 0.730152005574049324093416252031D+00 xtab(10) = 0.887062599768095299075157769304D+00 xtab(11) = 0.978228658146056992803938001123D+00 weight(1) = 0.556685671161736664827537204425D-01 weight(2) = 0.125580369464904624634694299224D+00 weight(3) = 0.186290210927734251426097641432D+00 weight(4) = 0.233193764591990479918523704843D+00 weight(5) = 0.262804544510246662180688869891D+00 weight(6) = 0.272925086777900630714483528336D+00 weight(7) = 0.262804544510246662180688869891D+00 weight(8) = 0.233193764591990479918523704843D+00 weight(9) = 0.186290210927734251426097641432D+00 weight(10) = 0.125580369464904624634694299224D+00 weight(11) = 0.556685671161736664827537204425D-01 else if ( norder == 12 ) then xtab(1) = - 0.981560634246719250690549090149D+00 xtab(2) = - 0.904117256370474856678465866119D+00 xtab(3) = - 0.769902674194304687036893833213D+00 xtab(4) = - 0.587317954286617447296702418941D+00 xtab(5) = - 0.367831498998180193752691536644D+00 xtab(6) = - 0.125233408511468915472441369464D+00 xtab(7) = 0.125233408511468915472441369464D+00 xtab(8) = 0.367831498998180193752691536644D+00 xtab(9) = 0.587317954286617447296702418941D+00 xtab(10) = 0.769902674194304687036893833213D+00 xtab(11) = 0.904117256370474856678465866119D+00 xtab(12) = 0.981560634246719250690549090149D+00 weight(1) = 0.471753363865118271946159614850D-01 weight(2) = 0.106939325995318430960254718194D+00 weight(3) = 0.160078328543346226334652529543D+00 weight(4) = 0.203167426723065921749064455810D+00 weight(5) = 0.233492536538354808760849898925D+00 weight(6) = 0.249147045813402785000562436043D+00 weight(7) = 0.249147045813402785000562436043D+00 weight(8) = 0.233492536538354808760849898925D+00 weight(9) = 0.203167426723065921749064455810D+00 weight(10) = 0.160078328543346226334652529543D+00 weight(11) = 0.106939325995318430960254718194D+00 weight(12) = 0.471753363865118271946159614850D-01 else if ( norder == 13 ) then xtab(1) = - 0.984183054718588149472829448807D+00 xtab(2) = - 0.917598399222977965206547836501D+00 xtab(3) = - 0.801578090733309912794206489583D+00 xtab(4) = - 0.642349339440340220643984606996D+00 xtab(5) = - 0.448492751036446852877912852128D+00 xtab(6) = - 0.230458315955134794065528121098D+00 xtab(7) = 0.0D+00 xtab(8) = 0.230458315955134794065528121098D+00 xtab(9) = 0.448492751036446852877912852128D+00 xtab(10) = 0.642349339440340220643984606996D+00 xtab(11) = 0.801578090733309912794206489583D+00 xtab(12) = 0.917598399222977965206547836501D+00 xtab(13) = 0.984183054718588149472829448807D+00 weight(1) = 0.404840047653158795200215922010D-01 weight(2) = 0.921214998377284479144217759538D-01 weight(3) = 0.138873510219787238463601776869D+00 weight(4) = 0.178145980761945738280046691996D+00 weight(5) = 0.207816047536888502312523219306D+00 weight(6) = 0.226283180262897238412090186040D+00 weight(7) = 0.232551553230873910194589515269D+00 weight(8) = 0.226283180262897238412090186040D+00 weight(9) = 0.207816047536888502312523219306D+00 weight(10) = 0.178145980761945738280046691996D+00 weight(11) = 0.138873510219787238463601776869D+00 weight(12) = 0.921214998377284479144217759538D-01 weight(13) = 0.404840047653158795200215922010D-01 else if ( norder == 14 ) then xtab(1) = - 0.986283808696812338841597266704D+00 xtab(2) = - 0.928434883663573517336391139378D+00 xtab(3) = - 0.827201315069764993189794742650D+00 xtab(4) = - 0.687292904811685470148019803019D+00 xtab(5) = - 0.515248636358154091965290718551D+00 xtab(6) = - 0.319112368927889760435671824168D+00 xtab(7) = - 0.108054948707343662066244650220D+00 xtab(8) = 0.108054948707343662066244650220D+00 xtab(9) = 0.319112368927889760435671824168D+00 xtab(10) = 0.515248636358154091965290718551D+00 xtab(11) = 0.687292904811685470148019803019D+00 xtab(12) = 0.827201315069764993189794742650D+00 xtab(13) = 0.928434883663573517336391139378D+00 xtab(14) = 0.986283808696812338841597266704D+00 weight(1) = 0.351194603317518630318328761382D-01 weight(2) = 0.801580871597602098056332770629D-01 weight(3) = 0.121518570687903184689414809072D+00 weight(4) = 0.157203167158193534569601938624D+00 weight(5) = 0.185538397477937813741716590125D+00 weight(6) = 0.205198463721295603965924065661D+00 weight(7) = 0.215263853463157790195876443316D+00 weight(8) = 0.215263853463157790195876443316D+00 weight(9) = 0.205198463721295603965924065661D+00 weight(10) = 0.185538397477937813741716590125D+00 weight(11) = 0.157203167158193534569601938624D+00 weight(12) = 0.121518570687903184689414809072D+00 weight(13) = 0.801580871597602098056332770629D-01 weight(14) = 0.351194603317518630318328761382D-01 else if ( norder == 15 ) then xtab(1) = - 0.987992518020485428489565718587D+00 xtab(2) = - 0.937273392400705904307758947710D+00 xtab(3) = - 0.848206583410427216200648320774D+00 xtab(4) = - 0.724417731360170047416186054614D+00 xtab(5) = - 0.570972172608538847537226737254D+00 xtab(6) = - 0.394151347077563369897207370981D+00 xtab(7) = - 0.201194093997434522300628303395D+00 xtab(8) = 0.0D+00 xtab(9) = 0.201194093997434522300628303395D+00 xtab(10) = 0.394151347077563369897207370981D+00 xtab(11) = 0.570972172608538847537226737254D+00 xtab(12) = 0.724417731360170047416186054614D+00 xtab(13) = 0.848206583410427216200648320774D+00 xtab(14) = 0.937273392400705904307758947710D+00 xtab(15) = 0.987992518020485428489565718587D+00 weight(1) = 0.307532419961172683546283935772D-01 weight(2) = 0.703660474881081247092674164507D-01 weight(3) = 0.107159220467171935011869546686D+00 weight(4) = 0.139570677926154314447804794511D+00 weight(5) = 0.166269205816993933553200860481D+00 weight(6) = 0.186161000015562211026800561866D+00 weight(7) = 0.198431485327111576456118326444D+00 weight(8) = 0.202578241925561272880620199968D+00 weight(9) = 0.198431485327111576456118326444D+00 weight(10) = 0.186161000015562211026800561866D+00 weight(11) = 0.166269205816993933553200860481D+00 weight(12) = 0.139570677926154314447804794511D+00 weight(13) = 0.107159220467171935011869546686D+00 weight(14) = 0.703660474881081247092674164507D-01 weight(15) = 0.307532419961172683546283935772D-01 else if ( norder == 16 ) then xtab(1) = - 0.989400934991649932596154173450D+00 xtab(2) = - 0.944575023073232576077988415535D+00 xtab(3) = - 0.865631202387831743880467897712D+00 xtab(4) = - 0.755404408355003033895101194847D+00 xtab(5) = - 0.617876244402643748446671764049D+00 xtab(6) = - 0.458016777657227386342419442984D+00 xtab(7) = - 0.281603550779258913230460501460D+00 xtab(8) = - 0.950125098376374401853193354250D-01 xtab(9) = 0.950125098376374401853193354250D-01 xtab(10) = 0.281603550779258913230460501460D+00 xtab(11) = 0.458016777657227386342419442984D+00 xtab(12) = 0.617876244402643748446671764049D+00 xtab(13) = 0.755404408355003033895101194847D+00 xtab(14) = 0.865631202387831743880467897712D+00 xtab(15) = 0.944575023073232576077988415535D+00 xtab(16) = 0.989400934991649932596154173450D+00 weight(1) = 0.271524594117540948517805724560D-01 weight(2) = 0.622535239386478928628438369944D-01 weight(3) = 0.951585116824927848099251076022D-01 weight(4) = 0.124628971255533872052476282192D+00 weight(5) = 0.149595988816576732081501730547D+00 weight(6) = 0.169156519395002538189312079030D+00 weight(7) = 0.182603415044923588866763667969D+00 weight(8) = 0.189450610455068496285396723208D+00 weight(9) = 0.189450610455068496285396723208D+00 weight(10) = 0.182603415044923588866763667969D+00 weight(11) = 0.169156519395002538189312079030D+00 weight(12) = 0.149595988816576732081501730547D+00 weight(13) = 0.124628971255533872052476282192D+00 weight(14) = 0.951585116824927848099251076022D-01 weight(15) = 0.622535239386478928628438369944D-01 weight(16) = 0.271524594117540948517805724560D-01 else if ( norder == 17 ) then xtab(1) = - 0.990575475314417335675434019941D+00 xtab(2) = - 0.950675521768767761222716957896D+00 xtab(3) = - 0.880239153726985902122955694488D+00 xtab(4) = - 0.781514003896801406925230055520D+00 xtab(5) = - 0.657671159216690765850302216643D+00 xtab(6) = - 0.512690537086476967886246568630D+00 xtab(7) = - 0.351231763453876315297185517095D+00 xtab(8) = - 0.178484181495847855850677493654D+00 xtab(9) = 0.0D+00 xtab(10) = 0.178484181495847855850677493654D+00 xtab(11) = 0.351231763453876315297185517095D+00 xtab(12) = 0.512690537086476967886246568630D+00 xtab(13) = 0.657671159216690765850302216643D+00 xtab(14) = 0.781514003896801406925230055520D+00 xtab(15) = 0.880239153726985902122955694488D+00 xtab(16) = 0.950675521768767761222716957896D+00 xtab(17) = 0.990575475314417335675434019941D+00 weight(1) = 0.241483028685479319601100262876D-01 weight(2) = 0.554595293739872011294401653582D-01 weight(3) = 0.850361483171791808835353701911D-01 weight(4) = 0.111883847193403971094788385626D+00 weight(5) = 0.135136368468525473286319981702D+00 weight(6) = 0.154045761076810288081431594802D+00 weight(7) = 0.168004102156450044509970663788D+00 weight(8) = 0.176562705366992646325270990113D+00 weight(9) = 0.179446470356206525458265644262D+00 weight(10) = 0.176562705366992646325270990113D+00 weight(11) = 0.168004102156450044509970663788D+00 weight(12) = 0.154045761076810288081431594802D+00 weight(13) = 0.135136368468525473286319981702D+00 weight(14) = 0.111883847193403971094788385626D+00 weight(15) = 0.850361483171791808835353701911D-01 weight(16) = 0.554595293739872011294401653582D-01 weight(17) = 0.241483028685479319601100262876D-01 else if ( norder == 18 ) then xtab(1) = - 0.991565168420930946730016004706D+00 xtab(2) = - 0.955823949571397755181195892930D+00 xtab(3) = - 0.892602466497555739206060591127D+00 xtab(4) = - 0.803704958972523115682417455015D+00 xtab(5) = - 0.691687043060353207874891081289D+00 xtab(6) = - 0.559770831073947534607871548525D+00 xtab(7) = - 0.411751161462842646035931793833D+00 xtab(8) = - 0.251886225691505509588972854878D+00 xtab(9) = - 0.847750130417353012422618529358D-01 xtab(10) = 0.847750130417353012422618529358D-01 xtab(11) = 0.251886225691505509588972854878D+00 xtab(12) = 0.411751161462842646035931793833D+00 xtab(13) = 0.559770831073947534607871548525D+00 xtab(14) = 0.691687043060353207874891081289D+00 xtab(15) = 0.803704958972523115682417455015D+00 xtab(16) = 0.892602466497555739206060591127D+00 xtab(17) = 0.955823949571397755181195892930D+00 xtab(18) = 0.991565168420930946730016004706D+00 weight(1) = 0.216160135264833103133427102665D-01 weight(2) = 0.497145488949697964533349462026D-01 weight(3) = 0.764257302548890565291296776166D-01 weight(4) = 0.100942044106287165562813984925D+00 weight(5) = 0.122555206711478460184519126800D+00 weight(6) = 0.140642914670650651204731303752D+00 weight(7) = 0.154684675126265244925418003836D+00 weight(8) = 0.164276483745832722986053776466D+00 weight(9) = 0.169142382963143591840656470135D+00 weight(10) = 0.169142382963143591840656470135D+00 weight(11) = 0.164276483745832722986053776466D+00 weight(12) = 0.154684675126265244925418003836D+00 weight(13) = 0.140642914670650651204731303752D+00 weight(14) = 0.122555206711478460184519126800D+00 weight(15) = 0.100942044106287165562813984925D+00 weight(16) = 0.764257302548890565291296776166D-01 weight(17) = 0.497145488949697964533349462026D-01 weight(18) = 0.216160135264833103133427102665D-01 else if ( norder == 19 ) then xtab(1) = - 0.992406843843584403189017670253D+00 xtab(2) = - 0.960208152134830030852778840688D+00 xtab(3) = - 0.903155903614817901642660928532D+00 xtab(4) = - 0.822714656537142824978922486713D+00 xtab(5) = - 0.720966177335229378617095860824D+00 xtab(6) = - 0.600545304661681023469638164946D+00 xtab(7) = - 0.464570741375960945717267148104D+00 xtab(8) = - 0.316564099963629831990117328850D+00 xtab(9) = - 0.160358645640225375868096115741D+00 xtab(10) = 0.0D+00 xtab(11) = 0.160358645640225375868096115741D+00 xtab(12) = 0.316564099963629831990117328850D+00 xtab(13) = 0.464570741375960945717267148104D+00 xtab(14) = 0.600545304661681023469638164946D+00 xtab(15) = 0.720966177335229378617095860824D+00 xtab(16) = 0.822714656537142824978922486713D+00 xtab(17) = 0.903155903614817901642660928532D+00 xtab(18) = 0.960208152134830030852778840688D+00 xtab(19) = 0.992406843843584403189017670253D+00 weight(1) = 0.194617882297264770363120414644D-01 weight(2) = 0.448142267656996003328381574020D-01 weight(3) = 0.690445427376412265807082580060D-01 weight(4) = 0.914900216224499994644620941238D-01 weight(5) = 0.111566645547333994716023901682D+00 weight(6) = 0.128753962539336227675515784857D+00 weight(7) = 0.142606702173606611775746109442D+00 weight(8) = 0.152766042065859666778855400898D+00 weight(9) = 0.158968843393954347649956439465D+00 weight(10) = 0.161054449848783695979163625321D+00 weight(11) = 0.158968843393954347649956439465D+00 weight(12) = 0.152766042065859666778855400898D+00 weight(13) = 0.142606702173606611775746109442D+00 weight(14) = 0.128753962539336227675515784857D+00 weight(15) = 0.111566645547333994716023901682D+00 weight(16) = 0.914900216224499994644620941238D-01 weight(17) = 0.690445427376412265807082580060D-01 weight(18) = 0.448142267656996003328381574020D-01 weight(19) = 0.194617882297264770363120414644D-01 else if ( norder == 20 ) then xtab(1) = - 0.993128599185094924786122388471D+00 xtab(2) = - 0.963971927277913791267666131197D+00 xtab(3) = - 0.912234428251325905867752441203D+00 xtab(4) = - 0.839116971822218823394529061702D+00 xtab(5) = - 0.746331906460150792614305070356D+00 xtab(6) = - 0.636053680726515025452836696226D+00 xtab(7) = - 0.510867001950827098004364050955D+00 xtab(8) = - 0.373706088715419560672548177025D+00 xtab(9) = - 0.227785851141645078080496195369D+00 xtab(10) = - 0.765265211334973337546404093988D-01 xtab(11) = 0.765265211334973337546404093988D-01 xtab(12) = 0.227785851141645078080496195369D+00 xtab(13) = 0.373706088715419560672548177025D+00 xtab(14) = 0.510867001950827098004364050955D+00 xtab(15) = 0.636053680726515025452836696226D+00 xtab(16) = 0.746331906460150792614305070356D+00 xtab(17) = 0.839116971822218823394529061702D+00 xtab(18) = 0.912234428251325905867752441203D+00 xtab(19) = 0.963971927277913791267666131197D+00 xtab(20) = 0.993128599185094924786122388471D+00 weight(1) = 0.176140071391521183118619623519D-01 weight(2) = 0.406014298003869413310399522749D-01 weight(3) = 0.626720483341090635695065351870D-01 weight(4) = 0.832767415767047487247581432220D-01 weight(5) = 0.101930119817240435036750135480D+00 weight(6) = 0.118194531961518417312377377711D+00 weight(7) = 0.131688638449176626898494499748D+00 weight(8) = 0.142096109318382051329298325067D+00 weight(9) = 0.149172986472603746787828737002D+00 weight(10) = 0.152753387130725850698084331955D+00 weight(11) = 0.152753387130725850698084331955D+00 weight(12) = 0.149172986472603746787828737002D+00 weight(13) = 0.142096109318382051329298325067D+00 weight(14) = 0.131688638449176626898494499748D+00 weight(15) = 0.118194531961518417312377377711D+00 weight(16) = 0.101930119817240435036750135480D+00 weight(17) = 0.832767415767047487247581432220D-01 weight(18) = 0.626720483341090635695065351870D-01 weight(19) = 0.406014298003869413310399522749D-01 weight(20) = 0.176140071391521183118619623519D-01 else if ( norder == 32 ) then xtab(1) = - 0.997263861849481563544981128665D+00 xtab(2) = - 0.985611511545268335400175044631D+00 xtab(3) = - 0.964762255587506430773811928118D+00 xtab(4) = - 0.934906075937739689170919134835D+00 xtab(5) = - 0.896321155766052123965307243719D+00 xtab(6) = - 0.849367613732569970133693004968D+00 xtab(7) = - 0.794483795967942406963097298970D+00 xtab(8) = - 0.732182118740289680387426665091D+00 xtab(9) = - 0.663044266930215200975115168663D+00 xtab(10) = - 0.587715757240762329040745476402D+00 xtab(11) = - 0.506899908932229390023747474378D+00 xtab(12) = - 0.421351276130635345364119436172D+00 xtab(13) = - 0.331868602282127649779916805730D+00 xtab(14) = - 0.239287362252137074544603209166D+00 xtab(15) = - 0.144471961582796493485186373599D+00 xtab(16) = - 0.483076656877383162348125704405D-01 xtab(17) = 0.483076656877383162348125704405D-01 xtab(18) = 0.144471961582796493485186373599D+00 xtab(19) = 0.239287362252137074544603209166D+00 xtab(20) = 0.331868602282127649779916805730D+00 xtab(21) = 0.421351276130635345364119436172D+00 xtab(22) = 0.506899908932229390023747474378D+00 xtab(23) = 0.587715757240762329040745476402D+00 xtab(24) = 0.663044266930215200975115168663D+00 xtab(25) = 0.732182118740289680387426665091D+00 xtab(26) = 0.794483795967942406963097298970D+00 xtab(27) = 0.849367613732569970133693004968D+00 xtab(28) = 0.896321155766052123965307243719D+00 xtab(29) = 0.934906075937739689170919134835D+00 xtab(30) = 0.964762255587506430773811928118D+00 xtab(31) = 0.985611511545268335400175044631D+00 xtab(32) = 0.997263861849481563544981128665D+00 weight(1) = 0.701861000947009660040706373885D-02 weight(2) = 0.162743947309056706051705622064D-01 weight(3) = 0.253920653092620594557525897892D-01 weight(4) = 0.342738629130214331026877322524D-01 weight(5) = 0.428358980222266806568786466061D-01 weight(6) = 0.509980592623761761961632446895D-01 weight(7) = 0.586840934785355471452836373002D-01 weight(8) = 0.658222227763618468376500637069D-01 weight(9) = 0.723457941088485062253993564785D-01 weight(10) = 0.781938957870703064717409188283D-01 weight(11) = 0.833119242269467552221990746043D-01 weight(12) = 0.876520930044038111427714627518D-01 weight(13) = 0.911738786957638847128685771116D-01 weight(14) = 0.938443990808045656391802376681D-01 weight(15) = 0.956387200792748594190820022041D-01 weight(16) = 0.965400885147278005667648300636D-01 weight(17) = 0.965400885147278005667648300636D-01 weight(18) = 0.956387200792748594190820022041D-01 weight(19) = 0.938443990808045656391802376681D-01 weight(20) = 0.911738786957638847128685771116D-01 weight(21) = 0.876520930044038111427714627518D-01 weight(22) = 0.833119242269467552221990746043D-01 weight(23) = 0.781938957870703064717409188283D-01 weight(24) = 0.723457941088485062253993564785D-01 weight(25) = 0.658222227763618468376500637069D-01 weight(26) = 0.586840934785355471452836373002D-01 weight(27) = 0.509980592623761761961632446895D-01 weight(28) = 0.428358980222266806568786466061D-01 weight(29) = 0.342738629130214331026877322524D-01 weight(30) = 0.253920653092620594557525897892D-01 weight(31) = 0.162743947309056706051705622064D-01 weight(32) = 0.701861000947009660040706373885D-02 else if ( norder == 64 ) then xtab(1) = - 0.999305041735772139456905624346D+00 xtab(2) = - 0.996340116771955279346924500676D+00 xtab(3) = - 0.991013371476744320739382383443D+00 xtab(4) = - 0.983336253884625956931299302157D+00 xtab(5) = - 0.973326827789910963741853507352D+00 xtab(6) = - 0.961008799652053718918614121897D+00 xtab(7) = - 0.946411374858402816062481491347D+00 xtab(8) = - 0.929569172131939575821490154559D+00 xtab(9) = - 0.910522137078502805756380668008D+00 xtab(10) = - 0.889315445995114105853404038273D+00 xtab(11) = - 0.865999398154092819760783385070D+00 xtab(12) = - 0.840629296252580362751691544696D+00 xtab(13) = - 0.813265315122797559741923338086D+00 xtab(14) = - 0.783972358943341407610220525214D+00 xtab(15) = - 0.752819907260531896611863774886D+00 xtab(16) = - 0.719881850171610826848940217832D+00 xtab(17) = - 0.685236313054233242563558371031D+00 xtab(18) = - 0.648965471254657339857761231993D+00 xtab(19) = - 0.611155355172393250248852971019D+00 xtab(20) = - 0.571895646202634034283878116659D+00 xtab(21) = - 0.531279464019894545658013903544D+00 xtab(22) = - 0.489403145707052957478526307022D+00 xtab(23) = - 0.446366017253464087984947714759D+00 xtab(24) = - 0.402270157963991603695766771260D+00 xtab(25) = - 0.357220158337668115950442615046D+00 xtab(26) = - 0.311322871990210956157512698560D+00 xtab(27) = - 0.264687162208767416373964172510D+00 xtab(28) = - 0.217423643740007084149648748989D+00 xtab(29) = - 0.169644420423992818037313629748D+00 xtab(30) = - 0.121462819296120554470376463492D+00 xtab(31) = - 0.729931217877990394495429419403D-01 xtab(32) = - 0.243502926634244325089558428537D-01 xtab(33) = 0.243502926634244325089558428537D-01 xtab(34) = 0.729931217877990394495429419403D-01 xtab(35) = 0.121462819296120554470376463492D+00 xtab(36) = 0.169644420423992818037313629748D+00 xtab(37) = 0.217423643740007084149648748989D+00 xtab(38) = 0.264687162208767416373964172510D+00 xtab(39) = 0.311322871990210956157512698560D+00 xtab(40) = 0.357220158337668115950442615046D+00 xtab(41) = 0.402270157963991603695766771260D+00 xtab(42) = 0.446366017253464087984947714759D+00 xtab(43) = 0.489403145707052957478526307022D+00 xtab(44) = 0.531279464019894545658013903544D+00 xtab(45) = 0.571895646202634034283878116659D+00 xtab(46) = 0.611155355172393250248852971019D+00 xtab(47) = 0.648965471254657339857761231993D+00 xtab(48) = 0.685236313054233242563558371031D+00 xtab(49) = 0.719881850171610826848940217832D+00 xtab(50) = 0.752819907260531896611863774886D+00 xtab(51) = 0.783972358943341407610220525214D+00 xtab(52) = 0.813265315122797559741923338086D+00 xtab(53) = 0.840629296252580362751691544696D+00 xtab(54) = 0.865999398154092819760783385070D+00 xtab(55) = 0.889315445995114105853404038273D+00 xtab(56) = 0.910522137078502805756380668008D+00 xtab(57) = 0.929569172131939575821490154559D+00 xtab(58) = 0.946411374858402816062481491347D+00 xtab(59) = 0.961008799652053718918614121897D+00 xtab(60) = 0.973326827789910963741853507352D+00 xtab(61) = 0.983336253884625956931299302157D+00 xtab(62) = 0.991013371476744320739382383443D+00 xtab(63) = 0.996340116771955279346924500676D+00 xtab(64) = 0.999305041735772139456905624346D+00 weight(1) = 0.178328072169643294729607914497D-02 weight(2) = 0.414703326056246763528753572855D-02 weight(3) = 0.650445796897836285611736039998D-02 weight(4) = 0.884675982636394772303091465973D-02 weight(5) = 0.111681394601311288185904930192D-01 weight(6) = 0.134630478967186425980607666860D-01 weight(7) = 0.157260304760247193219659952975D-01 weight(8) = 0.179517157756973430850453020011D-01 weight(9) = 0.201348231535302093723403167285D-01 weight(10) = 0.222701738083832541592983303842D-01 weight(11) = 0.243527025687108733381775504091D-01 weight(12) = 0.263774697150546586716917926252D-01 weight(13) = 0.283396726142594832275113052002D-01 weight(14) = 0.302346570724024788679740598195D-01 weight(15) = 0.320579283548515535854675043479D-01 weight(16) = 0.338051618371416093915654821107D-01 weight(17) = 0.354722132568823838106931467152D-01 weight(18) = 0.370551285402400460404151018096D-01 weight(19) = 0.385501531786156291289624969468D-01 weight(20) = 0.399537411327203413866569261283D-01 weight(21) = 0.412625632426235286101562974736D-01 weight(22) = 0.424735151236535890073397679088D-01 weight(23) = 0.435837245293234533768278609737D-01 weight(24) = 0.445905581637565630601347100309D-01 weight(25) = 0.454916279274181444797709969713D-01 weight(26) = 0.462847965813144172959532492323D-01 weight(27) = 0.469681828162100173253262857546D-01 weight(28) = 0.475401657148303086622822069442D-01 weight(29) = 0.479993885964583077281261798713D-01 weight(30) = 0.483447622348029571697695271580D-01 weight(31) = 0.485754674415034269347990667840D-01 weight(32) = 0.486909570091397203833653907347D-01 weight(33) = 0.486909570091397203833653907347D-01 weight(34) = 0.485754674415034269347990667840D-01 weight(35) = 0.483447622348029571697695271580D-01 weight(36) = 0.479993885964583077281261798713D-01 weight(37) = 0.475401657148303086622822069442D-01 weight(38) = 0.469681828162100173253262857546D-01 weight(39) = 0.462847965813144172959532492323D-01 weight(40) = 0.454916279274181444797709969713D-01 weight(41) = 0.445905581637565630601347100309D-01 weight(42) = 0.435837245293234533768278609737D-01 weight(43) = 0.424735151236535890073397679088D-01 weight(44) = 0.412625632426235286101562974736D-01 weight(45) = 0.399537411327203413866569261283D-01 weight(46) = 0.385501531786156291289624969468D-01 weight(47) = 0.370551285402400460404151018096D-01 weight(48) = 0.354722132568823838106931467152D-01 weight(49) = 0.338051618371416093915654821107D-01 weight(50) = 0.320579283548515535854675043479D-01 weight(51) = 0.302346570724024788679740598195D-01 weight(52) = 0.283396726142594832275113052002D-01 weight(53) = 0.263774697150546586716917926252D-01 weight(54) = 0.243527025687108733381775504091D-01 weight(55) = 0.222701738083832541592983303842D-01 weight(56) = 0.201348231535302093723403167285D-01 weight(57) = 0.179517157756973430850453020011D-01 weight(58) = 0.157260304760247193219659952975D-01 weight(59) = 0.134630478967186425980607666860D-01 weight(60) = 0.111681394601311288185904930192D-01 weight(61) = 0.884675982636394772303091465973D-02 weight(62) = 0.650445796897836285611736039998D-02 weight(63) = 0.414703326056246763528753572855D-02 weight(64) = 0.178328072169643294729607914497D-02 else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LEGENDRE_SET - Fatal error!' write ( *, '(a,i6)' ) ' Illegal value of NORDER = ', norder write ( *, '(a)' ) ' Legal values are 1 to 20, 32 or 64.' stop end if return end subroutine map ( code, w, lda, n ) ! !******************************************************************************* ! !! MAP returns the interpolation matrix for any available element. ! ! ! Formula: ! ! Given data Q(J) associated with the nodes, the coefficients of ! the interpolating polynomial are ! ! A(I) = W(I,J) * Q(J) ! ! In other words, if we let PHI(I)(R,S) be the I-th basis function, ! evaluated at the point (R,S), and we let REXP(I) and SEXP(I) ! be the exponents of R and S in the I-th associated polynomial, ! then the interpolating polynomial P(R,S) has two forms: ! ! P(R,S) = ! ! = Sum ( 1 <= I <= N ) Q(I) * PHI(I)(R,S) ! = Sum ( 1 <= I <= N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 09 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! ! Output, real W(LDA,N), the interpolation matrix. ! ! Input, integer LDA, the leading dimension of the W array. ! LDA must be at least N. ! ! Output, integer N, the order of the element. ! implicit none ! integer lda integer, parameter :: maxn = 20 integer n ! real area character ( len = * ) code integer i integer info integer ipivot(maxn) integer j real r(maxn) integer rexp(maxn) real rfact real s(maxn) integer sexp(maxn) real sfact real w(lda,n) real work(maxn) ! ! Get the (R,S) location of the nodes. ! call node ( code, n, r, s, area ) if ( lda < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MAP - Fatal error!' write ( *, '(a)' ) ' LDA < N.' stop end if if ( maxn < n ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MAP - Fatal error!' write ( *, '(a)' ) ' Internal parameter MAXN exceeded.' write ( *, '(a)' ) ' MAXN < N.' stop end if ! ! Get the associated polynomials. ! call poly ( code, n, rexp, sexp ) ! ! Set up the Vandermonde matrix. ! Factors of the form 0**0 are to be understood as 1. ! do i = 1, n do j = 1, n if ( rexp(j) == 0 ) then rfact = 1.0E+00 else rfact = r(i)**rexp(j) end if if ( sexp(j) == 0 ) then sfact = 1.0E+00 else sfact = s(i)**sexp(j) end if w(i,j) = rfact * sfact end do end do ! ! Factor the Vandermonde matrix. ! call sge_fa ( w, lda, n, ipivot, info ) if ( info /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MAP - Fatal error!' write ( *, '(a)' ) ' The Vandermonde matrix is singular.' stop end if ! ! Invert the Vandermonde matrix. ! call sge_inv ( w, lda, n, ipivot, work ) return end subroutine map_test ( code ) ! !******************************************************************************* ! !! MAP_TEST tests the map routines. ! ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = 4 ) CODE, the code for the element. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! implicit none ! integer, parameter :: maxn = 20 integer, parameter :: lda = maxn ! character ( len = * ) code integer order_code integer i integer j integer n real w(maxn,maxn) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MAP_TEST' write ( *, '(a,a)' ) ' The interpolation matrix for element ', trim ( code ) write ( *, '(a)' ) ' ' n = order_code ( code ) call map ( code, w, lda, n ) do i = 1, n write ( *, '(7f9.3)' ) w(i,1:n) end do return end subroutine node ( code, n, r, s, area ) ! !******************************************************************************* ! !! NODE returns the basis nodes for any available element. ! ! ! Modified: ! ! 04 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element desired. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(N), S(N), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area character ( len = * ) code integer n real r(*) real s(*) ! if ( code == 'Q4' ) then call node_q4 ( n, r, s, area ) else if ( code == 'Q8' ) then call node_q8 ( n, r, s, area ) else if ( code == 'Q9' ) then call node_q9 ( n, r, s, area ) else if ( code == 'Q12' ) then call node_q12 ( n, r, s, area ) else if ( code == 'Q16' ) then call node_q16 ( n, r, s, area ) else if ( code == 'QL' ) then call node_ql ( n, r, s, area ) else if ( code == 'T3' ) then call node_t3 ( n, r, s, area ) else if ( code == 'T6' ) then call node_t6 ( n, r, s, area ) else if ( code == 'T10' ) then call node_t10 ( n, r, s, area ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'NODE - Fatal error!' write ( *, '(a,a)' ) ' Illegal value of CODE = ', trim ( code ) stop end if return end subroutine node_q4 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_Q4 returns the basis nodes for a 4 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 3-------4 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1-------2 ! | ! +--0---R---1--> ! ! Modified: ! ! 19 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(4), S(4), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(4) real s(4) ! n = 4 r(1:4) = (/ 0.0E+00, 1.0E+00, 0.0E+00, 1.0E+00 /) s(1:4) = (/ 0.0E+00, 0.0E+00, 1.0E+00, 1.0E+00 /) area = 1.0E+00 return end subroutine node_q8 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_Q8 returns the basis nodes for an 8 node quadrilateral. ! ! ! Comment: ! ! This element is known as the quadratic "serendipity" element. ! ! Diagram: ! ! | ! 1 6---7---8 ! | | | ! | | | ! S 4 5 ! | | | ! | | | ! 0 1---2---3 ! | ! +--0---R---1--> ! ! Modified: ! ! 19 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(8), S(8), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(8) real s(8) ! n = 8 r(1:8) = (/ 0.0E+00, 0.5E+00, 1.0E+00, 0.0E+00, & 1.0E+00, 0.0E+00, 0.5E+00, 1.0E+00 /) s(1:8) = (/ 0.0E+00, 0.0E+00, 0.0E+00, 0.5E+00, & 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00 /) area = 1.0E+00 return end subroutine node_q9 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_Q9 returns the basis nodes for a 9 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 7---8---9 ! | | : | ! | | : | ! S 4...5...6 ! | | : | ! | | : | ! 0 1---2---3 ! | ! +--0--R--1--> ! ! Modified: ! ! 19 September 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Input, real R(9), S(9), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(9) real s(9) ! n = 9 r(1:9) = (/ 0.0E+00, 0.5E+00, 1.0E+00, 0.0E+00, 0.5E+00, & 1.0E+00, 0.0E+00, 0.5E+00, 1.0E+00 /) s(1:9) = (/ 0.0E+00, 0.0E+00, 0.0E+00, 0.5E+00, 0.5E+00, & 0.5E+00, 1.0E+00, 1.0E+00, 1.0E+00 /) area = 1.0E+00 return end subroutine node_q12 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_Q12 returns the basis nodes for a 12 node quadrilateral. ! ! ! Comment: ! ! This element is known as the cubic "serendipity" element. ! ! Diagram: ! ! | ! 1 9-10-11-12 ! | | | ! | 7 8 ! S | | ! | 5 6 ! | | | ! 0 1--2--3--4 ! | ! +--0---R---1--> ! ! Modified: ! ! 04 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(12), S(12), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real, parameter :: a = 0.0E+00 real area real, parameter :: b = 1.0E+00 / 3.0E+00 real, parameter :: c = 2.0E+00 / 3.0E+00 real, parameter :: d = 1.0E+00 integer n real r(12) real s(12) ! n = 12 r(1:12) = (/ a, b, c, d, a, d, a, d, a, b, c, d /) s(1:12) = (/ a, a, a, a, b, b, c, c, d, d, d, d /) area = 1.0E+00 return end subroutine node_q16 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_Q16 returns the basis nodes for a 16 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 13--14--15--16 ! | | : : | ! | | : : | ! | 9..10..11..12 ! S | : : | ! | | : : | ! | 5...6...7...8 ! | | : : | ! | | : : | ! 0 1---2---3---4 ! | ! +--0-----R-----1--> ! ! Modified: ! ! 14 October 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(16), S(16), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer i integer j integer k integer n real r(16) real s(16) ! n = 16 k = 0 do i = 0, 3 do j = 0, 3 k = k + 1 r(k) = real ( j ) / 3.0E+00 s(k) = real ( i ) / 3.0E+00 end do end do area = 1.0E+00 return end subroutine node_ql ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_QL returns the basis nodes for a quadratic/linear. ! ! ! Diagram: ! ! | ! 1 4---5---6 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1---2---3 ! | ! +--0---R---1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(6), S(6), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(6) real s(6) ! n = 6 r(1:6) = (/ 0.0E+00, 0.5E+00, 1.0E+00, 0.0E+00, 0.5E+00, 1.0E+00 /) s(1:6) = (/ 0.0E+00, 0.0E+00, 0.0E+00, 1.0E+00, 1.0E+00, 1.0E+00 /) area = 1.0E+00 return end subroutine node_t3 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_T3 returns the basis nodes for the 3 node triangle. ! ! ! Diagram: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(3), S(3), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(3) real s(3) n = 3 r(1:3) = (/ 0.0E+00, 1.0E+00, 0.0E+00 /) s(1:3) = (/ 0.0E+00, 0.0E+00, 1.0E+00 /) area = 0.5E+00 return end subroutine node_t6 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_T6 returns the basis nodes for a 6 node triangle. ! ! ! Diagram: ! ! | ! 1 6 ! | |\ ! | | \ ! S 4 5 ! | | \ ! | | \ ! 0 1--2--3 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(6), S(6), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(6) real s(6) ! n = 6 r(1:6) = (/ 0.0E+00, 0.5E+00, 1.0E+00, 0.0E+00, 0.5E+00, 0.0E+00 /) s(1:6) = (/ 0.0E+00, 0.0E+00, 0.0E+00, 0.5E+00, 0.5E+00, 1.0E+00 /) area = 0.5E+00 return end subroutine node_t10 ( n, r, s, area ) ! !******************************************************************************* ! !! NODE_T10 returns the basis nodes for a 10 node triangle. ! ! ! Diagram: ! ! | ! 1 10 ! | |\ ! | | \ ! | 8 9 ! | | \ ! S | \ ! | 5 6 7 ! | | \ ! | | \ ! 0 1--2--3--4 ! | ! +--0----R---1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of nodes in the element. ! ! Output, real R(10), S(10), the coordinates of the basis nodes. ! ! Output, real AREA, the area of the element. ! implicit none ! real area integer n real r(10) real s(10) ! n = 10 r(1) = 0.0E+00 s(1) = 0.0E+00 r(2) = 1.0E+00 / 3.0E+00 s(2) = 0.0E+00 r(3) = 2.0E+00 / 3.0E+00 s(3) = 0.0E+00 r(4) = 1.0E+00 s(4) = 0.0E+00 r(5) = 0.0E+00 s(5) = 1.0E+00 / 3.0E+00 r(6) = 1.0E+00 / 3.0E+00 s(6) = 1.0E+00 / 3.0E+00 r(7) = 2.0E+00 / 3.0E+00 s(7) = 1.0E+00 / 3.0E+00 r(8) = 0.0E+00 s(8) = 2.0E+00 / 3.0E+00 r(9) = 1.0E+00 / 3.0E+00 s(9) = 2.0E+00 / 3.0E+00 r(10) = 0.0E+00 s(10) = 1.0E+00 area = 0.5E+00 return end function order_code ( code ) ! !******************************************************************************* ! !! ORDER_CODE returns the order for each element. ! ! ! List: ! ! CODE Order Definition ! ---- ----- ---------- ! Q4 4 4 node linear Lagrange/serendipity quadrilateral; ! Q8 8 8 node quadratic serendipity quadrilateral; ! Q9 9 9 node quadratic Lagrange quadrilateral; ! Q12 12 12 node cubic serendipity quadrilateral; ! Q16 16 16 node cubic Lagrange quadrilateral; ! QL 6 6 node linear/quadratic quadrilateral; ! T3 3 3 node linear triangle; ! T6 6 6 node quadratic triangle; ! T10 10 10 node cubic triangle. ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, the code for the element. ! ! Output, integer ORDER_CODE, the order of the element. ! implicit none ! character ( len = * ) code integer order_code ! if ( code == 'Q4' ) then order_code = 4 else if ( code == 'Q8' ) then order_code = 8 else if ( code == 'Q9' ) then order_code = 9 else if ( code == 'Q12' ) then order_code = 12 else if ( code == 'Q16' ) then order_code = 16 else if ( code == 'QL' ) then order_code = 6 else if ( code == 'T3' ) then order_code = 3 else if ( code == 'T6' ) then order_code = 6 else if ( code == 'T10' ) then order_code = 10 else order_code = 0 end if return end subroutine poly ( code, n, rexp, sexp ) ! !******************************************************************************* ! !! POLY returns the polynomials associated with any available element. ! ! ! Modified: ! ! 04 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element desired. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', 'T3', ! 'T6' and 'T10'. ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! character ( len = * ) code integer n integer rexp(*) integer sexp(*) ! if ( code == 'Q4' ) then call poly_q4 ( n, rexp, sexp ) else if ( code == 'Q8' ) then call poly_q8 ( n, rexp, sexp ) else if ( code == 'Q9' ) then call poly_q9 ( n, rexp, sexp ) else if ( code == 'Q12' ) then call poly_q12 ( n, rexp, sexp ) else if ( code == 'Q16' ) then call poly_q16 ( n, rexp, sexp ) else if ( code == 'QL' ) then call poly_ql ( n, rexp, sexp ) else if ( code == 'T3' ) then call poly_t3 ( n, rexp, sexp ) else if ( code == 'T6' ) then call poly_t6 ( n, rexp, sexp ) else if ( code == 'T10' ) then call poly_t10 ( n, rexp, sexp ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'POLY - Fatal error!' write ( *, '(a,a)' ) ' Illegal value of CODE = ', trim ( code ) stop end if return end subroutine poly_q4 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_Q4 returns the polynomials associated with a 4 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 3-----4 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1-----2 ! | ! +--0--R--1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 4 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:4) = (/ 0, 0, 1, 1 /) sexp(1:4) = (/ 0, 1, 0, 1 /) return end subroutine poly_q8 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_Q8 returns the polynomials associated with an 8 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 6---7---8 ! | | | ! | | | ! S 4 5 ! | | | ! | | | ! 0 1---2---3 ! | ! +--0--R--1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 8 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:8) = (/ 0, 0, 0, 1, 1, 1, 2, 2 /) sexp(1:8) = (/ 0, 1, 2, 0, 1, 2, 0, 2 /) return end subroutine poly_q9 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_Q9 returns the polynomials associated with a 9 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 7---8---9 ! | | : | ! | | : | ! S 4...5...6 ! | | : | ! | | : | ! 0 1---2---3 ! | ! +--0--R--1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 9 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:9) = (/ 0, 0, 0, 1, 1, 1, 2, 2, 2 /) sexp(1:9) = (/ 0, 1, 2, 0, 1, 2, 0, 1, 2 /) return end subroutine poly_q12 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_Q12 returns the polynomials associated with a 12 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 9--10--11--12 ! | | | ! | | | ! | 7 8 ! S | | ! | | | ! | 5 6 ! | | | ! | | | ! 0 1---2---3---4 ! | ! +--0-----R-----1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 12 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:12) = (/ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 3, 3 /) sexp(1:12) = (/ 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 0, 1 /) return end subroutine poly_q16 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_Q16 returns the polynomials associated with a 16 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 13--14--15--16 ! | | : : | ! | | : : | ! | 9..10..11..12 ! S | : : | ! | | : : | ! | 5...6...7...8 ! | | : : | ! | | : : | ! 0 1---2---3---4 ! | ! +--0-----R-----1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 16 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:16) = (/ 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /) sexp(1:16) = (/ 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3 /) return end subroutine poly_ql ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_QL returns the polynomials for a quadratic/linear quadrilateral. ! ! ! Diagram: ! ! | ! 1 4---5---6 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1---2---3 ! | ! +--0---R---1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 6 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:6) = (/ 0, 0, 1, 1, 2, 2 /) sexp(1:6) = (/ 0, 1, 0, 1, 0, 1 /) return end subroutine poly_t3 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_T3 returns the polynomials associated with a 3 node triangle. ! ! ! Diagram: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 3 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:3) = (/ 0, 0, 1 /) sexp(1:3) = (/ 0, 1, 0 /) return end subroutine poly_t6 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_T6 returns the polynomials associated with a 6 node triangle. ! ! ! Diagram: ! ! | ! 1 6 ! | |\ ! | | \ ! S 4 5 ! | | \ ! | | \ ! 0 1--2--3 ! | ! +--0--R--1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 6 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:6) = (/ 0, 0, 0, 1, 1, 2 /) sexp(1:6) = (/ 0, 1, 2, 0, 1, 0 /) return end subroutine poly_t10 ( n, rexp, sexp ) ! !******************************************************************************* ! !! POLY_T10 returns the polynomials associated with a 10 node triangle. ! ! ! Diagram: ! ! | ! 1 10 ! | |\ ! | | \ ! | 8 9 ! | | \ ! S | \ ! | 5 6 7 ! | | \ ! | | \ ! 0 1--2--3--4 ! | ! +--0----R---1--> ! ! Formula: ! ! Given coefficients A(I), the polynomial interpolant at (R,S) is ! ! P(R,S) = SUM ( I = 1 to N ) A(I) * R**REXP(I) * S**SEXP(I) ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer N, the number of polynomials. ! ! Output, integer REXP(N), SEXP(N), the powers of R and S associated ! with each polynomial. ! implicit none ! integer, parameter :: n_internal = 10 ! integer n integer rexp(n_internal) integer sexp(n_internal) ! n = n_internal rexp(1:10) = (/ 0, 0, 0, 0, 1, 1, 1, 2, 2, 3 /) sexp(1:10) = (/ 0, 1, 2, 3, 0, 1, 2, 0, 1, 0 /) return end subroutine r_swap ( x, y ) ! !******************************************************************************* ! !! R_SWAP switches two real values. ! ! ! Modified: ! ! 30 November 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real X, Y. On output, the values of X and ! Y have been interchanged. ! implicit none ! real x real y real z ! z = x x = y y = z return end subroutine s_l2norm ( psi_num, element_num, quad_num, element_area, & quad_weight, psi_quad, s_coef, l2norm ) ! !******************************************************************************* ! !! S_L2NORM computes the "big" L2 norm of a scalar function over a region. ! ! ! Discussion: ! ! It is assumed that a set of finite element basis functions PSI ! have been defined over a collection of elements that compose ! the region. Moreover, integrals over the region are assumed to ! be approximated by applying a fixed quadrature rule to all the ! elements. ! ! Finally, we assume that we have a scalar function S(X), which ! is represented as a linear combination of basis functions, and ! it is desired to determine the L2 norm of S. ! ! This routine estimates the integral ! ! Sqrt ( Integral ( X in Omega ) S(X)**2 dOmega ) ! ! using the finite element representation of S, and applying the ! given quadrature rule. ! ! It assumes that a (probably very large) array of data is available, ! recording the value of each basis function PSI in every element ! at every quadrature point. If this is true, then the computation ! becomes very simple. ! ! If your problem is small or sufficient memory is available, this ! may be an efficient computation. It requires that the value of ! all the basis functions be stored for all the elements and all ! the quadrature points. That particular information need only ! be computed once. ! ! Actually, no assumptions are made here about the dimension of the ! space, so this same code can handle problems in 1D, 2D, 3D and ! so on. ! ! Modified: ! ! 30 March 2002 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer PSI_NUM, the number of global element functions. ! ! Input, integer ELEMENT_NUM, the number of elements. ! ! Input, integer QUAD_NUM, the number of points in the quadrature rule. ! ! Input, real ELEMENT_AREA(ELEMENT_NUM), the area of each element. ! ! Input, real QUAD_WEIGHT(QUAD_NUM), the quadrature weights associated ! with the quadrature points. ! ! Input, real PSI_QUAD(PSI_NUM,ELEMENT_NUM,QUAD_NUM), the value of ! the I-th PSI function in the J-th element at the K-th quadrature ! point. ! ! Input, real S_COEF(PSI_NUM), the coefficients of the PSI functions ! associated with the scalar function S. ! ! Output, L2NORM, the L2 norm of the scalar function S over the region. ! implicit none ! integer element_num integer psi_num integer quad_num ! real element_area(element_num) integer j integer k real l2norm real psi_quad(psi_num,element_num,quad_num) real quad_weight(quad_num) real s(element_num,quad_num) real s_coef(psi_num) real t(quad_num) real u ! ! #1: Sum over all basis functions to get the value of S in each element ! at each quadrature point. ! ! The MATMUL function requires that one of its arguments be shaped ! like a vector, and one like a 2 dimensional matrix, so we have ! to insert a loop on the quadrature points. ! do j = 1, quad_num s(1:element_num,j) = matmul ( & s_coef(1:psi_num), psi_quad(1:psi_num,1:element_num,j) ) end do ! ! #2: Sum over all elements to get the value of S**2 weighted by its element ! area. SUM expects to see vector quantities, so we have a loop on ! quadrature points. ! do k = 1, quad_num t(k) = sum ( s(1:element_num,k)**2 * element_area(1:element_num) ) end do ! ! #3: Sum over all quadrature points weighted by the quadrature weight. ! u = dot_product ( t(1:quad_num), quad_weight(1:quad_num) ) l2norm = sqrt ( u ) return end subroutine serene ( type, ve, vn, vne, vnw, vs, vse, vsw, vw, vterp ) ! !******************************************************************************* ! !! SERENE interpolates data using a serendipity quadrilateral. ! ! ! Modified: ! ! 02 March 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = 2 ) TYPE, tells SERENE the geometry of the ! finite element that surrounds the point of interest. The options ! are displayed in the following table, which suggests the meaning ! of each option by its position: ! ! | | ! NW * N * NE ! | | ! -*-*-*-*-*- ! | | ! W * C * E ! | | ! -*-*-*-*-*- ! | | ! SW * S * SE ! | | ! ! Input, real VE, VN, VNE, VNW, VS, VSE, VSW, VW, ! are the values of the function at the nodes to the east, ! north, northeast, northwest, south, southeast, southwest and ! west of the point of interest. If the finite element is of ! type 'C', then all 8 values are needed. However, if the ! finite element is of type 'SE', for instance, then only three ! values are needed, namely VE, VN, and VNW, since these are ! the only node positions defined in such a finite element. ! ! Output, real VTERP, the interpolated value of the ! function at the point of interest. ! implicit none ! real eta real pe real pn real pne real pnw real ps real pse real psw real pw character ( len = 2 ) type real ve real vn real vne real vnw real vs real vse real vsw real vw real vterp real xsi ! ! To make this routine more general, simply pass in the values of XSI ! and ETA at which the interpolated value is desired. ! ! By setting XSI = ETA = 0, we are asking for the interpolated value ! at the center of the finite element. ! xsi = 0.0E+00 eta = 0.0E+00 ! ! 8 node center ! ! Polynomial space is spanned by: ! ! 1 ! x y ! x^2 xy y^2 ! x^2y xy^2 ! ! ! ^ 1 4--7--3 ! | ! ! ! E ! ! ! T 0 8 X 6 ! A ! ! ! | ! ! ! V -1 1--5--2 ! ! -1 0 1 ! ! <---XSI---> ! if ( type == 'C' ) then psw = - 0.25E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 - eta ) & * ( 1.0E+00 + xsi + eta ) pse = - 0.25E+00 * ( 1.0E+00 + xsi ) * ( 1.0E+00 - eta ) & * ( 1.0E+00 - xsi + eta ) pne = - 0.25E+00 * ( 1.0E+00 + xsi ) * ( 1.0E+00 + eta ) & * ( 1.0E+00 - xsi - eta ) pnw = - 0.25E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 + eta ) & * ( 1.0E+00 + xsi - eta ) ps = 0.50E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 + xsi ) & * ( 1.0E+00 - eta ) pe = 0.50E+00 * ( 1.0E+00 + xsi ) * ( 1.0E+00 + eta ) & * ( 1.0E+00 - eta ) pn = 0.50E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 + xsi ) & * ( 1.0E+00 + eta ) pw = 0.50E+00 * ( 1.0E+00 - xsi ) * ( 1.0E+00 + eta ) & * ( 1.0E+00 - eta ) ! ! 5 node side ! ! ^ 1 ! | ! E ! T 0 8 X 6 ! A ! ! ! | ! ! ! V -1 1--5--2 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'N' ) then psw = 0.5E+00 * ( xsi - 1.0E+00 ) * ( 1.0E+00 + xsi + eta ) pse = -0.5E+00 * ( xsi + 1.0E+00 ) * ( 1.0E+00 - xsi + eta ) ps = - ( xsi + 1.0E+00 ) * ( xsi - 1.0E+00 ) pe = 0.5E+00 * ( xsi + 1.0E+00 ) * ( eta + 1.0E+00 ) pw = -0.5E+00 * ( xsi - 1.0E+00 ) * ( eta + 1.0E+00 ) ! ! ^ 1 4--7 ! | ! ! E ! ! T 0 8 X ! A ! ! | ! ! V -1 1--5 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'E' ) then pse = 0.5E+00 * ( eta - 1.0E+00 ) * ( 1.0E+00 + xsi + eta ) pne = -0.5E+00 * ( eta + 1.0E+00 ) * ( 1.0E+00 + xsi - eta ) ps = -0.5E+00 * ( xsi + 1.0E+00 ) * ( eta - 1.0E+00 ) pn = 0.5E+00 * ( xsi + 1.0E+00 ) * ( eta + 1.0E+00 ) pw = - ( eta + 1.0E+00 ) * ( eta - 1.0E+00 ) ! ! 5 node side ! ! ^ 1 7--3 ! | ! ! E ! ! T 0 X 6 ! A ! ! | ! ! V -1 5--2 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'W' ) then pse = 0.5E+00 * ( eta - 1.0E+00 ) * ( 1.0E+00 - xsi + eta ) pne = - 0.5E+00 * ( eta + 1.0E+00 ) * ( 1.0E+00 - xsi - eta ) ps = 0.5E+00 * ( xsi - 1.0E+00 ) * ( eta - 1.0E+00 ) pe = - ( eta - 1.0E+00 ) * ( eta + 1.0E+00 ) pn = - 0.5E+00 * ( xsi - 1.0E+00 ) * ( eta + 1.0E+00 ) ! ! 5 node side ! ! ^ 1 4--7--3 ! | ! ! ! E ! ! ! T 0 8 X 6 ! A ! | ! V -1 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'S' ) then pne = - 0.5E+00 * ( xsi + 1.0E+00 ) * ( 1.0E+00 - xsi - eta ) pnw = 0.5E+00 * ( xsi - 1.0E+00 ) * ( 1.0E+00 + xsi - eta ) pe = - 0.5E+00 * ( eta - 1.0E+00 ) * ( xsi + 1.0E+00 ) pn = - ( xsi + 1.0E+00 ) * ( xsi - 1.0E+00 ) pw = 0.5E+00 * ( eta - 1.0E+00 ) * ( xsi - 1.0E+00 ) ! ! 3 node corner ! ! Polynomial space is spanned by: ! ! 1 ! x y ! ! ! ^ 1 ! | ! E ! T 0 8 X ! A ! ! | ! ! V -1 1--5 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'NE' ) then psw = - 1.0E+00 - xsi - eta ps = 1.0E+00 + xsi pw = 1.0E+00 + eta ! ! 3 node corner ! ! Polynomial space is spanned by: ! ! 1 ! x y ! ! ^ 1 ! | ! E ! T 0 X 6 ! A ! ! | ! ! V -1 5--2 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'NW' ) then pse = 1.0E+00 + xsi - eta ps = 1.0E+00 - xsi pe = 1.0E+00 + eta ! ! 3 node corner ! ! Polynomial space is spanned by: ! 1 ! x y ! ! ! ^ 1 4--7 ! | ! ! E ! ! T 0 8 X ! A ! | ! V -1 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'SE' ) then pnw = - 1.0E+00 - xsi + eta pn = 1.0E+00 + xsi pw = 1.0E+00 - eta ! ! 3 node corner ! ! Polynomial space is spanned by: ! ! 1 ! x y ! ! ^ 1 7--3 ! | ! ! E ! ! T 0 X 6 ! A ! | ! V -1 ! ! -1 0 1 ! ! <---XSI---> ! else if ( type == 'SW' ) then pne = - 1.0E+00 + xsi + eta pe = 1.0E+00 - eta pn = 1.0E+00 - xsi end if vterp = vsw * psw + vse * pse + vne * pne + vnw * pnw & + vs * ps + ve * pe + vn * pn + vw * pw return end subroutine sge_check ( lda, m, n, ierror ) ! !******************************************************************************* ! !! SGE_CHECK checks the dimensions of a general matrix. ! ! ! Modified: ! ! 11 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least M. ! ! Input, integer M, the number of rows of the matrix. ! M must be positive. ! ! Input, integer N, the number of columns of the matrix. ! N must be positive. ! ! Output, integer IERROR, reports whether any errors were detected. ! IERROR is set to 0 before the checks are made, and then: ! IERROR = IERROR + 1 if LDA is illegal; ! IERROR = IERROR + 2 if M is illegal; ! IERROR = IERROR + 4 if N is illegal. ! implicit none ! integer ierror integer lda integer m integer n ! ierror = 0 if ( lda < m ) then ierror = ierror + 1 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGE_CHECK - Illegal LDA = ', lda end if if ( m < 1 ) then ierror = ierror + 2 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGE_CHECK - Illegal M = ', m end if if ( n < 1 ) then ierror = ierror + 4 write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) 'SGE_CHECK - Illegal N = ', n end if return end subroutine sge_fa ( a, lda, n, pivot, info ) ! !******************************************************************************* ! !! SGE_FA factors a general matrix. ! ! ! Discussion: ! ! SGE_FA is a simplified version of the LINPACK routine SGEFA. ! ! Modified: ! ! 26 February 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(LDA,N), the matrix to be factored. ! On output, A contains an upper triangular matrix and the multipliers ! which were used to obtain it. The factorization can be written ! A = L * U, where L is a product of permutation and unit lower ! triangular matrices and U is upper triangular. ! ! Input, integer LDA, the leading dimension of the array. ! LDA must be at least N. ! ! Input, integer N, the order of the matrix. ! N must be positive. ! ! Output, integer PIVOT(N), a vector of pivot indices. ! ! Output, integer INFO, singularity flag. ! 0, no singularity detected. ! nonzero, the factorization failed on the INFO-th step. ! implicit none ! integer lda integer n ! real a(lda,n) integer i integer ierror integer info integer pivot(n) integer j integer k integer l real t ! ! Check the dimensions. ! call sge_check ( lda, n, n, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_FA - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions.' return end if info = 0 do k = 1, n-1 ! ! Find L, the index of the pivot row. ! l = k do i = k+1, n if ( abs ( a(i,k) ) > abs ( a(l,k) ) ) then l = i end if end do pivot(k) = l ! ! If the pivot index is zero, the algorithm has failed. ! if ( a(l,k) == 0.0E+00 ) then info = k write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_FA - Fatal error!' write ( *, '(a,i6)' ) ' Zero pivot on step ', info return end if ! ! Interchange rows L and K if necessary. ! if ( l /= k ) then call r_swap ( a(l,k), a(k,k) ) end if ! ! Normalize the values that lie below the pivot entry A(K,K). ! a(k+1:n,k) = -a(k+1:n,k) / a(k,k) ! ! Row elimination with column indexing. ! do j = k+1, n if ( l /= k ) then call r_swap ( a(l,j), a(k,j) ) end if a(k+1:n,j) = a(k+1:n,j) + a(k+1:n,k) * a(k,j) end do end do pivot(n) = n if ( a(n,n) == 0.0E+00 ) then info = n write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_FA - Fatal error!' write ( *, '(a,i6)' ) ' Zero pivot on step ', info end if return end subroutine sge_inv ( a, lda, n, pivot ) ! !******************************************************************************* ! !! SGE_INV computes the inverse of a matrix factored by SGE_FA. ! ! ! Note: ! ! SGE_INV is a simplified standalone version of the LINPACK routine ! SGEDI. ! ! Modified: ! ! 16 September 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, real A(LDA,N). ! On input, the factor information computed by SGE_FA. ! On output, the inverse matrix. ! ! Input, integer LDA, the leading dimension of the array A, ! which must be at least N. ! ! Input, integer N, the order of the matrix A. ! ! Input, integer PIVOT(N), the pivot vector from SGE_FA. ! implicit none ! integer lda integer n ! real a(lda,n) integer i integer ierror integer pivot(n) integer j integer k real temp real work(n) ! ! Check the dimensions. ! call sge_check ( lda, n, n, ierror ) if ( ierror /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SGE_INV - Fatal error!' write ( *, '(a)' ) ' Illegal dimensions.' return end if ! ! Compute Inverse(U). ! do k = 1, n a(k,k) = 1.0E+00 / a(k,k) a(1:k-1,k) = -a(1:k-1,k) * a(k,k) do j = k + 1, n temp = a(k,j) a(k,j) = 0.0E+00 a(1:k,j) = a(1:k,j) + temp * a(1:k,k) end do end do ! ! Form Inverse(U) * Inverse(L). ! do k = n - 1, 1, -1 work(k+1:n) = a(k+1:n,k) a(k+1:n,k) = 0.0E+00 do j = k + 1, n a(1:n,k) = a(1:n,k) + work(j) * a(1:n,j) end do if ( pivot(k) /= k ) then do i = 1, n call r_swap ( a(i,k), a(i,pivot(k)) ) end do end if end do return end subroutine shape ( code, n, r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE evaluates shape functions for any available element. ! ! ! Modified: ! ! 10 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! ! Input, integer N, the number of nodes in the element. ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(N), the basis functions at the point. ! ! Output, real DTDR(N), the R basis derivatives at the point. ! ! Output, real DTDS(N), the S basis derivatives at the point. ! implicit none ! integer n ! character ( len = * ) code real dtdr(n) real dtds(n) real r real s real t(n) ! if ( code == 'Q4' ) then call shape_q4 ( r, s, t, dtdr, dtds ) else if( code == 'Q8' ) then call shape_q8 ( r, s, t, dtdr, dtds ) else if ( code == 'Q9' ) then call shape_q9 ( r, s, t, dtdr, dtds ) else if ( code == 'Q12' ) then call shape_q12 ( r, s, t, dtdr, dtds ) else if ( code == 'Q16' ) then call shape_q16 ( r, s, t, dtdr, dtds ) else if ( code == 'QL' ) then call shape_ql ( r, s, t, dtdr, dtds ) else if ( code == 'T3' ) then call shape_t3 ( r, s, t, dtdr, dtds ) else if ( code == 'T6' ) then call shape_t6 ( r, s, t, dtdr, dtds ) else if ( code == 'T10' ) then call shape_t10 ( r, s, t, dtdr, dtds ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SHAPE - Fatal error!' write ( *, '(a,a)' ) ' Unrecognized code = ', trim ( code ) stop end if return end subroutine shape_test ( code ) ! !******************************************************************************* ! !! SHAPE_TEST verifies the shape function values at the basis nodes. ! ! ! Modified: ! ! 09 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) CODE, identifies the element to be used. ! Legal values include 'Q4', 'Q8', 'Q9', 'Q12', 'Q16', 'QL', ! 'T3', 'T6' and 'T10'. ! implicit none ! integer, parameter :: maxn = 16 ! real area character ( len = * ) code real dtdr(maxn) real dtds(maxn) integer i integer j integer n real r(maxn) real rsum real s(maxn) real ssum real t(maxn) ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'SHAPE_TEST' write ( *, '(a,a)' ) ' Verify shape functions of type ', trim ( code ) call node ( code, n, r, s, area ) write ( *, '(a)' ) ' ' write ( *, '(a,i6)' ) ' Number of nodes = ', n write ( *, '(a)' ) ' Basis function values at basis nodes' write ( *, '(a)' ) ' should form the identity matrix.' write ( *, '(a)' ) ' ' do i = 1, n call shape ( code, n, r(i), s(i), t, dtdr, dtds ) write ( *, '(10f7.3)' ) t(1:n) end do write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' R and S derivatives should sum to 0.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' dTdR sum, dTdS sum:' write ( *, '(a)' ) ' ' do i = 1, n call shape ( code, n, r(i), s(i), t, dtdr, dtds ) rsum = sum ( dtdr(1:n) ) ssum = sum ( dtds(1:n) ) write ( *, '(2f14.8)' ) rsum, ssum end do return end subroutine shape_q4 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_Q4 evaluates shape functions for a 4 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 3-----4 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(4), the basis functions at the point. ! ! Output, real DTDR(4), the R basis derivatives at the point. ! ! Output, real DTDS(4), the S basis derivatives at the point. ! implicit none ! real dtdr(4) real dtds(4) real r real s real t(4) ! t(1) = ( 1.0E+00 - r ) * ( 1.0E+00 - s ) t(2) = r * ( 1.0E+00 - s ) t(3) = ( 1.0E+00 - r ) * s t(4) = r * s dtdr(1) = - 1.0E+00 + s dtdr(2) = 1.0E+00 - s dtdr(3) = - s dtdr(4) = s dtds(1) = - 1.0E+00 + r dtds(2) = - r dtds(3) = 1.0E+00 - r dtds(4) = r return end subroutine shape_q8 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_Q8 evaluates shape functions for an 8 node quadrilateral. ! ! ! Comment: ! ! This element is known as the "serendipity" element. ! ! Diagram: ! ! | ! 1 6--7--8 ! | | | ! | | | ! S 4 5 ! | | | ! | | | ! 0 1--2--3 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(8), the basis functions at the point. ! ! Output, real DTDR(8), the R basis derivatives at the point. ! ! Output, real DTDS(8), the S basis derivatives at the point. ! implicit none ! real dtdr(8) real dtds(8) real r real s real t(8) ! t(1) = ( r - 1.0E+00 ) * ( s - 1.0E+00 ) & * ( 1.0E+00 - 2.0E+00 * r - 2.0E+00 * s ) t(2) = 4.0E+00 * r * ( r - 1.0E+00 ) * ( s - 1.0E+00 ) t(3) = r * ( s - 1.0E+00 ) & * ( 1.0E+00 - 2.0E+00 * r + 2.0E+00 * s ) t(4) = 4.0E+00 * ( r - 1.0E+00 ) * s * ( s - 1.0E+00 ) t(5) = - 4.0E+00 * r * s * ( s - 1.0E+00 ) t(6) = ( r - 1.0E+00 ) * s & * ( 2.0E+00 * r - 2.0E+00 * s + 1.0E+00 ) t(7) = - 4.0E+00 * r * ( r - 1.0E+00 ) * s t(8) = r * s & * ( 2.0E+00 * r + 2.0E+00 * s - 3.0E+00 ) dtdr(1) = ( s - 1.0E+00 ) * ( - 4.0E+00 * r - 2.0E+00 * s + 3.0E+00 ) dtdr(2) = 4.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) * ( s - 1.0E+00 ) dtdr(3) = ( s - 1.0E+00 ) * ( - 4.0E+00 * r + 2.0E+00 * s + 1.0E+00 ) dtdr(4) = 4.0E+00 * s * ( s - 1.0E+00 ) dtdr(5) = - 4.0E+00 * s * ( s - 1.0E+00 ) dtdr(6) = s * ( 4.0E+00 * r - 2.0E+00 * s - 1.0E+00 ) dtdr(7) = - 4.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) * s dtdr(8) = s * ( 4.0E+00 * r + 2.0E+00 * s - 3.0E+00 ) dtds(1) = ( r - 1.0E+00 ) * ( - 4.0E+00 * s - 2.0E+00 * r + 3.0E+00 ) dtds(2) = 4.0E+00 * r * ( r - 1.0E+00 ) dtds(3) = r * ( 4.0E+00 * s - 2.0E+00 * r - 1.0E+00 ) dtds(4) = 4.0E+00 * ( r - 1.0E+00 ) * ( 2.0E+00 * s - 1.0E+00 ) dtds(5) = - 4.0E+00 * r * ( 2.0E+00 * s - 1.0E+00 ) dtds(6) = ( r - 1.0E+00 ) * ( - 4.0E+00 * s + 2.0E+00 * r + 1.0E+00 ) dtds(7) = - 4.0E+00 * r * ( r - 1.0E+00 ) dtds(8) = r * ( 4.0E+00 * s + 2.0E+00 * r - 3.0E+00 ) return end subroutine shape_q9 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_Q9 evaluates shape functions for a 9 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 7--8--9 ! | | : | ! | | : | ! S 4..5..6 ! | | : | ! | | : | ! 0 1--2--3 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(9), the basis functions at the point. ! ! Output, real DTDR(9), the R basis derivatives at the point. ! ! Output, real DTDS(9), the S basis derivatives at the point. ! implicit none ! real dtdr(9) real dtds(9) real r real s real t(9) ! t(1) = 4.0E+00 * ( r - 1.0E+00 ) * ( r - 0.5E+00 ) * ( s - 1.0E+00 ) & * ( s - 0.5 ) t(2) = - 8.0E+00 * r * ( r - 1.0E+00 ) * ( s - 1.0E+00 ) * ( s - 0.5E+00 ) t(3) = 4.0E+00 * r * ( r - 0.5E+00 ) * ( s - 1.0E+00 ) * ( s - 0.5E+00 ) t(4) = - 8.0E+00 * ( r - 1.0E+00 ) * ( r - 0.5E+00 ) * s * ( s - 1.0E+00 ) t(5) = 16.0E+00 * r * ( r - 1.0E+00 ) * s * ( s - 1.0E+00 ) t(6) = - 8.0E+00 * r * ( r - 0.5E+00 ) * s * ( s - 1.0E+00 ) t(7) = 4.0E+00 * ( r - 1.0E+00 ) * ( r - 0.5E+00 ) * s * ( s - 0.5E+00 ) t(8) = - 8.0E+00 * r * ( r - 1.0E+00 ) * s * ( s - 0.5E+00 ) t(9) = 4.0E+00 * r * ( r - 0.5E+00 ) * s * ( s - 0.5E+00 ) dtdr(1) = 4.0E+00 * ( 2.0E+00 * r - 1.5E+00 ) * ( s - 1.0E+00 ) & * ( s - 0.5E+00 ) dtdr(2) = - 8.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) * ( s - 1.0E+00 ) & * ( s - 0.5E+00 ) dtdr(3) = 4.0E+00 * ( 2.0E+00 * r - 0.5E+00 ) * ( s - 1.0E+00 ) & * ( s - 0.5E+00 ) dtdr(4) = - 8.0E+00 * ( 2.0E+00 * r - 1.5E+00 ) * s * ( s - 1.0E+00 ) dtdr(5) = 16.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) * s * ( s - 1.0E+00 ) dtdr(6) = - 8.0E+00 * ( 2.0E+00 * r - 0.5E+00 ) * s * ( s - 1.0E+00 ) dtdr(7) = 4.0E+00 * ( 2.0E+00 * r - 1.5E+00 ) * s * ( s - 0.5E+00 ) dtdr(8) = - 8.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) * s * ( s - 0.5E+00 ) dtdr(9) = 4.0E+00 * ( 2.0E+00 * r - 0.5E+00 ) * s * ( s - 0.5E+00 ) dtds(1) = 4.0E+00 * ( r - 1.0E+00 ) * ( r - 0.5E+00 ) & * ( 2.0E+00 * s - 1.5E+00 ) dtds(2) = - 8.0E+00 * r * ( r - 1.0E+00 ) * ( 2.0E+00 * s - 1.5E+00 ) dtds(3) = 4.0E+00 * r * ( r - 0.5E+00 ) * ( 2.0E+00 * s - 1.5 ) dtds(4) = - 8.0E+00 * ( r - 1.0E+00 ) * ( r - 0.5E+00 ) & * ( 2.0E+00 * s - 1.0E+00 ) dtds(5) = 16.0E+00 * r * ( r - 1.0E+00 ) * ( 2.0E+00 * s - 1.0E+00 ) dtds(6) = - 8.0E+00 * r * ( r - 0.5E+00 ) * ( 2.0E+00 * s - 1.0E+00 ) dtds(7) = 4.0E+00 * ( r - 1.0E+00 ) * ( r - 0.5E+00 ) & * ( 2.0E+00 * s - 0.5E+00 ) dtds(8) = - 8.0E+00 * r * ( r - 1.0E+00 ) * ( 2.0E+00 * s - 0.5E+00 ) dtds(9) = 4.0E+00 * r * ( r - 0.5E+00 ) * ( 2.0E+00 * s - 0.5E+00 ) return end subroutine shape_q12 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_Q12 evaluates shape functions for a 12 node quadrilateral. ! ! ! Note: ! ! This routine is being worked on. ! ! Diagram: ! ! | ! 1 9-10-11-12 ! | | | ! | 7 8 ! S | | ! | 5 6 ! | | | ! 0 1--2--3--4 ! | ! +--0---R---1--> ! ! Modified: ! ! 12 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(12), the basis functions at the point. ! ! Output, real DTDR(12), the R basis derivatives at the point. ! ! Output, real DTDS(12), the S basis derivatives at the point. ! implicit none ! real a real b real c real corner real d real dcdr real dcds real dtdr(12) real dtds(12) real r real s real t(12) ! a = 0.0E+00 b = 1.0E+00 / 3.0E+00 c = 2.0E+00 / 3.0E+00 d = 1.0E+00 corner = 9.0E+00 * ( ( 2.0E+00 * r - 1.0E+00 )**2 & + ( 2.0E+00 * s - 1.0E+00 )**2 ) - 10.0E+00 t(1) = 0.125E+00 * ( r - d ) * ( s - d ) * corner t(2) = - 13.5E+00 * ( r - a ) * ( r - c ) * ( r - d ) * ( s - d ) t(3) = 13.5E+00 * ( r - a ) * ( r - b ) * ( r - d ) * ( s - d ) t(4) = - 0.125E+00 * ( r - a ) * ( s - d ) * corner t(5) = - 13.5E+00 * ( r - d ) * ( s - a ) * ( s - c ) * ( s - d ) t(6) = 13.5E+00 * ( r - a ) * ( s - a ) * ( s - c ) * ( s - d ) t(7) = 13.5E+00 * ( r - d ) * ( s - a ) * ( s - b ) * ( s - d ) t(8) = - 13.5E+00 * ( r - a ) * ( s - a ) * ( s - b ) * ( s - d ) t(9) = - 0.125E+00 * ( r - d ) * ( s - a ) * corner t(10) = 13.5E+00 * ( r - a ) * ( r - c ) * ( r - d ) * ( s - a ) t(11) = - 13.5E+00 * ( r - a ) * ( r - b ) * ( r - d ) * ( s - a ) t(12) = 0.125E+00 * ( r - a ) * ( s - a ) * corner dcdr = 36.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) dtdr(1) = 0.125 * ( s - d ) * ( ( r - d ) * dcdr + corner ) dtdr(2) = - 13.5E+00 * ( s - d ) * ( 3.0E+00 * r**2 & - 2.0E+00 * ( a + c + d ) * r + a * c + c * d + d * a ) dtdr(3) = 13.5E+00 * ( s - d ) * ( 3.0E+00 * r**2 & - 2.0E+00 * ( a + b + d ) * r + a * b + b * d + d * a ) dtdr(4) = - 0.125E+00 * ( s - d ) * ( ( r - a ) * dcdr + corner ) dtdr(5) = - 13.5E+00 * ( s - a ) * ( s - c ) * ( s - d ) dtdr(6) = 13.5E+00 * ( s - a ) * ( s - c ) * ( s - d ) dtdr(7) = 13.5E+00 * ( s - a ) * ( s - b ) * ( s - d ) dtdr(8) = - 13.5E+00 * ( s - a ) * ( s - b ) * ( s - d ) dtdr(9) = - 0.125E+00 * ( s - a ) * ( ( r - d ) * dcdr + corner ) dtdr(10) = 13.5E+00 * ( s - a ) * ( 3.0E+00 * r**2 & - 2.0E+00 * ( a + c + d ) * r + a * c + c * d + d * a ) dtdr(11) = - 13.5E+00 * ( s - a ) * ( 3.0E+00 * r**2 & - 2.0E+00 * ( a + b + d ) * r + a * b + b * d + d * a ) dtdr(12) = 0.125E+00 * ( s - a ) * ( ( r - a ) * dcdr + corner ) dcds = 36.0E+00 * ( 2.0E+00 * s - 1.0E+00 ) dtds(1) = 0.125E+00 * ( r - d ) * ( corner + ( s - d ) * dcds ) dtds(2) = - 13.5E+00 * ( r - a ) * ( r - c ) * ( r - d ) dtds(3) = 13.5E+00 * ( r - a ) * ( r - b ) * ( r - d ) dtds(4) = - 0.125E+00 * ( r - a ) * ( corner + ( s - d ) * dcds ) dtds(5) = - 13.5E+00 * ( r - d ) * ( 3.0E+00 * s**2 & - 2.0E+00 * ( a + c + d ) * s + a * c + c * d + d * a ) dtds(6) = 13.5E+00 * ( r - a ) * ( 3.0E+00 * s**2 & - 2.0E+00 * ( a + c + d ) * s + a * c + c * d + d * a ) dtds(7) = 13.5E+00 * ( r - d ) * ( 3.0E+00 * s**2 & - 2.0E+00 * ( a + b + d ) * s + a * b + b * d + d * a ) dtds(8) = - 13.5E+00 * ( r - a ) * ( 3.0E+00 * s**2 & - 2.0E+00 * ( a + b + d ) * s + a * b + b * d + d * a ) dtds(9) = - 0.125E+00 * ( r - d ) * ( corner + ( s - a ) * dcds ) dtds(10) = 13.5E+00 * ( r - a ) * ( r - c ) * ( r - d ) dtds(11) = - 13.5E+00 * ( r - a ) * ( r - b ) * ( r - d ) dtds(12) = 0.125E+00 * ( r - a ) * ( corner + ( s - a ) * dcds ) return end subroutine shape_q16 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_Q16 evaluates shape functions for a 16 node quadrilateral. ! ! ! Diagram: ! ! | ! 1 13--14--15--16 ! | | : : | ! | | : : | ! | 9..10..11..12 ! S | : : | ! | | : : | ! | 5...6...7...8 ! | | : : | ! | | : : | ! 0 1---2---3---4 ! | ! +--0-----R-----1--> ! ! Modified: ! ! 12 December 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(16), the basis functions at the point. ! ! Output, real DTDR(16), the R basis derivatives at the point. ! ! Output, real DTDS(16), the S basis derivatives at the point. ! implicit none ! real dabc real dabd real dacd real dbcd real dtdr(16) real dtds(16) real r real ra real rb real rc real rd real s real sa real sb real sc real sd real t(16) ! ra = r - 0.0E+00 rb = r - 1.0E+00 / 3.0E+00 rc = r - 2.0E+00 / 3.0E+00 rd = r - 1.0E+00 sa = s - 0.0E+00 sb = s - 1.0E+00 / 3.0E+00 sc = s - 2.0E+00 / 3.0E+00 sd = s - 1.0E+00 t(1) = ( 81.0E+00 / 4.0E+00 ) * rb * rc * rd * sb * sc * sd t(2) = - ( 243.0E+00 / 4.0E+00 ) * ra * rc * rd * sb * sc * sd t(3) = ( 243.0E+00 / 4.0E+00 ) * ra * rb * rd * sb * sc * sd t(4) = - ( 81.0E+00 / 4.0E+00 ) * ra * rb * rc * sb * sc * sd t(5) = - ( 243.0E+00 / 4.0E+00 ) * rb * rc * rd * sa * sc * sd t(6) = ( 729.0E+00 / 4.0E+00 ) * ra * rc * rd * sa * sc * sd t(7) = - ( 729.0E+00 / 4.0E+00 ) * ra * rb * rd * sa * sc * sd t(8) = ( 243.0E+00 / 4.0E+00 ) * ra * rb * rc * sa * sc * sd t(9) = ( 243.0E+00 / 4.0E+00 ) * rb * rc * rd * sa * sb * sd t(10) = - ( 729.0E+00 / 4.0E+00 ) * ra * rc * rd * sa * sb * sd t(11) = ( 729.0E+00 / 4.0E+00 ) * ra * rb * rd * sa * sb * sd t(12) = - ( 243.0E+00 / 4.0E+00 ) * ra * rb * rc * sa * sb * sd t(13) = - ( 81.0E+00 / 4.0E+00 ) * rb * rc * rd * sa * sb * sc t(14) = ( 243.0E+00 / 4.0E+00 ) * ra * rc * rd * sa * sb * sc t(15) = - ( 243.0E+00 / 4.0E+00 ) * ra * rb * rd * sa * sb * sc t(16) = ( 81.0E+00 / 4.0E+00 ) * ra * rb * rc * sa * sb * sc dbcd = 3.0E+00 * r**2 - 4.0E+00 * r + 11.0E+00 / 9.0E+00 dacd = 3.0E+00 * r**2 - 10.0E+00 * r / 3.0E+00 + 2.0E+00 / 3.0E+00 dabd = 3.0E+00 * r**2 - 8.0E+00 * r / 3.0E+00 + 1.0E+00 / 3.0E+00 dabc = 3.0E+00 * r**2 - 2.0E+00 * r + 2.0E+00 / 9.0E+00 dtdr( 1) = ( 81.0E+00 / 4.0E+00 ) * dbcd * sb * sc * sd dtdr( 2) = - ( 243.0E+00 / 4.0E+00 ) * dacd * sb * sc * sd dtdr( 3) = ( 243.0E+00 / 4.0E+00 ) * dabd * sb * sc * sd dtdr( 4) = - ( 81.0E+00 / 4.0E+00 ) * dabc * sb * sc * sd dtdr( 5) = - ( 243.0E+00 / 4.0E+00 ) * dbcd * sa * sc * sd dtdr( 6) = ( 729.0E+00 / 4.0E+00 ) * dacd * sa * sc * sd dtdr( 7) = - ( 729.0E+00 / 4.0E+00 ) * dabd * sa * sc * sd dtdr( 8) = ( 243.0E+00 / 4.0E+00 ) * dabc * sa * sc * sd dtdr( 9) = ( 243.0E+00 / 4.0E+00 ) * dbcd * sa * sb * sd dtdr(10) = - ( 729.0E+00 / 4.0E+00 ) * dacd * sa * sb * sd dtdr(11) = ( 729.0E+00 / 4.0E+00 ) * dabd * sa * sb * sd dtdr(12) = - ( 243.0E+00 / 4.0E+00 ) * dabc * sa * sb * sd dtdr(13) = - ( 81.0E+00 / 4.0E+00 ) * dbcd * sa * sb * sc dtdr(14) = ( 243.0E+00 / 4.0E+00 ) * dacd * sa * sb * sc dtdr(15) = - ( 243.0E+00 / 4.0E+00 ) * dabd * sa * sb * sc dtdr(16) = ( 81.0E+00 / 4.0E+00 ) * dabc * sa * sb * sc dbcd = 3.0E+00 * s**2 - 4.0E+00 * s + 11.0E+00 / 9.0E+00 dacd = 3.0E+00 * s**2 - 10.0E+00 * s / 3.0E+00 + 2.0E+00 / 3.0E+00 dabd = 3.0E+00 * s**2 - 8.0E+00 * s / 3.0E+00 + 1.0E+00 / 3.0E+00 dabc = 3.0E+00 * s**2 - 2.0E+00 * s + 2.0E+00 / 9.0E+00 dtds( 1) = ( 81.0E+00 / 4.0E+00 ) * rb * rc * rd * dbcd dtds( 2) = - ( 243.0E+00 / 4.0E+00 ) * ra * rc * rd * dbcd dtds( 3) = ( 243.0E+00 / 4.0E+00 ) * ra * rb * rd * dbcd dtds( 4) = - ( 81.0E+00 / 4.0E+00 ) * ra * rb * rc * dbcd dtds( 5) = - ( 243.0E+00 / 4.0E+00 ) * rb * rc * rd * dacd dtds( 6) = ( 729.0E+00 / 4.0E+00 ) * ra * rc * rd * dacd dtds( 7) = - ( 729.0E+00 / 4.0E+00 ) * ra * rb * rd * dacd dtds( 8) = ( 243.0E+00 / 4.0E+00 ) * ra * rb * rc * dacd dtds( 9) = ( 243.0E+00 / 4.0E+00 ) * rb * rc * rd * dabd dtds(10) = - ( 729.0E+00 / 4.0E+00 ) * ra * rc * rd * dabd dtds(11) = ( 729.0E+00 / 4.0E+00 ) * ra * rb * rd * dabd dtds(12) = - ( 243.0E+00 / 4.0E+00 ) * ra * rb * rc * dabd dtds(13) = - ( 81.0E+00 / 4.0E+00 ) * rb * rc * rd * dabc dtds(14) = ( 243.0E+00 / 4.0E+00 ) * ra * rc * rd * dabc dtds(15) = - ( 243.0E+00 / 4.0E+00 ) * ra * rb * rd * dabc dtds(16) = ( 81.0E+00 / 4.0E+00 ) * ra * rb * rc * dabc return end subroutine shape_ql ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_QL evaluates shape functions for a 6 node quadratic/linear. ! ! ! Diagram: ! ! | ! 1 4--5--6 ! | | | ! | | | ! S | | ! | | | ! | | | ! 0 1--2--3 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(6), the basis functions at the point. ! ! Output, real DTDR(6), the R basis derivatives at the point. ! ! Output, real DTDS(6), the S basis derivatives at the point. ! implicit none ! real dtdr(6) real dtds(6) real r real s real t(6) ! t(1) = - 2.0E+00 * ( r - 0.5E+00 ) * ( r - 1.0E+00 ) * ( s - 1.0E+00 ) t(2) = 4.0E+00 * r * ( r - 1.0E+00 ) * ( s - 1.0E+00 ) t(3) = - 2.0E+00 * r * ( r - 0.5E+00 ) * ( s - 1.0E+00 ) t(4) = 2.0E+00 * ( r - 0.5E+00 ) * ( r - 1.0E+00 ) * s t(5) = - 4.0E+00 * r * ( r - 1.0E+00 ) * s t(6) = 2.0E+00 * r * ( r - 0.5E+00 ) * s dtdr(1) = 2.0E+00 * ( - 2.0E+00 * r + 1.5E+00 ) * ( s - 1.0E+00 ) dtdr(2) = 4.0E+00 * ( 2.0E+00 * r - 1.0E+00 ) * ( s - 1.0E+00 ) dtdr(3) = 2.0E+00 * ( - 2.0E+00 * r + 0.5E+00 ) * ( s - 1.0E+00 ) dtdr(4) = 2.0E+00 * ( 2.0E+00 * r - 1.5E+00 ) * s dtdr(5) = 4.0E+00 * ( - 2.0E+00 * r + 1.0E+00 ) * s dtdr(6) = 2.0E+00 * ( 2.0E+00 * r - 0.5E+00 ) * s dtds(1) = - 2.0E+00 * ( r - 0.5E+00 ) * ( r - 1.0E+00 ) dtds(2) = 4.0E+00 * r * ( r - 1.0E+00 ) dtds(3) = - 2.0E+00 * r * ( r - 0.5E+00 ) dtds(4) = 2.0E+00 * ( r - 0.5E+00 ) * ( r - 1.0E+00 ) dtds(5) = - 4.0E+00 * r * ( r - 1.0E+00 ) dtds(6) = 2.0E+00 * r * ( r - 0.5E+00 ) return end subroutine shape_t3 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_T3 evaluates shape functions for a 3 node triangle. ! ! ! Diagram: ! ! | ! 1 3 ! | |\ ! | | \ ! S | \ ! | | \ ! | | \ ! 0 1-----2 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(3), the basis functions at the point. ! ! Output, real DTDR(3), the R basis derivatives at the point. ! ! Output, real DTDS(3), the S basis derivatives at the point. ! implicit none ! real dtdr(3) real dtds(3) real r real s real t(3) ! t(1) = 1.0E+00 - r - s t(2) = r t(3) = s dtdr(1) = -1.0E+00 dtdr(2) = 1.0E+00 dtdr(3) = 0.0E+00 dtds(1) = -1.0E+00 dtds(2) = 0.0E+00 dtds(3) = 1.0E+00 return end subroutine shape_t6 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_T6 evaluates shape functions for a 6 node triangle. ! ! ! Diagram: ! ! | ! 1 6 ! | |\ ! | | \ ! S 4 5 ! | | \ ! | | \ ! 0 1--2--3 ! | ! +--0--R--1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(6), the basis functions at the point. ! ! Output, real DTDR(6), the R basis derivatives at the point. ! ! Output, real DTDS(6), the S basis derivatives at the point. ! implicit none ! real dtdr(6) real dtds(6) real r real s real t(6) ! t(1) = 2.0E+00 * ( 1.0E+00 - r - s ) * ( 0.5E+00 - r - s ) t(2) = 4.0E+00 * r * ( 1.0E+00 - r - s ) t(3) = 2.0E+00 * r * ( r - 0.5E+00 ) t(4) = 4.0E+00 * s * ( 1.0E+00 - r - s ) t(5) = 4.0E+00 * r * s t(6) = 2.0E+00 * s * ( s - 0.5E+00 ) dtdr(1) = - 3.0E+00 + 4.0E+00 * r + 4.0E+00 * s dtdr(2) = 4.0E+00 - 8.0E+00 * r - 4.0E+00 * s dtdr(3) = - 1.0E+00 + 4.0E+00 * r dtdr(4) = - 4.0E+00 * s dtdr(5) = 4.0E+00 * s dtdr(6) = 0.0E+00 dtds(1) = - 3.0E+00 + 4.0E+00 * r + 4.0E+00 * s dtds(2) = - 4.0E+00 * r dtds(3) = 0.0E+00 dtds(4) = 4.0E+00 - 4.0E+00 * r - 8.0E+00 * s dtds(5) = 4.0E+00 * r dtds(6) = - 1.0E+00 + 4.0E+00 * s return end subroutine shape_t10 ( r, s, t, dtdr, dtds ) ! !******************************************************************************* ! !! SHAPE_T10 evaluates shape functions for a 10 node triangle. ! ! ! Diagram: ! ! | ! 1 10 ! | |\ ! | | \ ! | 8 9 ! | | \ ! S | \ ! | 5 6 7 ! | | \ ! | | \ ! 0 1--2--3--4 ! | ! +--0----R---1--> ! ! Modified: ! ! 05 January 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, real R, S, the reference coordinates of a point. ! ! Output, real T(10), the basis functions at the point. ! ! Output, real DTDR(10), the R basis derivatives at the point. ! ! Output, real DTDS(10), the S basis derivatives at the point. ! implicit none ! real a real b real c real dtdr(10) real dtds(10) real r real s real t(10) ! a = 1.0E+00 / 3.0E+00 b = 2.0E+00 / 3.0E+00 c = 1.0E+00 t( 1) = 4.5E+00 * ( a - r - s ) * ( b - r - s ) * ( c - r - s ) t( 2) = 13.5E+00 * r * ( b - r - s ) * ( c - r - s ) t( 3) = - 13.5E+00 * r * ( a - r ) * ( c - r - s ) t( 4) = 4.5E+00 * r * ( a - r ) * ( b - r ) t( 5) = 13.5E+00 * s * ( b - r - s ) * ( c - r - s ) t( 6) = 27.0E+00 * r * s * ( c - r - s ) t( 7) = 13.5E+00 * r * s * ( r - a ) t( 8) = 13.5E+00 * s * ( s - a ) * ( c - r - s ) t( 9) = 13.5E+00 * r * s * ( s - a ) t(10) = 4.5E+00 * s * ( a - s ) * ( b - s ) dtdr( 1) = 4.5E+00 * ( ( a - s ) * ( 2.0E+00 * r - c - b + 2.0E+00 * s ) & - ( s - b ) * ( s - c ) - 2.0E+00 * ( 2.0E+00 * s - b - c ) * r & - 3.0E+00 * r**2 ) dtdr( 2) = 13.5E+00 * ( & ( s - b ) * ( s - c ) + 2.0E+00 * ( 2.0E+00 * s - b - c ) * r & + 3.0E+00 * r**2 ) dtdr( 3) = - 13.5E+00 * ( a * ( c - s ) + 2.0E+00 * ( s - a - c ) * r & + 3.0E+00 * r**2 ) dtdr( 4) = 4.5E+00 * ( a * b - 2.0E+00 * ( a + b ) * r + 3.0E+00 * r**2 ) dtdr( 5) = 13.5E+00 * s * ( 2.0E+00 * s - b - c + 2.0E+00 * r ) dtdr( 6) = 27.0E+00 * s * ( c - s - 2.0E+00 * r ) dtdr( 7) = 13.5E+00 * s * ( 2.0E+00 * r - a ) dtdr( 8) = - 13.5E+00 * s * ( s - a ) dtdr( 9) = 13.5E+00 * s * ( s - a) dtdr(10) = 0.0E+00 dtds( 1) = 4.5E+00 * ( ( a - r ) * ( 2.0E+00 * s - c - b + 2.0E+00 * r ) & - ( r - b ) * ( r - c ) - 2.0E+00 * ( 2.0E+00 * r - b - c ) * s & - 3.0E+00 * s**2 ) dtds( 2) = 13.5E+00 * r * ( 2.0E+00 * s + 2.0E+00 * r - b - c ) dtds( 3) = 13.5E+00 * r * ( a - r ) dtds( 4) = 0.0E+00 dtds( 5) = 13.5E+00 * ( ( r - b ) * ( r - c ) + & 2.0E+00 * ( 2.0E+00 * r - b - c ) * s + 3.0E+00 * s**2 ) dtds( 6) = 27.0E+00 * r * ( c - r - 2.0E+00 * s ) dtds( 7) = 13.5E+00 * r * ( r - a ) dtds( 8) = - 13.5E+00 * ( a * ( c - r ) + 2.0E+00 * ( r - c - a ) * s & + 3.0E+00 * s**2 ) dtds( 9) = 13.5E+00 * r * ( 2.0E+00 * s - a) dtds(10) = 4.5E+00 * ( a * b - 2.0E+00 * ( a + b ) * s + 3.0E+00 * s**2 ) return end subroutine timestamp ( ) ! !******************************************************************************* ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! ! Example: ! ! May 31 2001 9:45:54.872 AM ! ! Modified: ! ! 31 May 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none ! character ( len = 8 ) ampm integer d character ( len = 8 ) date integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s character ( len = 10 ) time integer values(8) integer y character ( len = 5 ) zone ! call date_and_time ( date, time, zone, values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(a,1x,i2,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & trim ( month(m) ), d, y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end