genxyz.F90

Go to the documentation of this file.

!-------------------------------------------------------------------
!-------------------------------------------------------------------
subroutine setdef()
  use kinds, only : dp
  use size_m, only : nid, nelt
  use input, only : param, ifmvbd, if3d, ccurve, xc, yc, zc
  use mesh, only : ifdfrm
  implicit none

  real(DP) :: XCC(8),YCC(8),ZCC(8)
  integer :: INDX(8)
  REAL(DP) :: VEC(3,12)
  LOGICAL :: IFVCHK

  integer :: ie, iedg, i, i1, j, i2
  real(DP) :: veclen

! Corner notation:

!                4+-----+3    ^ Y
!                /     /|     |
!               /     / |     |
!             8+-----+7 +2    +----> X
!              |     | /     /
!              |     |/     /
!             5+-----+6    Z


  DO ie=1,nelt
      ifdfrm(ie)= .false. 
  END DO

  IF (ifmvbd) return

!   Force IFDFRM=.true. for all elements (for timing purposes only)

  IF (param(59) /= 0 .AND. nid == 0) &
  write(6,*) 'NOTE: All elements deformed , param(59) ^=0'
  IF (param(59) /= 0) return

!   Check against cases which won't allow for savings in HMHOLTZ

  indx(1)=1
  indx(2)=2
  indx(3)=4
  indx(4)=3
  indx(5)=5
  indx(6)=6
  indx(7)=8
  indx(8)=7

!   Check for deformation (rotation is acceptable).

  DO 500 ie=1,nelt
  
      vec = 0._dp
      IF (if3d) THEN
          DO 100 iedg=1,8
              IF(ccurve(iedg,ie) /= ' ') THEN
                  ifdfrm(ie)= .true. 
                  goto 500
              ENDIF
          100 END DO
      
          DO 105 i=1,8
              xcc(i)=xc(indx(i),ie)
              ycc(i)=yc(indx(i),ie)
              zcc(i)=zc(indx(i),ie)
          105 END DO
      
          DO 110 i=1,4
              vec(1,i)=xcc(2*i)-xcc(2*i-1)
              vec(2,i)=ycc(2*i)-ycc(2*i-1)
              vec(3,i)=zcc(2*i)-zcc(2*i-1)
          110 END DO
      
          i1=4
          DO i=0,1
              DO j=0,1
                  i1=i1+1
                  i2=4*i+j+3
                  vec(1,i1)=xcc(i2)-xcc(i2-2)
                  vec(2,i1)=ycc(i2)-ycc(i2-2)
                  vec(3,i1)=zcc(i2)-zcc(i2-2)
              enddo
          enddo
      
          i1=8
          DO 130 i=5,8
              i1=i1+1
              vec(1,i1)=xcc(i)-xcc(i-4)
              vec(2,i1)=ycc(i)-ycc(i-4)
              vec(3,i1)=zcc(i)-zcc(i-4)
          130 END DO
      
          DO 140 i=1,12
              veclen = vec(1,i)**2 + vec(2,i)**2 + vec(3,i)**2
              veclen = sqrt(veclen)
              vec(1,i)=vec(1,i)/veclen
              vec(2,i)=vec(2,i)/veclen
              vec(3,i)=vec(3,i)/veclen
          140 END DO
      
      !        Check the dot product of the adjacent edges to see that it is zero.
      
          ifdfrm(ie)= .false. 
          IF (  ifvchk(vec,1,5, 9)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,1,6,10)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,2,5,11)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,2,6,12)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,3,7, 9)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,3,8,10)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,4,7,11)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,4,8,12)  ) ifdfrm(ie)= .true. 
      
      !      Check the 2D case....
      
      ELSE
      
          DO 200 iedg=1,4
              IF(ccurve(iedg,ie) /= ' ') THEN
                  ifdfrm(ie)= .true. 
                  goto 500
              ENDIF
          200 END DO
      
          DO 205 i=1,4
              xcc(i)=xc(indx(i),ie)
              ycc(i)=yc(indx(i),ie)
          205 END DO
      
          vec(1,1)=xcc(2)-xcc(1)
          vec(1,2)=xcc(4)-xcc(3)
          vec(1,3)=xcc(3)-xcc(1)
          vec(1,4)=xcc(4)-xcc(2)
          vec(1,5)=0.0
          vec(2,1)=ycc(2)-ycc(1)
          vec(2,2)=ycc(4)-ycc(3)
          vec(2,3)=ycc(3)-ycc(1)
          vec(2,4)=ycc(4)-ycc(2)
          vec(2,5)=0.0
      
          DO 220 i=1,4
              veclen = vec(1,i)**2 + vec(2,i)**2
              veclen = sqrt(veclen)
              vec(1,i)=vec(1,i)/veclen
              vec(2,i)=vec(2,i)/veclen
          220 END DO
      
      !        Check the dot product of the adjacent edges to see that it is zero.
      
          ifdfrm(ie)= .false. 
          IF (  ifvchk(vec,1,3,5)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,1,4,5)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,2,3,5)  ) ifdfrm(ie)= .true. 
          IF (  ifvchk(vec,2,4,5)  ) ifdfrm(ie)= .true. 
      ENDIF
  500 END DO
  return
end subroutine setdef

LOGICAL FUNCTION ifvchk(VEC,I1,I2,I3)
  use kinds, only : dp
  implicit none

  real(DP) :: vec(3,12)
  integer :: i1, i2, i3
  LOGICAL :: IFTMP

  real(DP) :: epsm, dot1, dot2, dot3, dot

  iftmp= .false. 
  epsm=1.0e-06

  dot1=vec(1,i1)*vec(1,i2)+vec(2,i1)*vec(2,i2)+vec(3,i1)*vec(3,i2)
  dot2=vec(1,i2)*vec(1,i3)+vec(2,i2)*vec(2,i3)+vec(3,i2)*vec(3,i3)
  dot3=vec(1,i1)*vec(1,i3)+vec(2,i1)*vec(2,i3)+vec(3,i1)*vec(3,i3)

  dot1=abs(dot1)
  dot2=abs(dot2)
  dot3=abs(dot3)
  dot=dot1+dot2+dot3
  IF (dot > epsm) iftmp= .true. 

  ifvchk=iftmp
  return
END FUNCTION ifvchk

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
subroutine gencoor ()!xm3,ym3,zm3)
  use size_m, only : nx1, ny1, nz1
  use geom, only : ifgmsh3, xm1, ym1, zm1
  implicit none

!max  real(DP) :: XM3(LX3,LY3,LZ3,1),YM3(LX3,LY3,LZ3,1),ZM3(LX3,LY3,LZ3,1)

!   Select appropriate mesh
  IF ( ifgmsh3 ) THEN
    write(*,*) "Oops: IFGMSH3"
!max        CALL GENXYZ (XM3,YM3,ZM3,NX3,NY3,NZ3)
  ELSE
      CALL genxyz(xm1,ym1,zm1,nx1,ny1,nz1)
  ENDIF

  return
end subroutine gencoor

!-----------------------------------------------------------------------
subroutine genxyz (xml,yml,zml,nxl,nyl,nzl)
  use kinds, only : dp
  use size_m, only : lx1, nelt, ndim
  use input, only : ccurve, ifaxis
  implicit none

  integer :: nxl, nyl, nzl
  real(DP) :: xml(nxl,nyl,nzl,1),yml(nxl,nyl,nzl,1),zml(nxl,nyl,nzl,1)

  real(DP) :: h(lx1,3,2), zgml(lx1,3)

  character(1) :: ccv
  integer :: ie, nfaces, iface, isid


#ifdef MOAB
! already read/initialized vertex positions
  if (ifmoab) return
#endif

!   Initialize geometry arrays with bi- triquadratic deformations
  call linquad(xml,yml,zml,nxl,nyl,nzl)

  do ie=1,nelt

      call setzgml(zgml,ie,nxl,nyl,nzl,ifaxis)
      call sethmat(h,zgml,nxl,nyl,nzl)

  !        Deform surfaces - general 3D deformations
  !                        - extruded geometry deformations
      nfaces = 2*ndim
      do iface=1,nfaces
          ccv = ccurve(iface,ie)
          if (ccv == 's') then
            write(*,*) "Oops: ccv = 's'"
!max            call sphsrf(xml,yml,zml,iface,ie,nxl,nyl,nzl,work)
          endif
          if (ccv == 'e') then
            write(*,*) "Oops: ccv = 'e'" 
!max            call gensrf(xml,yml,zml,iface,ie,nxl,nyl,nzl,zgml)
          endif
      enddo

      do isid=1,8
          ccv = ccurve(isid,ie)
          if (ccv == 'C') then
            write(*,*) "Oops: ccv = 'C'"
!max          call arcsrf(xml,yml,zml,nxl,nyl,nzl,ie,isid)
          endif
      enddo

  enddo

  return
end subroutine genxyz

!-----------------------------------------------------------------------
subroutine sethmat(h,zgml,nxl,nyl,nzl)
  use kinds, only : dp
  use size_m, only : lx1
  use input, only : if3d
  implicit none

  real(DP) :: h(lx1,3,2),zgml(lx1,3)
  integer :: nxl, nyl, nzl

  integer :: ix, iy, iz

  do 10 ix=1,nxl
      h(ix,1,1)=(1.0-zgml(ix,1))*0.5
      h(ix,1,2)=(1.0+zgml(ix,1))*0.5
  10 END DO
  do 20 iy=1,nyl
      h(iy,2,1)=(1.0-zgml(iy,2))*0.5
      h(iy,2,2)=(1.0+zgml(iy,2))*0.5
  20 END DO
  if (if3d) then
      do 30 iz=1,nzl
          h(iz,3,1)=(1.0-zgml(iz,3))*0.5
          h(iz,3,2)=(1.0+zgml(iz,3))*0.5
      30 END DO
  else
      h(:,3,:) = 1._dp
  endif

  return
end subroutine sethmat

!-----------------------------------------------------------------------
subroutine setzgml (zgml,e,nxl,nyl,nzl,ifaxl)
  use kinds, only : dp
  use size_m, only : lx1, nx1, nx3, ny1, nz1
  use geom, only : ifgmsh3, ifrzer
  use wz_m, only : zgm1
  implicit none

  real(DP), intent(out) :: zgml(lx1,3)
  integer :: e, nxl, nyl, nzl
  logical :: ifaxl

  integer :: k
  real(DP) :: zam1

  zgml = 0._dp

  if (nxl == 3 .AND. .NOT. ifaxl) then
      do k=1,3
          zgml(1,k) = -1
          zgml(2,k) =  0
          zgml(3,k) =  1
      enddo
  elseif (ifgmsh3 .AND. nxl == nx3) then
    write(*,*) "Oops: IFGMSH3"
#if 0
      call copy(zgml(1,1),zgm3(1,1),nx3)
      call copy(zgml(1,2),zgm3(1,2),ny3)
      call copy(zgml(1,3),zgm3(1,3),nz3)
      if (ifaxl .AND. ifrzer(e)) call copy(zgml(1,2),zam3,ny3)
#endif
  elseif (nxl == nx1 .and. nyl == ny1 .and. nzl == nz1) then
      call copy(zgml(1,1),zgm1(1,1),nxl)
      call copy(zgml(1,2),zgm1(1,2),nyl)
      call copy(zgml(1,3),zgm1(1,3),nzl)
      if (ifaxl .AND. ifrzer(e)) call copy(zgml(1,2),zam1,ny1)
  else
      call exitti('ABORT setzgml! $',nxl)
  endif

  return
  end subroutine setzgml

!-----------------------------------------------------------------------
REAL(DP) FUNCTION dot(V1,V2,N)
  use kinds, only : dp
  implicit none

  integer :: n
  real(DP) :: V1(n),V2(n)
 
  real(DP) :: sum
  integer :: i

  sum = 0
  DO 100 i=1,n
      sum = sum + v1(i)*v2(i)
  100 END DO
  dot = sum
  return
END FUNCTION dot

!-----------------------------------------------------------------------
subroutine linquad(xl,yl,zl,nxl,nyl,nzl)
  use kinds, only : dp
  use size_m, only : ndim, nelt, lx1
  use input, only : ccurve, ifaxis
  implicit none

  integer :: nxl, nyl, nzl
  real(DP) :: xl(nxl*nyl*nzl,*),yl(nxl*nyl*nzl,*),zl(nxl*nyl*nzl,*)

  integer :: e, nedge, k
  logical :: ifmid

  nedge = 4 + 8*(ndim-2)

  do e=1,nelt ! Loop over all elements

      ifmid = .false. 
      do k=1,nedge
          if (ccurve(k,e) == 'm') ifmid = .true. 
      enddo

      if (lx1 == 2) ifmid = .false. 
      if (ifmid) then
        write(*,*) "Oops: ifmid"
!max            call xyzquad(xl(1,e),yl(1,e),zl(1,e),nxl,nyl,nzl,e)
      else
          call xyzlin(xl(1,e),yl(1,e),zl(1,e),nxl,nyl,nzl,e,ifaxis)
      endif
  enddo

  return
end subroutine linquad

!-----------------------------------------------------------------------
subroutine xyzlin(xl,yl,zl,nxl,nyl,nzl,e,ifaxl)
  use kinds, only : dp
  use size_m, only : lx1, ly1, lz1, ndim
  use input, only : xc, yc, zc
  implicit none

  integer :: nxl, nyl, nzl, e
  real(DP) :: xl(nxl,nyl,nzl),yl(nxl,nyl,nzl),zl(nxl,nyl,nzl)
  logical :: ifaxl ! local ifaxis specification

! Preprocessor Corner notation:      Symmetric Corner notation:

!         4+-----+3    ^ s                    3+-----+4    ^ s
!         /     /|     |                      /     /|     |
!        /     / |     |                     /     / |     |
!      8+-----+7 +2    +----> r            7+-----+8 +2    +----> r
!       |     | /     /                     |     | /     /
!       |     |/     /                      |     |/     /
!      5+-----+6    t                      5+-----+6    t

  integer, save :: indx(8) = (/ 1,2,4,3,5,6,8,7 /)

  integer, parameter :: ldw=4*lx1*ly1*lz1
  real(DP) :: xcb(2,2,2),ycb(2,2,2),zcb(2,2,2),w(ldw)

!  real(DP) :: zgml, jx,jy,jz,jxt,jyt,jzt, zlin
  real(DP) :: zgml(lx1,3),jx (lx1*2)
  real(DP) :: jxt(lx1*2),jyt(lx1*2),jzt(lx1*2),zlin(2)

  integer :: i, k, ndim2

  call setzgml(zgml,e,nxl,nyl,nzl,ifaxl)

  zlin(1) = -1
  zlin(2) =  1

  k = 1
  do i=1,nxl
      call fd_weights_full(zgml(i,1),zlin,1,0,jxt(k))
      call fd_weights_full(zgml(i,2),zlin,1,0,jyt(k))
      call fd_weights_full(zgml(i,3),zlin,1,0,jzt(k))
      k=k+2
  enddo
  call transpose(jx,nxl,jxt,2)

  ndim2 = 2**ndim
  ! Convert prex notation to lexicographical
  xcb = reshape(xc(indx(1:ndim2),e), (/2,2,2/))
  ycb = reshape(yc(indx(1:ndim2),e), (/2,2,2/))
  zcb = reshape(zc(indx(1:ndim2),e), (/2,2,2/))

!   Map R-S-T space into physical X-Y-Z space.

! NOTE:  Assumes nxl=nyl=nzl !

  call tensr3(xl,nxl,xcb,2,jx,jyt,jzt,w)
  call tensr3(yl,nxl,ycb,2,jx,jyt,jzt,w)
  call tensr3(zl,nxl,zcb,2,jx,jyt,jzt,w)

  return
end subroutine xyzlin
!-----------------------------------------------------------------------