xref: /phasta/phSolver/common/genint.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine genint
c
c----------------------------------------------------------------------
c
c This subroutine inputs the integration information.
c
c----------------------------------------------------------------------
c
      include "common.h"

      real*8, allocatable :: tmpQpt (:,:), tmpQwt (:)
      real*8, allocatable :: tmpQptb(:,:), tmpQwtb(:)

c
c.... compute the shape function parameters
c


c
c.... get quadrature data for interior and boundary elements
c
c
c  Tets
c
      if (nen.eq.4) then
         nshape  = (ipord+1)*(ipord+2)*(ipord+3)/6
         nshapeb = (ipord+1)*(ipord+2)/2
      endif

      select case (intg(1,1))
      case (1)
         nint(1) = 1
      case (2)
         nint(1) = 4
      case (3)
         nint(1) = 16
      case (4)
         nint(1) = 29
      end select
      select case (intg(2,1))
      case (1)
         nintb(1) = 1
      case (2)
         nintb(1) = 3
      case (3)
         nintb(1) = 6
      case (4)
         nintb(1) = 12
      end select

      allocate (tmpQpt (4,nint(1)))
      allocate (tmpQwt (nint(1)))
      allocate (tmpQptb(4,nintb(1)))
      allocate (tmpQwtb(nintb(1)))

      call symtet(nint(1),tmpQpt,tmpQwt,nerr) ! interior elements
      Qpt(1,1:4,1:nint(1)) = tmpQpt(1:4,1:nint(1))
      Qwt(1,1:nint(1))     = tmpQwt(1:nint(1))


      call symtri(nintb(1),tmpQptb,tmpQwtb,nerr) ! boundary elements
      Qptb(1,1:4,1:nintb(1)) = tmpQptb(1:4,1:nintb(1))
      Qwtb(1,1:nintb(1))     = tmpQwtb(1:nintb(1))

      deallocate (tmpQpt)
      deallocate (tmpQwt)
      deallocate (tmpQptb)
      deallocate (tmpQwtb)

c
c.... adjust quadrature weights to be consistent with the
c     design of tau.
c
      Qwt(1,:) = (four/three)*Qwt(1,:)
      Qwtb(1,:) = two*Qwtb(1,:)
c
c     Hexes now
c
      if (nen.eq.8) then
         nshape = nen
         if ( ipord .gt. 1 ) then
            nshape = nshape + 12*(ipord - 1)
         endif

         if ( ipord .gt. 3 ) then
            nshape = nshape + 3*(ipord -2)*(ipord - 3)
         endif

         if ( ipord .gt. 5 ) then
            nshape = nshape + (ipord - 3)*(ipord - 4)*(ipord - 5)/6
         endif

         nshapeb = nenb

         if ( ipord .gt. 1 ) then
            nshapeb = nshapeb + 4*(ipord - 1)
         endif

         if ( ipord .gt. 3 ) then
            nshapeb = nshapeb +(ipord -2)*(ipord - 3)/2
         endif
      endif

      select case (intg(1,1))
      case (1)
         nint(2) = 1
      case (2)
         nint(2) = 8
      case (3)
         nint(2) =27
      case (4)
         nint(2) = 64
      case (5)
         nint(2) = 125
      end select

      select case (intg(2,1))
      case (1)
         nintb(2) = 1
      case (2)
         nintb(2) = 4
      case(3)
         nintb(2) = 9
      case (4)
         nintb(2) = 16
      case (5)
         nintb(2) = 25
      end select


      allocate (tmpQpt (4,nint(2)))
      allocate (tmpQwt (nint(2)))
      allocate (tmpQptb(4,nintb(2)))
      allocate (tmpQwtb(nintb(2)))

      call symhex(nint(2),tmpQpt,tmpQwt,nerr)
      Qpt(2,1:4,1:nint(2)) = tmpQpt(1:4,1:nint(2))
      Qwt(2,1:nint(2)) = tmpQwt(1:nint(2))

      call symquad(nintb(2),tmpQptb,tmpQwtb,nerr)
      Qptb(2,1:4,1:nintb(2)) = tmpQptb(1:4,1:nintb(2))
      Qwtb(2,1:nintb(2)) = tmpQwtb(1:nintb(2))

      deallocate (tmpQpt)
      deallocate (tmpQwt)
      deallocate (tmpQptb)
      deallocate (tmpQwtb)

      if (nen.eq.5) then        ! for pyramid

         nshape = nen           ! counting the nodal functions

         if ( ipord .gt. 1 ) then
            nshape = nshape + 8*(ipord - 1) ! counting the edge functions
         endif

         if ( ipord .gt. 2 ) then
                                ! counting the triangular face functions
            nshape = nshape + 2*(ipord -1)*(ipord - 2)
         endif

         if ( ipord .gt. 3 ) then
                                ! counting the quadrilateral face functions
            nshape = nshape + (ipord - 2)*(ipord - 3)/2
         endif

         if ( ipord .gt. 5 ) then
            nshape = nshape + (ipord - 3)*(ipord - 4)*(ipord - 5)/6
         endif

         nshapeb = nenb
                                ! assume only the quadrilateral face on boundary
         if ( ipord .gt. 1 ) then
            nshapeb = nshapeb + 4*(ipord - 1)
         endif
         if ( ipord .gt. 3 ) then
            nshape = nshape +(ipord - 2)*(ipord - 3)/2
         endif

      endif

      select case (intg(1,1))
      case (1)
         nint(5) = 1
      case (2)
         nint(5) = 8
      case (3)
         nint(5) = 27
      case (4)
         nint(5) = 64
      case (5)
         nint(5) = 125
      end select
      select case (intg(2,1))
      case (1)
         nintb(5) = 1
      case (2)
         nintb(5) = 4
      case (3)
         nintb(5) = 9
      case (4)
         nintb(5) = 16
      case (5)
         nintb(5) = 25
      end select
      select case (intg(2,1))
      case (1)
         nintb(6) = 1
      case (2)
         nintb(6) = 4
      case (3)
         nintb(6) = 9
      case (4)
         nintb(6) = 12
      case (5)
         nintb(6) = 12 ! tri face boundary integration rules not coded
      end select

      allocate (tmpQpt (4,nint(5)))
      allocate (tmpQwt (nint(5)))

      call sympyr(nint(5),tmpQpt,tmpQwt,nerr) ! interior elements
      Qpt(5,1:4,1:nint(5)) = tmpQpt(1:4,1:nint(5))
      Qwt(5,1:nint(5)) = tmpQwt(1:nint(5))

      deallocate (tmpQpt)
      deallocate (tmpQwt)

      allocate (tmpQptb(4,nintb(5)))
      allocate (tmpQwtb(nintb(5)))

      call symquad (nintb(5),tmpQptb,tmpQwtb,nerr) ! quad boundary elements
      Qptb(5,1:4,1:nintb(5)) = tmpQptb(1:4,1:nintb(5))
      Qwtb(5,1:nintb(5)) = tmpQwtb(1:nintb(5))

      deallocate (tmpQptb)
      deallocate (tmpQwtb)

      allocate (tmpQptb(4,nintb(6)))
      allocate (tmpQwtb(nintb(6)))

      call symtripyr (nintb(6),tmpQptb,tmpQwtb,nerr) ! tri boundary elements
      Qptb(6,1:4,1:nintb(6)) = tmpQptb(1:4,1:nintb(6))
      Qwtb(6,1:nintb(6)) = tmpQwtb(1:nintb(6))

      deallocate (tmpQptb)
      deallocate (tmpQwtb)

      if (nen.eq.6) then
c
c     later for hierarchic basis nshape formulas will be inserted
c
         nshape = nen

         if ( ipord .gt. 1 ) then
            nshape = nshape + 9*(ipord - 1)
         endif

         if ( ipord .gt. 2 ) then
            nshape = nshape + (ipord -1)*(ipord - 2)
         endif

         if ( ipord .gt. 3 ) then
            nshape = nshape + 3*(ipord - 2)*(ipord - 3)/2
         endif

         if ( ipord .gt. 4 ) then
            nshape = nshape + (ipord - 2)*(ipord - 3)*(ipord - 4)/6
         endif

         nshapeb = nenb
c
         if (nenb .eq. 3) then
            if ( ipord .gt. 1 ) then
               nshapeb = nshapeb + 3*(ipord - 1)
            endif
            if ( ipord .gt. 2 ) then
               nshape = nshape +(ipord - 1)*(ipord - 2)/2
            endif
         endif
c
         if (nenb .eq. 4) then
            if ( ipord .gt. 1 ) then
               nshapeb = nshapeb + 4*(ipord - 1)
            endif
            if ( ipord .gt. 2 ) then
               nshape = nshape +(ipord - 2)*(ipord - 3)/2
            endif
         endif
      endif

      select case (intg(1,1))
      case (2)
         nint(3) = 6
      case (3)
         nint(3) = 18
      case (4)
         nint(3) = 48
      end select
      select case (intg(2,1))
      case (2)
         nintb(3) = 3
         nintb(4) = 4
      case (3)
         nintb(3) = 6
         nintb(4) = 9
      case (4)
         nintb(3) = 12
         nintb(4) = 16
      end select

      allocate (tmpQpt (4,nint(3)))
      allocate (tmpQwt (nint(3)))

      call symwdg(nint(3),tmpQpt,tmpQwt,nerr) ! interior elements
      Qpt(3,1:4,1:nint(3)) = tmpQpt(1:4,1:nint(3))
      Qwt(3,1:nint(3)) = tmpQwt(1:nint(3))

      deallocate (tmpQpt)
      deallocate (tmpQwt)

      allocate (tmpQptb(4,nintb(3)))
      allocate (tmpQwtb(nintb(3)))

      call symtri(nintb(3),tmpQptb,tmpQwtb,nerr) ! boundary elements
      Qptb(3,1:2,2:nintb(3)) = tmpQptb(1:2,1:nintb(3)-1)
      Qptb(3,1:2,1) = tmpQptb(1:2,nintb(3))
c
c     wedges want the third entry to be zeta=-1 (not t=1-r-s)
c     4th entry not used
c
      Qptb(3,3:4,1:nintb(3)) = -1
c$$$  Qptb(3,1:4,1:nintb(3)) = tmpQptb(1:4,1:nintb(3))
      Qwtb(3,1:nintb(3)) = tmpQwtb(1:nintb(3))

      deallocate (tmpQptb)
      deallocate (tmpQwtb)

      allocate (tmpQptb(4,nintb(4)))
      allocate (tmpQwtb(nintb(4)))

      call symquadw(nintb(4),tmpQptb,tmpQwtb,nerr) ! boundary elements
      Qptb(4,1:4,1:nintb(4)) = tmpQptb(1:4,1:nintb(4))
      Qwtb(4,1:nintb(4)) = tmpQwtb(1:nintb(4))

      deallocate (tmpQptb)
      deallocate (tmpQwtb)

c
c.... return
c
      return
      end