connect1.F90

Go to the documentation of this file.

!-----------------------------------------------------------------------
subroutine setup_topo()
  use kinds, only : dp, i8
  use size_m, only : nid, ndim, nx1, ny1, nz1, nelv, nelt, nfield
  use size_m, only : lx1, ly1, lz1, lelv
  use input, only : ifflow, ifmvbd, ifheat, param
  use mesh, only : vertex
  use mvgeom, only : wmult
  use parallel, only : nelgv, nelgt, gsh_fld, nelg
  use soln, only : vmult, tmult
  use tstep, only : ifield
  use zper, only : ifgtp
  implicit none

  integer(i8), allocatable :: glo_num(:)

  integer :: nxl, nyl, nzl, mfield, ncrnr, ntotv, ntott
  real(DP) :: vmltmax, ivmltmax
  real(DP), external :: glmax

  allocate(glo_num(lx1*ly1*lz1*lelv))

  if(nid == 0) write(6,*) 'setup mesh topology'

!   Initialize key arrays for Direct Stiffness SUM.

  nxl=3
  nyl=3
  nzl=1+2*(ndim-2)

  call initds
  call dsset(nx1,ny1,nz1)
  call setedge

!=================================================
!     Establish (global) domain topology
!=================================================
!     .Generate topologically correct mesh data.
!     .Set up element centers, face centers, etc.
!     .Check  right handedness of elements.
!     .Check  element boundary conditions.
!     .Establish Element-Element rotations
!     .Construct the element to processor map and

  call genxyzl
  call setside
  if (param(75) < 1) call verify

  CALL setcdof
!max  IF (IFAXIS            ) CALL SETRZER
!max    IF (IFMVBD            ) CALL CBCMESH
#if 0
  IF (ifmodel .AND. ifkeps) CALL cbcturb
#endif
  if (param(75) < 1) CALL chkaxcb

!========================================================================
!     Set up element-processor mapping and establish global numbering
!========================================================================

  mfield=2
  if (ifflow) mfield=1
  if (ifmvbd) mfield=0

  ncrnr = 2**ndim

  if (nelgv == nelgt) then
      if (ifgtp) then
        write(*,*) "Oops: ifgtp"
!max          call gen_gtp_vertex    (vertex, ncrnr)
      else
          call get_vert
      endif
      call setupds(gsh_fld(1),nx1,ny1,nz1,nelv,nelgv,vertex,glo_num)
      gsh_fld(2)=gsh_fld(1)

  !        call gs_counter(glo_num,gsh_fld(1))

  else

  !        For conjugate heat transfer, it is assumed that fluid
  !        elements are listed both globally and locally with lower
  !        element numbers than the solid elements.
  !        We currently assume that there is at least one fluid elem.
  !        per processor.
  

      call get_vert
  !        call outmati(vertex,4,nelv,'vrtx V')
      call setupds(gsh_fld(1),nx1,ny1,nz1,nelv,nelgv,vertex,glo_num)

  !        call get_vert  (vertex, ncrnr, nelgt, '.mp2')  !  LATER !
  !        call outmati(vertex,4,nelt,'vrtx T')
      call setupds(gsh_fld(2),nx1,ny1,nz1,nelt,nelgt,vertex,glo_num)

  
  !        Feb 20, 2012:  It appears that we do not need this restriction: (pff)
  !        check if there is a least one fluid element on each processor
  !        do iel = 1,nelt
  !           ieg = lglel(iel)
  !           if (ieg.le.nelgv) goto 101
  !        enddo
  !        if(nid.eq.0) write(6,*)
  !    &     'ERROR: each domain must contain at least one fluid element!'
  !        call exitt
  ! 101   continue

  endif

!max     if (ifmvbd) call setup_mesh_dssum ! Set up dssum for mesh

!========================================================================
!     Set up multiplicity and direct stiffness arrays for each IFIELD
!========================================================================

  ntotv = nx1*ny1*nz1*nelv
  ntott = nx1*ny1*nz1*nelt

  if (ifheat) then
      ifield = 2
      allocate(tmult(nx1,ny1,nz1,nelt,1))
      tmult(:,:,:,:,1) = 1._dp
      call dssum(tmult)
      tmult(:,:,:,:,1) = 1._dp / tmult(:,:,:,:,1)
  endif

  if (ifflow) then
    if (ifheat) then
      vmult => tmult(:,:,:,:,1)
    else
      ifield = 1
      allocate(vmult(nx1,ny1,nz1,nelv))
      vmult = 1._dp
      call dssum(vmult)
      vmltmax=glmax(vmult,ntotv)
      ivmltmax=vmltmax
      if (nid == 0) write(6,*) ivmltmax,' max multiplicity'
      vmult = 1._dp / vmult
    endif
  endif
  if ( .NOT. ifflow) call copy(vmult,tmult,ntott)
  if (ifmvbd)  call copy(wmult,vmult,ntott)
  do ifield=3,nfield                  ! Additional pass. scalrs.
      if (nelg(ifield) == nelgv) then
          gsh_fld(ifield) = gsh_fld(1)
          call copy(tmult(1,1,1,1,ifield-1),vmult,ntotv)
      else
          gsh_fld(ifield) = gsh_fld(2)
          call copy(tmult(1,1,1,1,ifield-1),tmult,ntott)
      endif
  enddo

  if(nid == 0) then
      write(6,*) 'done :: setup mesh topology'
      write(6,*) ' '
  endif

  return
end subroutine setup_topo

!-----------------------------------------------------------------------
subroutine initds
  use size_m
  use topol, only : group, eface1, eface, nomlis
  implicit none

  integer :: j, idim, iface

!   Nominal ordering for direct stiffness summation of faces
  j=0
  DO idim=1,ndim
      DO iface=1,2
          j=j+1
          nomlis(iface,idim)=j
      enddo
  END DO

!   Assign Ed's numbering scheme to PF's scheme.

  eface(1)=4
  eface(2)=2
  eface(3)=1
  eface(4)=3
  eface(5)=5
  eface(6)=6

!   Assign inverse of Ed's numbering scheme to PF's scheme.

  eface1(1)=3
  eface1(2)=2
  eface1(3)=4
  eface1(4)=1
  eface1(5)=5
  eface1(6)=6

!   Assign group designation to each face to determine ordering of indices.

  group(1)=0
  group(2)=1
  group(3)=1
  group(4)=0
  group(5)=0
  group(6)=1

  RETURN
end subroutine initds

!-----------------------------------------------------------------------
subroutine setedge()
  use size_m, only : ndim
  use topol, only : iedgef, iedge, invedg, skpdat, group
  implicit none

  integer :: ITMP(3,3,3)
  INTEGER :: ORDER

  integer :: nxl, nyl, nzl, nxy, nxyz, nfaces
  integer :: i3d, i, i3, iz, i2, iy, i1, ix
  integer :: j, j1, j2, jskip2
  integer :: jf2, js2, jskip1, jf1, js1, iface, image, ii
  integer :: iedgef_flat(48)
  nxl=3
  nyl=3
  nzl=1+2*(ndim-2)
  nxy   =nxl*nyl
  nxyz  =nxl*nyl*nzl
  nfaces=2*ndim

!----------------------------------------------------------------------
!     Set up edge arrays (temporary - required only for defining DS)
!----------------------------------------------------------------------
!     Fill corners - 1 through 8.

  i3d=1
  IF (ndim == 2) i3d=0

  i=0
  DO i3=0,i3d
      iz=1+(nzl-1)*i3
      DO i2=0,1
          iy=1+(nyl-1)*i2
          DO i1=0,1
              ix=1+(nxl-1)*i1
              i=i+1
              iedge(i)=ix+nxl*(iy-1)+nxy*(iz-1)
          enddo
      enddo
  END DO

!   Fill X-direction edges.

  DO i3=0,i3d
      iz=1+(nzl-1)*i3
      DO i2=0,1
          iy=1+(nyl-1)*i2
          DO ix=2,nxl-1
              i=i+1
              iedge(i)=ix+nxl*(iy-1)+nxy*(iz-1)
          enddo
      enddo
  enddo

!   Fill Y-direction edges.

  DO i3=0,i3d
      iz=1+(nzl-1)*i3
      DO i1=0,1
          ix=1+(nxl-1)*i1
          DO iy=2,nyl-1
              i=i+1
              iedge(i)=ix+nxl*(iy-1)+nxy*(iz-1)
          enddo
      enddo
  enddo

!   Fill Z-direction edges.

  IF (ndim == 3) THEN
      DO i2=0,1
          iy=1+(nyl-1)*i2
          DO i1=0,1
              ix=1+(nxl-1)*i1
              DO iz=2,nzl-1
                  i=i+1
                  iedge(i)=ix+nxl*(iy-1)+nxy*(iz-1)
              enddo
          enddo
      enddo
  ENDIF

  invedg = 0
  DO ii=1,20
      ix=iedge(ii)
      invedg(ix)=ii
  END DO


!   GENERAL FACE, GENERAL ROTATION EDGE NUMBERS.
#if 0

  IF (ndim == 3) THEN
  
  !        Pack 3-D edge numbering:
  
  !        Fill temporary array with local index numbers:
      itmp = reshape( (/(i,i=1,nxyz)/), (/nxl, nyl, nzl/)) 
  
  !        Two sets are required, the base cube and the image cube
  !        which is being summed with it.
  
      DO 1000 image=0,1
      
      !        Pack edges for each face, no rotation.
      
          DO 500 iface=1,nfaces
              js1    = skpdat(1,iface)
              jf1    = skpdat(2,iface)
              jskip1 = skpdat(3,iface)
              js2    = skpdat(4,iface)
              jf2    = skpdat(5,iface)
              jskip2 = skpdat(6,iface)
          
          !           Choose proper indexing order according to face type and image.
          
              order = (-1)**(group(iface)+image)
              IF (order == 1) THEN
              
              !              Forward ordering:
              
              !            +-------------+    ^ v1
              !            | --------->| |    |
              !            | ^    2    | |    +-->
              !            | |         | |      v2
              !            | |1       3| |
              !            | |    4    V |
              !            | |<--------- |
              !            F-------------I     F is fiducial node.
              
              !                                I is location of fiducial node for
              !                                     image face.
              
              !           Load edge 1:
              
                  j=0
                  j2=js2
                  DO 100 j1=js1,jf1-jskip1,jskip1
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  100 END DO
              
              !           Load edge 2:
              
                  j=0
                  j1=jf1
                  DO 200 j2=js2,jf2-jskip2,jskip2
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  200 END DO
              
              
              !           Load edge 3:
              
                  j=0
                  j2=jf2
                  DO 300 j1=jf1,js1+jskip1,-jskip1
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  300 END DO
              
              !           Load edge 4:
              
                  j=0
                  j1=js1
                  DO 400 j2=jf2,js2+jskip2,-jskip2
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  400 END DO
              ! insert
              iedgef(:,:,:,image) = reshape(iedgef_flat, (/2,4,6/)) 

              ELSE
              
              !           Reverse ordering:
              
              !            +-------------+
              !            | |<--------- |       ^ v2
              !            | |    2    ^ |       |
              !            | |         | |    <--+
              !            | |3       1| |     v1
              !            | V    4    | |
              !            | --------->| |
              !            I-------------F     F is fiducial node.
              
              !                                I is location of fiducial node for
              !                                     image face.
              
              !           Load edge 1:
              
                  j=0
                  j1=js1
                  write(*,*) iface, js2, jf2, jskip2, image
                  DO 105 j2=js2,jf2-jskip2,jskip2
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  105 END DO
              
              !           Load edge 2:
              
                  j=0
                  j2=jf2
                  DO 205 j1=js1,jf1-jskip1,jskip1
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  205 END DO
              
              !           Load edge 3:
              
                  j=0
                  j1=jf1
                  DO 305 j2=jf2,js2+jskip2,-jskip2
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  305 END DO
              
              !           Load edge 4:
              
                  j=0
                  j2=js2
                  DO 405 j1=jf1,js1+jskip1,-jskip1
                      j=j+1
                      iedgef_flat(j+(iface-1)*8)=j1 + 3*(j2-1) !ITMP(J1,J2,1)
                  405 END DO
              ENDIF
              ! insert
              iedgef(:,:,:,image) = reshape(iedgef_flat, (/2,4,6/)) 
          
          500 END DO
      1000 END DO
  ELSE
  
  !        Load edge information for 2-D case
  
      iedgef(1,1,1,0) = nxy - nxl + 1
      iedgef(1,2,1,0) = 1
      iedgef(1,1,2,0) = nxl
      iedgef(1,2,2,0) = nxy
      iedgef(1,1,3,0) = 1
      iedgef(1,2,3,0) = nxl
      iedgef(1,1,4,0) = nxy
      iedgef(1,2,4,0) = nxy - nxl + 1
  
      iedgef(1,1,1,1) = 1
      iedgef(1,2,1,1) = nxy - nxl + 1
      iedgef(1,1,2,1) = nxy
      iedgef(1,2,2,1) = nxl
      iedgef(1,1,3,1) = nxl
      iedgef(1,2,3,1) = 1
      iedgef(1,1,4,1) = nxy - nxl + 1
      iedgef(1,2,4,1) = nxy
  ENDIF
#endif

  RETURN
end subroutine setedge

!-----------------------------------------------------------------------
subroutine dsset(nx,ny,nz)
  use topol, only : ixcn, skpdat, eskip, nedg
  implicit none

  integer :: nx, ny, nz

  INTEGER, save :: NXO = 0, NYO = 0, NZO = 0
  integer :: ic, icz, icy, icx, nxy, ied, iedm

!   Check if element surface counters are already set from last call...
  IF (nxo == nx .AND. nyo == ny .AND. nzo == nz) RETURN

!   else, proceed....

  nxo = nx
  nyo = ny
  nzo = nz

!   Establish corner to elemental node number mappings

  ic=0
  DO icz=0,1
      DO icy=0,1
          DO icx=0,1
          !       Supress vectorization to
          !        IF(ICX.EQ.0)DUMMY=0
          !        IF(ICX.EQ.1)DUMMY=1
          !        DUMMY2=DUMMY2+DUMMY
              ic=ic+1
              ixcn(ic)= 1 + (nx-1)*icx + nx*(ny-1)*icy + nx*ny*(nz-1)*icz
          enddo
      enddo
  enddo

!   Assign indices for direct stiffness summation of arbitrary faces.


!   Y-Z Planes (Faces 1 and 2)

  skpdat(1,1)=1
  skpdat(2,1)=nx*(ny-1)+1
  skpdat(3,1)=nx
  skpdat(4,1)=1
  skpdat(5,1)=ny*(nz-1)+1
  skpdat(6,1)=ny

  skpdat(1,2)=1             + (nx-1)
  skpdat(2,2)=nx*(ny-1)+1   + (nx-1)
  skpdat(3,2)=nx
  skpdat(4,2)=1
  skpdat(5,2)=ny*(nz-1)+1
  skpdat(6,2)=ny

!   X-Z Planes (Faces 3 and 4)

  skpdat(1,3)=1
  skpdat(2,3)=nx
  skpdat(3,3)=1
  skpdat(4,3)=1
  skpdat(5,3)=ny*(nz-1)+1
  skpdat(6,3)=ny

  skpdat(1,4)=1           + nx*(ny-1)
  skpdat(2,4)=nx          + nx*(ny-1)
  skpdat(3,4)=1
  skpdat(4,4)=1
  skpdat(5,4)=ny*(nz-1)+1
  skpdat(6,4)=ny

!   X-Y Planes (Faces 5 and 6)

  skpdat(1,5)=1
  skpdat(2,5)=nx
  skpdat(3,5)=1
  skpdat(4,5)=1
  skpdat(5,5)=ny
  skpdat(6,5)=1

  skpdat(1,6)=1           + nx*ny*(nz-1)
  skpdat(2,6)=nx          + nx*ny*(nz-1)
  skpdat(3,6)=1
  skpdat(4,6)=1
  skpdat(5,6)=ny
  skpdat(6,6)=1

!   Set up skip indices for each of the 12 edges

!       Note that NXY = NX*NY even for 2-D since
!       this branch does not apply to the 2D case anyway.

!   ESKIP(*,1) = start location
!   ESKIP(*,2) = end
!   ESKIP(*,3) = stride

  nxy=nx*ny
  eskip( 1,1) = ixcn(1) + 1
  eskip( 1,2) = ixcn(2) - 1
  eskip( 1,3) = 1
  eskip( 2,1) = ixcn(3) + 1
  eskip( 2,2) = ixcn(4) - 1
  eskip( 2,3) = 1
  eskip( 3,1) = ixcn(5) + 1
  eskip( 3,2) = ixcn(6) - 1
  eskip( 3,3) = 1
  eskip( 4,1) = ixcn(7) + 1
  eskip( 4,2) = ixcn(8) - 1
  eskip( 4,3) = 1
  eskip( 5,1) = ixcn(1) + nx
  eskip( 5,2) = ixcn(3) - nx
  eskip( 5,3) = nx
  eskip( 6,1) = ixcn(2) + nx
  eskip( 6,2) = ixcn(4) - nx
  eskip( 6,3) = nx
  eskip( 7,1) = ixcn(5) + nx
  eskip( 7,2) = ixcn(7) - nx
  eskip( 7,3) = nx
  eskip( 8,1) = ixcn(6) + nx
  eskip( 8,2) = ixcn(8) - nx
  eskip( 8,3) = nx
  eskip( 9,1) = ixcn(1) + nxy
  eskip( 9,2) = ixcn(5) - nxy
  eskip( 9,3) = nxy
  eskip(10,1) = ixcn(2) + nxy
  eskip(10,2) = ixcn(6) - nxy
  eskip(10,3) = nxy
  eskip(11,1) = ixcn(3) + nxy
  eskip(11,2) = ixcn(7) - nxy
  eskip(11,3) = nxy
  eskip(12,1) = ixcn(4) + nxy
  eskip(12,2) = ixcn(8) - nxy
  eskip(12,3) = nxy

!   Load reverse direction edge arrays for reverse mappings...

  DO 20 ied=1,12
      iedm=-ied
      eskip(iedm,1) =  eskip(ied,2)
      eskip(iedm,2) =  eskip(ied,1)
      eskip(iedm,3) = -eskip(ied,3)
  20 END DO

!   Compute offset for global edge vector given current element
!   dimensions.....

!   NGSPED(ITE,ICMP) = number of global (ie, distinct) special edges
!                      of type ITE (1,2, or 3)  for field ICMP.

!                      ITE = 1 implies an "X" edge
!                      ITE = 2 implies an "Y" edge
!                      ITE = 3 implies an "Z" edge

!   Set up number of nodes along each of the 3 types of edges
!   (endpoints excluded).

  nedg(1)=nx-2
  nedg(2)=ny-2
  nedg(3)=nz-2

  RETURN
end subroutine dsset

!-----------------------------------------------------------------------
subroutine genxyzl()
  use kinds, only : dp
  use size_m, only : ndim, nelt
  use input, only : xc, yc, zc
  use scratch, only : xml, yml, zml
  implicit none

  real(DP) :: XCB(2,2,2),YCB(2,2,2),ZCB(2,2,2), H(3,3,2)
  integer :: indx(8)

  integer :: nxl, nyl, nzl, ntot3
  integer :: ix, iy, iz, ixt, iyt, izt, ie, ndim2
  nxl=3
  nyl=3
  nzl=1+2*(ndim-2)
  ntot3=nxl*nyl*nzl*nelt

! 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

  DO 10 ix=1,nxl
      h(ix,1,1)=0.5*float(3-ix)
      h(ix,1,2)=0.5*float(ix-1)
  10 END DO
  DO 20 iy=1,nyl
      h(iy,2,1)=0.5*float(3-iy)
      h(iy,2,2)=0.5*float(iy-1)
  20 END DO
  DO 30 iz=1,nzl
      h(iz,3,1)=0.5*float(3-iz)
      h(iz,3,2)=0.5*float(iz-1)
  30 END DO

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

  xml = 0._dp; yml = 0._dp; zml = 0._dp
  xcb = 0._dp; ycb = 0._dp; zcb = 0._dp

  DO 5000 ie=1,nelt
  
      ndim2 = 2**ndim
      xcb = reshape(xc(indx(1:ndim2), ie), (/2,2,2/))
      ycb = reshape(yc(indx(1:ndim2), ie), (/2,2,2/))
      zcb = reshape(zc(indx(1:ndim2), ie), (/2,2,2/))
 
  !        Map R-S-T space into physical X-Y-Z space.
  
      DO izt=1,ndim-1
          DO iyt=1,2
              DO ixt=1,2
              
                  DO iz=1,nzl
                      DO iy=1,nyl
                          DO ix=1,nxl
                              xml(ix,iy,iz,ie)=xml(ix,iy,iz,ie)+ &
                              h(ix,1,ixt)*h(iy,2,iyt)*h(iz,3,izt)*xcb(ixt,iyt,izt)
                              yml(ix,iy,iz,ie)=yml(ix,iy,iz,ie)+ &
                              h(ix,1,ixt)*h(iy,2,iyt)*h(iz,3,izt)*ycb(ixt,iyt,izt)
                              zml(ix,iy,iz,ie)=zml(ix,iy,iz,ie)+ &
                              h(ix,1,ixt)*h(iy,2,iyt)*h(iz,3,izt)*zcb(ixt,iyt,izt)
                          enddo
                      enddo
                  enddo
              enddo
          enddo
      enddo
  
  5000 END DO
  RETURN
end subroutine genxyzl

!-----------------------------------------------------------------------
subroutine verify
  implicit none 
  call verrhe
  return
end subroutine verify

!-----------------------------------------------------------------------
subroutine setside
  use kinds, only : dp
  use size_m, only : nelt, ndim
  use input, only : if3d, xc, yc, zc
  use scratch, only : xyz, side
  use topol, only : indx, icface
  implicit none

  integer :: ie, j, ivtx, nfaces, ncrnr, icrn, ifac, icr1, ivt1, idim
  integer, external :: mod1
  real(DP) :: avwght

!   SIDE(i,IFACE,IE) -  Physical (xyz) location of element side midpoint.
!                       i=1,2,3 gives x,y,z value, respectively.
!                       i=4  gives average dimension of face for setting
!                            tolerances.
  indx(1)=1
  indx(2)=2
  indx(3)=4
  indx(4)=3
  indx(5)=5
  indx(6)=6
  indx(7)=8
  indx(8)=7

!   Flip vertex array structure

!   write(6,*) nelv,nelt,if3d
  xyz = 0._dp
  if (if3d) then
      do ie=1,nelt
          do j=1,8
              ivtx = indx(j)
              xyz(1,ivtx,ie) = xc(j,ie)
              xyz(2,ivtx,ie) = yc(j,ie)
              xyz(3,ivtx,ie) = zc(j,ie)
          !           write(6,1) ie,j,ivtx,xc(j,ie),yc(j,ie),zc(j,ie),' xcz'
          !           write(6,1) ie,j,ivtx,(xyz(k,ivtx,ie),k=1,3),' vtx'
          !   1       format(3i5,1p3e12.4,a4)
          enddo
      enddo
  else
      do ie=1,nelt
          do j=1,4
              ivtx = indx(j)
              xyz(1,ivtx,ie) = xc(j,ie)
              xyz(2,ivtx,ie) = yc(j,ie)
              xyz(3,ivtx,ie) = 0.0
          enddo
      enddo
  endif

!   Compute location of center and "diameter" of each element side.

  nfaces=ndim*2
  ncrnr =2**(ndim-1)
  side = 0._dp
  DO icrn=1,ncrnr
      DO ifac=1,nfaces
          ivtx = icface(icrn,ifac)
          icr1 = ncrnr+(icrn-1)
          icr1 = mod1(icr1,ncrnr)
          ivt1 = icface(icr1,ifac)
          DO 400 ie=1,nelt
              DO 300 idim=1,ndim
                  side(idim,ifac,ie)=side(idim,ifac,ie)+xyz(idim,ivtx,ie)
                  side(   4,ifac,ie)=side(   4,ifac,ie)+ &
                  ( xyz(idim,ivtx,ie)-xyz(idim,ivt1,ie) )**2
              300 END DO
              side(4,ifac,ie)=sqrt( side(4,ifac,ie) )
          400 END DO
      enddo
  END DO
  avwght=1.0/float(ncrnr)
  side = side*avwght

!   call exitt
  RETURN
end subroutine setside

!-----------------------------------------------------------------------
subroutine verrhe()
  use kinds, only : dp
  use size_m, only : nid, nelt
  use input, only : if3d
  use parallel, only : lglel, nelgt
  use scratch, only : xyz
  implicit none

  integer :: ie, ieg
  real(DP) :: v1, v2, v3, v4, v5, v6, v7, v8
  LOGICAL :: IFCSTT
  real(DP), external :: volum0

  ifcstt= .true. 
  IF ( .NOT. if3d) THEN
    write(*,*) "Oops: not if3d"
#if 0
      DO 1000 ie=1,nelt
      
      !        CRSS2D(A,B,O) = (A-O) X (B-O)
      
          c1=crss2d(xyz(1,2,ie),xyz(1,3,ie),xyz(1,1,ie))
          c2=crss2d(xyz(1,4,ie),xyz(1,1,ie),xyz(1,2,ie))
          c3=crss2d(xyz(1,1,ie),xyz(1,4,ie),xyz(1,3,ie))
          c4=crss2d(xyz(1,3,ie),xyz(1,2,ie),xyz(1,4,ie))
      
          IF (c1 <= 0.0 .OR. c2 <= 0.0 .OR. &
          c3 <= 0.0 .OR. c4 <= 0.0 ) THEN
          
              ieg=lglel(ie)
              WRITE(6,800) ieg,c1,c2,c3,c4
              call exitt
              800 FORMAT(/,2x,'WARNINGa: Detected non-right-handed element.', &
              /,2x,'Number',i8,'  C1-4:',4e12.4)
              ifcstt= .false. 
          !           CALL QUERY(IFYES,'Proceed                                 ')
          !           IF (.NOT.IFYES) GOTO 9000
          ENDIF
      1000 END DO
#endif 
  !     Else 3-D:
  ELSE
      DO 2000 ie=1,nelt
      
      !        VOLUM0(A,B,C,O) = (A-O)X(B-O).(C-O)
      
          v1= volum0(xyz(1,2,ie),xyz(1,3,ie),xyz(1,5,ie),xyz(1,1,ie))
          v2= volum0(xyz(1,4,ie),xyz(1,1,ie),xyz(1,6,ie),xyz(1,2,ie))
          v3= volum0(xyz(1,1,ie),xyz(1,4,ie),xyz(1,7,ie),xyz(1,3,ie))
          v4= volum0(xyz(1,3,ie),xyz(1,2,ie),xyz(1,8,ie),xyz(1,4,ie))
          v5=-volum0(xyz(1,6,ie),xyz(1,7,ie),xyz(1,1,ie),xyz(1,5,ie))
          v6=-volum0(xyz(1,8,ie),xyz(1,5,ie),xyz(1,2,ie),xyz(1,6,ie))
          v7=-volum0(xyz(1,5,ie),xyz(1,8,ie),xyz(1,3,ie),xyz(1,7,ie))
          v8=-volum0(xyz(1,7,ie),xyz(1,6,ie),xyz(1,4,ie),xyz(1,8,ie))
      
          IF (v1 <= 0.0 .OR. v2 <= 0.0 .OR. &
          v3 <= 0.0 .OR. v4 <= 0.0 .OR. &
          v5 <= 0.0 .OR. v6 <= 0.0 .OR. &
          v7 <= 0.0 .OR. v8 <= 0.0    ) THEN
          
              ieg=lglel(ie)
              WRITE(6,1800) ieg,v1,v2,v3,v4,v5,v6,v7,v8
              call exitt
              1800 FORMAT(/,2x,'WARNINGb: Detected non-right-handed element.', &
              /,2x,'Number',i8,'  V1-8:',4e12.4 &
              /,2x,'      ',4x,'       ',4e12.4)
              ifcstt= .false. 
          ENDIF
      2000 END DO
  ENDIF

!   Print out results from right-handed check

  IF ( .NOT. ifcstt) WRITE(6,2001)

!   Check consistency accross all processors.

  CALL gllog(ifcstt, .false. )

  IF ( .NOT. ifcstt) THEN
      IF (nid == 0) WRITE(6,2003) nelgt
      call exitt
  ELSE
      IF (nid == 0) WRITE(6,2002) nelgt
  ENDIF

  2001 FORMAT(//,'  Elemental geometry not right-handed, ABORTING' &
  ,' in routine VERRHE.')
  2002 FORMAT('   Right-handed check complete for',i8,' elements. OK.')
  2003 FORMAT('   Right-handed check failed for',i8,' elements.' &
  ,'   Exiting in routine VERRHE.')
  RETURN
end subroutine verrhe

!-----------------------------------------------------------------------
real(DP) FUNCTION volum0(P1,P2,P3,P0)
  use kinds, only : dp
  implicit none
  REAL(DP) :: P1(3),P2(3),P3(3),P0(3)
  real(DP) :: u1, u2, u3, v1, v2, v3, w1, w2, w3, cross1, cross2, cross3
  u1=p1(1)-p0(1)
  u2=p1(2)-p0(2)
  u3=p1(3)-p0(3)

  v1=p2(1)-p0(1)
  v2=p2(2)-p0(2)
  v3=p2(3)-p0(3)

  w1=p3(1)-p0(1)
  w2=p3(2)-p0(2)
  w3=p3(3)-p0(3)

  cross1 = u2*v3-u3*v2
  cross2 = u3*v1-u1*v3
  cross3 = u1*v2-u2*v1

  volum0  = w1*cross1 + w2*cross2 + w3*cross3
           
  RETURN
END FUNCTION volum0

!-----------------------------------------------------------------------
subroutine facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
  implicit none
  integer :: kx1, ky1, kz1, kx2, ky2, kz2, nx, ny, nz, iface

  kx1=1
  ky1=1
  kz1=1
  kx2=nx
  ky2=ny
  kz2=nz
  IF (iface == 1) ky2=1
  IF (iface == 2) kx1=nx
  IF (iface == 3) ky1=ny
  IF (iface == 4) kx2=1
  IF (iface == 5) kz2=1
  IF (iface == 6) kz1=nz
  return
end subroutine facind

!-----------------------------------------------------------------------
subroutine facev(a,ie,iface,val,nx,ny,nz)
  use kinds, only : dp
  use size_m, only : lelt
  implicit none
  integer :: ie, iface, nx, ny, nz
  real(DP) :: a(nx,ny,nz,lelt), val

  integer :: kx1, ky1, kz1, kx2, ky2, kz2
  integer :: ix, iy, iz

  CALL facind(kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
  DO iz=kz1,kz2
      DO iy=ky1,ky2
          DO ix=kx1,kx2
              a(ix,iy,iz,ie)=val
          enddo
      enddo
  enddo
  RETURN
end subroutine facev
!-----------------------------------------------------------------------