utilities/RocfracPrep/read_patran.f90

Go to the documentation of this file.

!*********************************************************************
!* Illinois Open Source License                                      *
!*                                                                   *
!* University of Illinois/NCSA                                       * 
!* Open Source License                                               *
!*                                                                   *
!* Copyright@2008, University of Illinois.  All rights reserved.     *
!*                                                                   *
!*  Developed by:                                                    *
!*                                                                   *
!*     Center for Simulation of Advanced Rockets                     *
!*                                                                   *
!*     University of Illinois                                        *
!*                                                                   *
!*     www.csar.uiuc.edu                                             *
!*                                                                   *
!* Permission is hereby granted, free of charge, to any person       *
!* obtaining a copy of this software and associated documentation    *
!* files (the "Software"), to deal with the Software without         *
!* restriction, including without limitation the rights to use,      *
!* copy, modify, merge, publish, distribute, sublicense, and/or      *
!* sell copies of the Software, and to permit persons to whom the    *
!* Software is furnished to do so, subject to the following          *
!* conditions:                                                       *
!*                                                                   *
!*                                                                   *
!* @ Redistributions of source code must retain the above copyright  * 
!*   notice, this list of conditions and the following disclaimers.  *
!*                                                                   * 
!* @ Redistributions in binary form must reproduce the above         *
!*   copyright notice, this list of conditions and the following     *
!*   disclaimers in the documentation and/or other materials         *
!*   provided with the distribution.                                 *
!*                                                                   *
!* @ Neither the names of the Center for Simulation of Advanced      *
!*   Rockets, the University of Illinois, nor the names of its       *
!*   contributors may be used to endorse or promote products derived * 
!*   from this Software without specific prior written permission.   *
!*                                                                   *
!* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,   *
!* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES   *
!* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND          *
!* NONINFRINGEMENT.  IN NO EVENT SHALL THE CONTRIBUTORS OR           *
!* COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER       * 
!* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,   *
!* ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE    *
!* USE OR OTHER DEALINGS WITH THE SOFTWARE.                          *
!*********************************************************************
!* Please acknowledge The University of Illinois Center for          *
!* Simulation of Advanced Rockets in works and publications          *
!* resulting from this software or its derivatives.                  *
!*********************************************************************
SUBROUTINE read_patran(numvertx2d,dhmin,nprocs)
  
  USE linked_list
  USE meshdata
  USE linked_list2
  
  IMPLICIT NONE
  
  INTEGER :: numvertx2d
  REAL*8 :: dhmin

! Local
  REAL*8 :: xx,yy,zz,size1,size2,size3,size4,size5,size6
  REAL*8 :: dt_courant
  REAL*8, DIMENSION(1:6) :: disflagvalue, tmpval
  INTEGER :: numdisflag
  
  INTEGER :: id
  INTEGER :: i, j
  INTEGER :: itype
  INTEGER :: n1,n2,n3,n4,n5,n6,n7,n8
  INTEGER :: iface
  REAL*8 :: value
  REAL*8 :: press
!-- Tempory holding array of partitioned metis  arrays
  INTEGER, ALLOCATABLE, DIMENSION(:) :: npart
  INTEGER :: edgecut
  INTEGER :: nprocs
  INTEGER :: ip
!
  INTEGER :: iaux,iaux1,iaux2
  INTEGER, ALLOCATABLE, DIMENSION(:) :: imin,nmin,imax,nmax,ninc
  INTEGER, ALLOCATABLE, DIMENSION(:) :: imap
  INTEGER, ALLOCATABLE, DIMENSION(:,:) :: lm4node
  INTEGER :: nnd1, nnd2, nnd3, nnd4, nnd5, nnd6, nnd7, nnd8, nnd9, nnd10
  INTEGER :: numnp4node
  integer, dimension(1:6) :: disflag
  INTEGER :: ifaceid, itmp
  
  
  INTEGER :: procid
  
  
  LOGICAL :: hex8, tet4, tet10, fileexist
  LOGICAL :: fileexist1, fileexist2, fileexist3, fileexist4, fileexist5, fileexist6
  integer :: numgeombodies,iunit,ibody

  INTEGER, POINTER, DIMENSION(:) :: numnp_glb,numel_glb,nummat_glb

  INTEGER :: ielemtest, icnt, inodetest, minel, maxel,maxnd, minnd, iel10

  INTEGER, POINTER, DIMENSION(:) :: epart_loc 

  INTEGER, POINTER, dimension(:) :: procdist
  
  TYPE(link_ptr_type)  :: link

  INTEGER :: accumel, accumnd, accumnd_old, accumel_old

  INTEGER :: ierror, iv, kc

  size1 = 0.d0
  size2 = 0.d0
  size3 = 0.d0
  size4 = 0.d0
  size5 = 0.d0
  size6 = 0.d0

  ierror = 0 
  
  hex8  = .false.
  tet4  = .false.
  tet10 = .false.
  overlaymesh = .false.
  interactmesh = .false.
  
  print*,'MESH OPTION:'
  print*,'  READING PATRAN MESH'
  print*,' '
  
  itmp   = 0
  maxnumbc_str = 0
  maxnumbc_mm  = 0
  maxnumbc_th  = 0

  numndhistory = 0
  
  ALLOCATE(numelhex2d(1:5,1:nprocs))
  ALLOCATE(numeltet2d(1:5,1:nprocs))
  
  numelhex2d(:,:) = 0
  numeltet2d(:,:) = 0

  fileexist1 = .false.
  fileexist2 = .false.
  fileexist3 = .false.
  fileexist4 = .false.
  fileexist5 = .false.
  fileexist6 = .false.

  INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.1.out', exist=fileexist1)
  INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.2.out', exist=fileexist2)
  INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.3.out', exist=fileexist3)
  INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.4.out', exist=fileexist4)
  INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.5.out', exist=fileexist5)
  INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.6.out', exist=fileexist6)

  numgeombodies = 1

  IF(fileexist6)THEN
     numgeombodies = 6
     OPEN(106,file='Modin/'//prefx(1:prefx_lngth)//'.6.out')
     OPEN(105,file='Modin/'//prefx(1:prefx_lngth)//'.5.out')
     OPEN(104,file='Modin/'//prefx(1:prefx_lngth)//'.4.out')
     OPEN(103,file='Modin/'//prefx(1:prefx_lngth)//'.3.out')
     OPEN(102,file='Modin/'//prefx(1:prefx_lngth)//'.2.out')
     OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.1.out')
  ELSE IF (fileexist5)THEN
     numgeombodies = 5
     OPEN(105,file='Modin/'//prefx(1:prefx_lngth)//'.5.out')
     OPEN(104,file='Modin/'//prefx(1:prefx_lngth)//'.4.out')
     OPEN(103,file='Modin/'//prefx(1:prefx_lngth)//'.3.out')
     OPEN(102,file='Modin/'//prefx(1:prefx_lngth)//'.2.out')
     OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.1.out')
  ELSE IF (fileexist4)THEN
     numgeombodies = 4
     OPEN(104,file='Modin/'//prefx(1:prefx_lngth)//'.4.out')
     OPEN(103,file='Modin/'//prefx(1:prefx_lngth)//'.3.out')
     OPEN(102,file='Modin/'//prefx(1:prefx_lngth)//'.2.out')
     OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.1.out')
  ELSE IF (fileexist3)THEN
     numgeombodies = 3
     OPEN(103,file='Modin/'//prefx(1:prefx_lngth)//'.3.out')
     OPEN(102,file='Modin/'//prefx(1:prefx_lngth)//'.2.out')
     OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.1.out')
  ELSE IF (fileexist2) THEN
     numgeombodies = 2
     OPEN(102,file='Modin/'//prefx(1:prefx_lngth)//'.2.out')
     OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.1.out')
  ELSE IF (fileexist1) THEN
     OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.1.out')
  ELSE
     fileexist = .false.
     
     INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.pat', exist=fileexist)
     
     IF (.NOT. fileexist)THEN
        print*,'File '//'Modin/'//prefx(1:prefx_lngth)//'.pat NOT found'
        print*,'...Trying '//'Modin/'//prefx(1:prefx_lngth)//'.out'
        
        fileexist = .false.
        INQUIRE (file='Modin/'//prefx(1:prefx_lngth)//'.out', exist=fileexist)
     
        IF (.NOT. fileexist)THEN
           print*,'ERROR: File '//'Modin/'//prefx(1:prefx_lngth)//'.out NOT found'
           print*,'ERROR: No Patran Neutral file Found, stopping'
           stop
        ELSE
           print*,'File '//'Modin/'//prefx(1:prefx_lngth)//'.out found'

           OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.out', form='formatted')
        ENDIF
     ELSE
        OPEN(101,file='Modin/'//prefx(1:prefx_lngth)//'.pat', form='formatted')
     ENDIF

  ENDIF

  print*,'No. of Mesh Segment Files Found =', numgeombodies    

  ALLOCATE(numnp_glb(1:numgeombodies),numel_glb(1:numgeombodies),nummat_glb(1:numgeombodies))

  DO ibody = 1, numgeombodies
     iunit = 100 + ibody


!
! - Packet Type 25: Title Card, Format(i2,8i8)
     READ(iunit,'()')
     READ(iunit,'()')

! -  Packet Type 26: Summary Data, Format(i2,8i8)
! -      26 ID IV KC N1 N2 N3 N4 N5
! -  N1 = number of nodes
! -  N2 = number of elements
! -  N3 = number of materials
! -  N4 = number of Element Properties
! -  N5 = number of Coordinate Frames

     READ(iunit,*) i,i,i,i,numnp_glb(ibody),numel_glb(ibody),nummat_glb(ibody)
  
     IF(nummat_glb(ibody).LE.0) nummat_glb(ibody) = 1

! - part of title card
     READ(iunit,'()')

  ENDDO

  numnp_prmry = sum(numnp_glb)
  numelv_prmry = sum(numel_glb)
  nummat = maxval(nummat_glb)
  
  print*,' No. of elements = ', numelv_prmry
  print*,' No. of nodes = ', numnp_prmry
  print*,' No. of Materials = ', nummat
  print*,' No. of vertex=',numvertx

  inodetest = 0
  ielemtest = 0

  iaux = 0

  ALLOCATE(epart(1:numelv_prmry))

  accumnd = 0

  accumel = 0

  DO ibody = 1, numgeombodies
     iunit = 100 + ibody
     numnp4node = 0
!
! --  Read Nodal coordinates

     IF(ibody.EQ.1)THEN
        ALLOCATE(coor(1:3,1:numnp_prmry),lmelv_prmry(1:numvertx,1:numelv_prmry))
        
!!$        IF(numvertx.NE.4)THEN
!!$           ALLOCATE(lm4node( 1:4,1:numelv_prmry))
!!$           ALLOCATE(imap(1:numnp_prmry))
!!$           imap(:) = 0
!!$           numnp4node = 0
!!$        ELSE
!!$           numnp4node = numnp_prmry
!!$        ENDIF
     ENDIF

     iel10 = 0

!      ConvertUnit = .0254d0 ! for titan, units in engles
      
     DO i = 1, numnp_glb(ibody)

! - Type 01: Node Data,  format(i2,8i8)
!        1 ID IV KC
!    ID = node id
!    IV = 0 n/a
!    KC = 2
        READ(iunit,*) j,id
        id = id + accumnd
! - format(3e16.9)
! - caresian coordinate of Nodes x, y, z
        READ(iunit,'(3e16.9)') coor(1:3,id)
        coor(1:3,id) = coor(1:3,id)*convertunit
! - ICF GTYPE NDF CONFIG CID PSPC
        READ(iunit,'()')

        inodetest = inodetest + 1
     ENDDO
     accumnd_old = accumnd
     accumnd =  accumnd + numnp_glb(ibody)
!
! --  Read element connectivity array
     IF(ibody.EQ.1)THEN
        ALLOCATE(matid(1:numelv_prmry))
        ALLOCATE(eltypeid(1:numelv_prmry))
     ENDIF
! flages for element type id
!
!    4 = 4 node tetrahedral
!    8 = 8 node hexahedral
!   10 = 10 node tetrahedral
!   12 = 12 node wedge* 
!   15 = 15 node wedge*
!   20 = 20 node hexahdral*
!
!    * Not yet implemented

     minel = 1e9
     maxel = 0

     element: DO i = 1, numel_glb(ibody)

! - Type 02: Element Data, format(i2,8i8)
!        2 ID IV KC N1
!     ID = element ID
!     IV = Shape (5 = tet)
        READ(iunit,*) j, id
        ielemtest = ielemtest + 1
        accumel_old = accumel
        id = accumel + id


!     NODES CONFIG CID CEID 
        READ(iunit,'(3i8)') numvertx,j, matid(id)
        IF(numvertx.NE.4)THEN
           IF(numvertx.NE.10)THEN
              IF(numvertx.NE.8)THEN
                 print*,'ERROR: FOUND UNSUPPORTED ELEMENT TYPE'
                 print*,'Found element with ', numvertx,' nodes '
                 print*,'Stopping'
                 stop
              ENDIF
           ENDIF
        ENDIF

        eltypeid(id) = numvertx
     
        IF(matid(id).LE.0) matid(id) = 1
!     LNODES
!        LNODES = Element corner nodes followed by additional nodes
        READ(iunit,*) lmelv_prmry(1:numvertx,id) 
        lmelv_prmry(1:numvertx,id) = lmelv_prmry(1:numvertx,id) + accumnd_old
        
        iel10 = iel10 + 1
        IF(numvertx.EQ.10.AND.i.EQ.1.AND.nprocs.GT.1)THEN
           ALLOCATE(lm4node( 1:4,1:numel_glb(ibody)))
           IF(ibody.EQ.1) ALLOCATE(imap(1:numnp_prmry))
           imap(:) = 0
           print*,'Re-setting imap'
        ENDIF
        IF(numvertx.EQ.4.AND.i.EQ.1.AND.nprocs.GT.1)THEN
           ALLOCATE(lm4node( 1:4,1:numel_glb(ibody)))
           IF(ibody.EQ.1) ALLOCATE(imap(1:numnp_prmry))
           imap(:) = 0
           print*,'Re-setting imap'
        ENDIF
        IF(numvertx.EQ.10.AND.nprocs.GT.1)THEN
           DO j = 1, 4
              IF(imap(lmelv_prmry(j,id)).EQ.0)THEN
                 numnp4node = numnp4node + 1
                 imap(lmelv_prmry(j,id)) = numnp4node
              ENDIF
              iaux = iaux + 1
              lm4node(j,iel10) = imap(lmelv_prmry(j,id))
           ENDDO

        ELSE IF(numvertx.EQ.4.AND.nprocs.GT.1)THEN
           DO j = 1, 4
              IF(imap(lmelv_prmry(j,id)).EQ.0)THEN
                 numnp4node = numnp4node + 1
                 imap(lmelv_prmry(j,id)) = numnp4node
              ENDIF
              iaux = iaux + 1
              lm4node(j,iel10) = imap(lmelv_prmry(j,id))
           ENDDO
        ENDIF
        maxel = sum(numel_glb(1:ibody)) ! SUM(          MAX(MaxEl,id)
        minel = accumel_old + 1 !MIN(,id)


!
! -- Find the size of the smallest element
!
        xx = coor(1,lmelv_prmry(1,id)) - coor(1,lmelv_prmry(2,id))
        yy = coor(2,lmelv_prmry(1,id)) - coor(2,lmelv_prmry(2,id))
        zz = coor(3,lmelv_prmry(1,id)) - coor(3,lmelv_prmry(2,id))
        size1 = sqrt(xx*xx+yy*yy+zz*zz)
        xx = coor(1,lmelv_prmry(2,id)) - coor(1,lmelv_prmry(3,id))
        yy = coor(2,lmelv_prmry(2,id)) - coor(2,lmelv_prmry(3,id))
        zz = coor(3,lmelv_prmry(2,id)) - coor(3,lmelv_prmry(3,id))
        size2 = sqrt(xx*xx+yy*yy+zz*zz)
        xx = coor(1,lmelv_prmry(3,id)) - coor(1,lmelv_prmry(1,id))
        yy = coor(2,lmelv_prmry(3,id)) - coor(2,lmelv_prmry(1,id))
        zz = coor(3,lmelv_prmry(3,id)) - coor(3,lmelv_prmry(1,id))
        size3 = sqrt(xx*xx+yy*yy+zz*zz)
        xx = coor(1,lmelv_prmry(4,id)) - coor(1,lmelv_prmry(1,id))
        yy = coor(2,lmelv_prmry(4,id)) - coor(2,lmelv_prmry(1,id))
        zz = coor(3,lmelv_prmry(4,id)) - coor(3,lmelv_prmry(1,id))
        size4 = sqrt(xx*xx+yy*yy+zz*zz)
        xx = coor(1,lmelv_prmry(4,id)) - coor(1,lmelv_prmry(2,id))
        yy = coor(2,lmelv_prmry(4,id)) - coor(2,lmelv_prmry(2,id))
        zz = coor(3,lmelv_prmry(4,id)) - coor(3,lmelv_prmry(2,id))
        size5 = sqrt(xx*xx+yy*yy+zz*zz)
        xx = coor(1,lmelv_prmry(4,id)) - coor(1,lmelv_prmry(3,id))
        yy = coor(2,lmelv_prmry(4,id)) - coor(2,lmelv_prmry(3,id))
        zz = coor(3,lmelv_prmry(4,id)) - coor(3,lmelv_prmry(3,id))
        size6 = sqrt(xx*xx+yy*yy+zz*zz)
        dhmin = min(size1,size2,size3,size4,size5,size6,dhmin)
     ENDDO element

     print*,'Done reading Element'
     accumel = accumel + numel_glb(ibody)

     IF(nprocs.GT.1)THEN
        IF(numvertx.EQ.10)THEN
           ALLOCATE(npart(1:numnp4node),stat=ierror)
           IF(ierror .NE.0)THEN
              ! Space for npart could not be allcoated
              print*,'Program could not allocate space for npart'
              print*,'Dimensions 1 to ', numnp4node
              stop
           ENDIF
              
!!$           IF(ibody.EQ.4) THEN
!!$              WRITE(44,*) NumEL_glb(ibody),numnp4node,nprocs
!!$              MaxNd = MAXVAL(lm4node)
!!$              MinNd = MINVAL(lm4node)
!!$              PRINT*,'MaxNd =', MaxNd
!!$              PRINT*,'MinND =', MinNd
!!$              PRINT*,'MinNDlocation',MINLOC(lm4node)
!!$              PRINT*,'Maximum renumbered node =', numnp4node
!!$              PRINT*,'Number of elements in section',NumEL_glb(ibody)
!!$              PRINT*,'MinEl:MaxEl',MinEl,MaxEl
!!$              DO i = 1, NumEL_glb(ibody)
!!$                 WRITE(44,*) lm4node(:,i)
!!$                 IF(ANY(lm4node(:,i) .EQ. 0)) THEN
!!$                    PRINT*,'zero connectivity'
!!$                    STOP
!!$                 ENDIF
!!$              ENDDO
!!$              CLOSE(44)
!!$           ENDIF
!           PRINT*,'METIS for elements ',MinEL,':',MaxEl


!!$           CALL METIS_PartMeshDual(NumEL_glb(ibody),numnp4node, &
!!$                lm4node,2,1,nprocs,edgecut,epart(MinEl:MaxEl),npart)

           print*,'METIS for elements ',minel,':',maxel,numel_glb(ibody)

           ALLOCATE(epart_loc(1:numel_glb(ibody)),stat=ierror)

            IF(ierror .NE.0)THEN
              ! Space for npart could not be allcoated
              print*,'Program could not allocate space for epart_loc'
              print*,'Dimensions 1 to ', numel_glb(ibody)
              stop
           ENDIF

           print*,'CALLING METIS'

           CALL metis_partmeshnodal(numel_glb(ibody),numnp4node, &
                lm4node,2,1,nprocs,edgecut,epart_loc,npart)

           epart(minel:maxel) = epart_loc(1:numel_glb(ibody))
           print*,minel,maxel,1,numel_glb(ibody)
           deallocate(epart_loc)

           print*,'Finished Metis'
           DEALLOCATE(npart)
           DEALLOCATE(lm4node)

        ELSE IF(numvertx.EQ.4)THEN
           ALLOCATE(npart(1:numnp_glb(ibody)),stat=ierror)
           IF(ierror .NE.0)THEN
              ! Space for npart could not be allcoated
              print*,'Program could not allocate space for npart'
              print*,'Dimensions 1 to ', numnp4node
              stop
           ENDIF

!!$           CALL METIS_PartMeshDual(NumEL_glb(ibody),numnp4node, &
!!$                lm4node,2,1,nprocs,edgecut,epart(MinEl:MaxEl),npart)

           print*,'METIS for elements ',minel,':',maxel,numel_glb(ibody)

           ALLOCATE(epart_loc(1:numel_glb(ibody)),stat=ierror)

            IF(ierror .NE.0)THEN
              ! Space for npart could not be allcoated
              print*,'Program could not allocate space for epart_loc'
              print*,'Dimensions 1 to ', numel_glb(ibody)
              stop
           ENDIF

           print*,'CALLING METIS',numel_glb(ibody),numnp_glb(ibody)

           CALL metis_partmeshnodal(numel_glb(ibody),numnp_glb(ibody), &
                lm4node,2,1,nprocs,edgecut,epart_loc,npart)

           epart(minel:maxel) = epart_loc(1:numel_glb(ibody))
           print*,minel,maxel,1,numel_glb(ibody)
           deallocate(epart_loc)

           print*,'Finished Metis'
           DEALLOCATE(npart)
           DEALLOCATE(lm4node)


        ENDIF
     ELSE
        epart = 1
     ENDIF
  ENDDO

!!$
!!$  OPEN(45,file='LoadBalanceStats.out')
!!$  ALLOCATE(ProcDist(1:nprocs))
!!$  ProcDist(:) = 0
!!$  DO i = 1, numelv_prmry
!!$     IF(epart(i).EQ.0)THEN
!!$        print*,'epart = 0',i
!!$        stop
!!$     endif
!!$     ProcDist(epart(i)) = ProcDist(epart(i)) + 1
!!$  ENDDO
!!$
!!$  DO i=1, nprocs
!!$     write(45,*) i, ProcDist(i)
!!$  enddo
!!$  PRINT*,'*** Load Balancing Stats ***'
!!$  PRINT*,'Maximum number of elements on a processor =', MAXVAL(ProcDist(:))
!!$  PRINT*,'Minimum number of elements on a processor =', MINVAL(ProcDist(:))
!!$  close(45)
!!$  DEALLOCATE(ProcDist)



  IF(inodetest.NE.numnp_prmry)THEN
     print*,'ERROR number of combined nodes not equal total nodes'
     print*,'Sum of Nodes =', numnp_prmry
     print*,'Actual Number Read =', inodetest
     stop
  ENDIF

  IF(ielemtest.NE.numelv_prmry)THEN
     print*,'ERROR number of combined elements not equal total elements'
     print*,'Sum of Elements =', numelv_prmry
     print*,'Actual Number Read =', ielemtest
     stop
  ENDIF

! -- Partition the finite element mesh

  print*,' REQUESTED NUMBER OF PARTITIONS =',nprocs

!
!     Call METIS using the partition the mesh
!


  print*,'CALLING METIS'

  IF(nprocs.GT.1)THEN
!-ALT       CALL METIS_PartMeshNodal(numel_old_z,nn,
!-ALT$          elmnts,1,1,nprocs,edgecut,epart_p,npart_p)
     IF(numvertx.EQ.8)THEN
        ALLOCATE(npart(1:numnp_prmry))
        CALL metis_partmeshdual(numelv_prmry,numnp_prmry, &
             lmelv_prmry,3,1,nprocs,edgecut,epart,npart)
        DEALLOCATE(npart)
     ENDIF
!!$     ELSE IF(numvertx.EQ.4)THEN
!!$        ALLOCATE(npart(1:numnp_prmry))
!!$        CALL METIS_PartMeshDual(numelv_prmry,numnp_prmry, &
!!$             lmelv_prmry,2,1,nprocs,edgecut,epart,npart)
!!$        DEALLOCATE(npart)
!!$     ENDIF
  ELSE
     epart(:) = 1
  ENDIF

  print*,' Called METIS'
  OPEN(45,file='LoadBalanceStats.out')
  ALLOCATE(procdist(1:nprocs))
  procdist(:) = 0
  DO i = 1, numelv_prmry
     IF(epart(i).EQ.0)THEN
        print*,'epart = 0',i
        stop
     endif
     procdist(epart(i)) = procdist(epart(i)) + 1
  ENDDO

  DO i=1, nprocs
     write(45,*) i, procdist(i)
  enddo
  print*,'*** Load Balancing Stats ***'
  print*,'Maximum number of elements on a processor =', maxval(procdist(:))
  print*,'Minimum number of elements on a processor =', minval(procdist(:))
  close(45)
  DEALLOCATE(procdist)


! create communication file ! Note: should be able to remove restriction

  ALLOCATE(procndlist(1:numnp_prmry,1:maxnumberofprocstosharenode) )
  ALLOCATE(numprocpernd(1:numnp_prmry))
  ALLOCATE(numelperproc(1:nprocs))
  ALLOCATE(numndperproc(1:nprocs))
  
  numndperproc(:) = 0
  numelperproc(:) = 0

  print*,'Calling NewCommList'
  CALL newcommlist(numelv_prmry, numnp_prmry, nprocs, numvertx,  &
       lmelv_prmry, epart,numprocpernd,procndlist,maxnumberofprocstosharenode,numelperproc, &
       numndperproc)

  print*,'  .. completed'

! -- create a view of the partitioned mesh
! -- PMVIS software
! -- Command: pmvis.sgi.bin -n pmvis.nod -c pmvis.ele -p pmvis.part -o 1

  iopmvis = 0
  
  IF(iopmvis.EQ.1)THEN
     
     OPEN(13,file='pmvis.nod')
     DO i = 1, numnp_prmry
        WRITE(13,'(3e16.9)') coor(1:3,i)
     ENDDO
     CLOSE(13)
     OPEN(13,file='pmvis.ele')
     DO i = 1, numelv_prmry
        WRITE(13,'(4i10)') lmelv_prmry(1:4,i)
     ENDDO
     CLOSE(13)
     
     OPEN(13,file='pmvis.part')
     DO i = 1, numelv_prmry
        WRITE(13,'(1i10)') epart(i)
     ENDDO
     CLOSE(13)
  ENDIF
!
!
! -- Continue to read in the rest of the input parameters

  numbc_prmry = 0
  numbc_prmry_mm = 0
  numbc_prmry_ht = 0
  numserbc = 0
  
!!$  ALLOCATE(numel_2d(1:nprocs,1:2))
!!$  numel_2d(:,:) = 0
  

  ALLOCATE(indsburnflg(1:numnp_prmry))
  
  indsburnflg(:) = 0

! Initialize list for boundary conditions

  nullify(bc_structural_head,bc_structural_tail)
  nullify(bc_meshmotion_head,bc_meshmotion_tail)
  nullify(bc_thermal_head,bc_thermal_tail)
  
  nullify(surfmesh_tri3_ov1_head,surfmesh_tri3_ov1_tail)
  nullify(surfmesh_tri3_ov2_head,surfmesh_tri3_ov2_tail)
  nullify(surfmesh_tri3_s_head,surfmesh_tri3_s_tail)
  nullify(surfmesh_tri3_sf_head,surfmesh_tri3_sf_tail)
  nullify(surfmesh_tri6_s_head,surfmesh_tri6_s_tail)
  nullify(surfmesh_tri6_sf_head,surfmesh_tri6_sf_tail)
  nullify(surfmesh_hex8_s_head,surfmesh_hex8_s_tail)
  nullify(surfmesh_hex8_sf_head,surfmesh_hex8_sf_tail)

! number of boundary conditions per processors

  ALLOCATE(numbc_structural(1:nprocs))
  ALLOCATE(numbc_meshmotion(1:nprocs))
  ALLOCATE(numbc_thermal(1:nprocs))
  ALLOCATE(numbc_flag(1:nprocs))
  numbc_structural(:) = 0
  numbc_meshmotion(:) = 0
  numbc_thermal(:)    = 0 
  numbc_flag(:)       = 0
  ALLOCATE(bc_flag(1:3,1:numnp_prmry))
 
  bc_flag(:,:) = 0
  

  print*,'FINISH READING PATRAN FILE'
  print*,'Reading Boundary Conditions'
  accumel = 0
  accumnd = 0
  DO ibody = 1, numgeombodies
     iunit = 100 + ibody


     DO
        READ(iunit,*) itype,id,iv,kc
        IF(itype.EQ.99) EXIT
        
        IF(itype.EQ.4)THEN
           READ(iunit,'()')

! Marks of IntFaceFlag :
!
!   SolidFluid Interface = 1
!   NotASolidFluid Interface = 0
         
        ELSE IF(itype.EQ. 3) THEN
           DO i = 1, 20
              READ(iunit,'()')
           ENDDO
        ELSE IF(itype .EQ. 6)THEN      ! pressure loading
           id = id + accumel
           ip = epart(id)
           READ(iunit,'(i1,i1,i1,i6,8i1)')i,i,i,i,n1,n2,n3,n4,n5,n6,n7,n8
           READ(iunit,*) press
           intfaceflag = abs(int(press))

           IF(.NOT.(interactmesh))THEN
              IF(intfaceflag.GE.0.AND.intfaceflag.LE.2)  interactmesh = .true.
           ENDIF

! tet elements
           
           IF(n1.EQ.1.AND.n2.EQ.1.AND.n3.EQ.1   &
                .AND.n4.EQ.0.AND.n5.EQ.0.AND.n6.EQ.0.AND.n7.EQ.0.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.4)THEN
                 nnd1 = lmelv_prmry(1,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(2,id)
                 !itmp = itmp + 1
                 ifaceid = 1
              ELSE
                 nnd1 = lmelv_prmry(1,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(2,id)
                 nnd4 = lmelv_prmry(7,id)
                 nnd5 = lmelv_prmry(6,id)
                 nnd6 = lmelv_prmry(5,id)
                 ifaceid = 1
              ENDIF
           ELSE IF(n1.EQ.1.AND.n2.EQ.1.AND.n4.EQ.1  &
                .AND.n3.EQ.0.AND.n5.EQ.0.AND.n6.EQ.0.AND.n7.EQ.0.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.4)THEN
                 nnd1 = lmelv_prmry(1,id)
                 nnd2 = lmelv_prmry(2,id)
                 nnd3 = lmelv_prmry(4,id)
                 !itmp = itmp + 1
                 ifaceid = 2
              ELSE
                 nnd1 = lmelv_prmry(1,id)
                 nnd2 = lmelv_prmry(2,id)
                 nnd3 = lmelv_prmry(4,id)
                 nnd4 = lmelv_prmry(5,id)
                 nnd5 = lmelv_prmry(9,id)
                 nnd6 = lmelv_prmry(8,id)
                 ifaceid = 2
              ENDIF
           ELSE IF(n2.EQ.1.AND.n3.EQ.1.AND.n4.EQ.1  &
                .AND.n1.EQ.0.AND.n5.EQ.0.AND.n6.EQ.0.AND.n7.EQ.0.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.4)THEN
                 nnd1 = lmelv_prmry(2,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(4,id)
                 !itmp = itmp + 1
                 ifaceid = 3
              ELSE
                 nnd1 = lmelv_prmry(2,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(4,id)
                 nnd4 = lmelv_prmry(6,id)
                 nnd5 = lmelv_prmry(10,id)
                 nnd6  = lmelv_prmry(9,id)
                 ifaceid = 3
              ENDIF
           ELSE IF(n1.EQ.1.AND.n3.EQ.1.AND.n4.EQ.1  &
                .AND.n2.EQ.0.AND.n5.EQ.0.AND.n6.EQ.0.AND.n7.EQ.0.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.4)THEN
                 nnd1 = lmelv_prmry(4,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(1,id)
                 itmp = itmp + 1
                 ifaceid = 4
              ELSE
                 nnd1 = lmelv_prmry(4,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(1,id)
                 nnd4 = lmelv_prmry(10,id)
                 nnd5 = lmelv_prmry(7,id)
                 nnd6 = lmelv_prmry(8,id)
                 ifaceid = 4
              ENDIF
! Hex elements
           ELSE IF(n1.EQ.1.AND.n2.EQ.1.AND.n3.EQ.1.AND.n4.EQ.1  &
                .AND.n5.EQ.0.AND.n6.EQ.0.AND.n7.EQ.0.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.8)THEN
                 nnd1 = lmelv_prmry(4,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(2,id)
                 nnd4 = lmelv_prmry(1,id)
                 itmp = itmp + 1
                 ifaceid = 1
              ENDIF
           ELSE IF(n1.EQ.0.AND.n2.EQ.0.AND.n3.EQ.0.AND.n4.EQ.0  &
                .AND.n5.EQ.1.AND.n6.EQ.1.AND.n7.EQ.1.AND.n8.EQ.1)THEN
              IF(eltypeid(id).EQ.8)THEN
                 nnd1 = lmelv_prmry(5,id)
                 nnd2 = lmelv_prmry(6,id)
                 nnd3 = lmelv_prmry(7,id)
                 nnd4 = lmelv_prmry(8,id)
                 itmp = itmp + 1
                 ifaceid = 2
              ENDIF
           ELSE IF(n1.EQ.1.AND.n2.EQ.1.AND.n3.EQ.0.AND.n4.EQ.0  &
                .AND.n5.EQ.1.AND.n6.EQ.1.AND.n7.EQ.0.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.8)THEN
                 nnd1 = lmelv_prmry(1,id)
                 nnd2 = lmelv_prmry(2,id)
                 nnd3 = lmelv_prmry(6,id)
                 nnd4 = lmelv_prmry(5,id)
                 itmp = itmp + 1
                 ifaceid = 3
              ENDIF
           ELSE IF(n1.EQ.0.AND.n2.EQ.1.AND.n3.EQ.1.AND.n4.EQ.0  &
                .AND.n5.EQ.0.AND.n6.EQ.1.AND.n7.EQ.1.AND.n8.EQ.0)THEN
              IF(eltypeid(id).EQ.8)THEN
                 nnd1 = lmelv_prmry(2,id)
                 nnd2 = lmelv_prmry(3,id)
                 nnd3 = lmelv_prmry(7,id)
                 nnd4 = lmelv_prmry(6,id)
                 itmp = itmp + 1
                 ifaceid = 4
              ENDIF
           ELSE IF(n1.EQ.0.AND.n2.EQ.0.AND.n3.EQ.1.AND.n4.EQ.1  &
                .AND.n5.EQ.0.AND.n6.EQ.0.AND.n7.EQ.1.AND.n8.EQ.1)THEN
              IF(eltypeid(id).EQ.8)THEN
                 nnd1 = lmelv_prmry(3,id)
                 nnd2 = lmelv_prmry(4,id)
                 nnd3 = lmelv_prmry(8,id)
                 nnd4 = lmelv_prmry(7,id)
                 itmp = itmp + 1
                 ifaceid = 5
              ENDIF
           ELSE IF(n1.EQ.1.AND.n2.EQ.0.AND.n3.EQ.0.AND.n4.EQ.1  &
                .AND.n5.EQ.1.AND.n6.EQ.0.AND.n7.EQ.0.AND.n8.EQ.1)THEN
              IF(eltypeid(id).EQ.8)THEN
                 nnd1 = lmelv_prmry(5,id)
                 nnd2 = lmelv_prmry(8,id)
                 nnd3 = lmelv_prmry(4,id)
                 nnd4 = lmelv_prmry(1,id)
                 itmp = itmp + 1
                 ifaceid = 6
              ENDIF
           ELSE
              print*,'Error in pressure face numbering'
              stop
           ENDIF
        
! BC flags
! ignitable (=1)
! non-ignitable (=0)
! noninteracting (=2)


! 3-node triangle, (via 4-node tetrahedral)
        
           IF(eltypeid(id).EQ.4)THEN

              ALLOCATE(surfmesh_tri3_item)
              
              surfmesh_tri3_item%ElemData(1) = nnd1
              surfmesh_tri3_item%ElemData(2) = nnd2
              surfmesh_tri3_item%ElemData(3) = nnd3
              surfmesh_tri3_item%ElemData(4) = id
              
              IF(intfaceflag.EQ.1)THEN ! Fluid Structure interaction, with burning
                 CALL add_surfmesh_tri3(surfmesh_tri3_item, surfmesh_tri3_sf_head, surfmesh_tri3_sf_tail)
                 numeltet2d(1,ip) = numeltet2d(1,ip) + 1
                 
              ELSE IF(intfaceflag.EQ.2)THEN ! fluid Structure interaction, without burning
                 CALL add_surfmesh_tri3(surfmesh_tri3_item, surfmesh_tri3_sf_nonignt_head, surfmesh_tri3_sf_nonignt_tail)
                 numeltet2d(3,ip) = numeltet2d(3,ip) + 1
              ELSE IF(intfaceflag.EQ.0)THEN ! No fluid structure interaction
                 CALL add_surfmesh_tri3(surfmesh_tri3_item, surfmesh_tri3_s_head, surfmesh_tri3_s_tail)
                 numeltet2d(2,ip) = numeltet2d(2,ip) + 1
              ELSE IF(intfaceflag.EQ.10)THEN ! No fluid structure interaction for non-matching meshes
                 overlaymesh = .true.
                 surfmesh_tri3_item%ElemData(5) =  ifaceid
                 CALL add_surfmesh_tri3(surfmesh_tri3_item, surfmesh_tri3_ov1_head, surfmesh_tri3_ov1_tail)
                 numeltet2d(4,ip) = numeltet2d(4,ip) + 1
              ELSE IF(intfaceflag.EQ.100)THEN ! No fluid structure interaction for non-matching meshes    
                 overlaymesh = .true.    
                 surfmesh_tri3_item%ElemData(5) =  ifaceid         
                 CALL add_surfmesh_tri3(surfmesh_tri3_item, surfmesh_tri3_ov2_head, surfmesh_tri3_ov2_tail)
                 numeltet2d(5,ip) = numeltet2d(5,ip) + 1

              ENDIF
              tet4 = .true.

! 6-node triangle, (via 10-node tetrahedral)

           ELSE IF(eltypeid(id).EQ.10)THEN
              
              ALLOCATE(surfmesh_tri6_item)
              
              surfmesh_tri6_item%ElemData(1) = nnd1
              surfmesh_tri6_item%ElemData(2) = nnd2
              surfmesh_tri6_item%ElemData(3) = nnd3
              surfmesh_tri6_item%ElemData(4) = nnd4
              surfmesh_tri6_item%ElemData(5) = nnd5
              surfmesh_tri6_item%ElemData(6) = nnd6
              surfmesh_tri6_item%ElemData(7) = id
              
              IF(intfaceflag.EQ.1)THEN
                 CALL add_surfmesh_tri6(surfmesh_tri6_item, surfmesh_tri6_sf_head, surfmesh_tri6_sf_tail)
                 numeltet2d(1,ip) = numeltet2d(1,ip) + 1
              ELSEIF(intfaceflag.EQ.2)THEN
                 CALL add_surfmesh_tri6(surfmesh_tri6_item, surfmesh_tri6_sf_nonignt_head, surfmesh_tri6_sf_nonignt_tail)
                 numeltet2d(3,ip) = numeltet2d(3,ip) + 1
              ELSE IF(intfaceflag.EQ.0) THEN
                 CALL add_surfmesh_tri6(surfmesh_tri6_item, surfmesh_tri6_s_head, surfmesh_tri6_s_tail)
                 numeltet2d(2,ip) = numeltet2d(2,ip) + 1
              ENDIF
              tet10 = .true.

! 4-node quad, (via 8-node hexahedral)
          
           ELSE IF(eltypeid(id).EQ.8)THEN
              
              ALLOCATE(surfmesh_hex8_item)
              
              surfmesh_hex8_item%ElemData(1) = nnd1
              surfmesh_hex8_item%ElemData(2) = nnd2
              surfmesh_hex8_item%ElemData(3) = nnd3
              surfmesh_hex8_item%ElemData(4) = nnd4
              surfmesh_hex8_item%ElemData(5) = id
              IF(intfaceflag.EQ.1)THEN
                 CALL add_surfmesh_hex8(surfmesh_hex8_item, surfmesh_hex8_sf_head, surfmesh_hex8_sf_tail)
                 numelhex2d(1,ip) = numelhex2d(1,ip) + 1
              ELSE IF(intfaceflag.EQ.2)THEN
                 CALL add_surfmesh_hex8(surfmesh_hex8_item, surfmesh_hex8_sf_nonignt_head, surfmesh_hex8_sf_nonignt_tail)
                 numelhex2d(3,ip) = numelhex2d(3,ip) + 1
              ELSE IF(intfaceflag.EQ.0) THEN
                 CALL add_surfmesh_hex8(surfmesh_hex8_item, surfmesh_hex8_s_head, surfmesh_hex8_s_tail)
                 numelhex2d(2,ip) = numelhex2d(2,ip) + 1
              ELSE
                 print*,'ERROR, interface flag type not found',intfaceflag
                 stop
              ENDIF
              
              hex8 = .true.
           ENDIF
        

!        numel_2d(ip,IntFaceFlag) = numel_2d(ip,IntFaceFlag) + 1
        

        ELSE IF(itype .EQ. 8)THEN ! boundary conditions
           id = id + accumnd
! 
!     CID = Coordinate frame ID
!     ICOMP = 6 displacement component flags (0 or 1)
!     note the flag is passed to node number one (i.e. fixed displacement)

           READ(iunit,'(I8,6I1)') iaux, disflag(1:6)
           numdisflag = sum(disflag(1:6))
           

           numserbc = numserbc + numdisflag

           READ(iunit,*) tmpval(1:numdisflag)
           
           icnt = 0
           do i = 1, 6
              IF(disflag(i).eq.1)THEN
                 icnt = icnt + 1
                 disflagvalue(i) = tmpval(icnt)
              endif
           enddo

           IF(icnt.NE.numdisflag)THEN
              print*,'ERROR in processing boundary conditions'
              stop
           ENDIF
        
! add to the number of boundary conditions per processor  

           IF(disflag(1).eq.1 )THEN ! structural boundary condition
              DO i = 1,numprocpernd(id)
                 procid = procndlist(id,i)
                 numbc_structural(procid) = numbc_structural(procid) + 1
                 IF(bc_flag(1,id).EQ.0) numbc_flag(procid) = numbc_flag(procid) + 1 
              ENDDO

              bc_flag(1,id) = bc_flag(1,id) + int(disflagvalue(1))
              
              maxnumbc_str = max(maxnumbc_str,int(disflagvalue(1)))
           endif

           IF(disflag(2).eq.1 )THEN ! temperature boundary condition
              DO i = 1,numprocpernd(id)
                 procid = procndlist(id,i)
                 numbc_thermal(procid) = numbc_thermal(procid) + 1
                 IF(bc_flag(2,id).EQ.0) numbc_flag(procid) = numbc_flag(procid) + 1 
              enddo
              
              bc_flag(2,id) = bc_flag(2,id) + 100*int(disflagvalue(2))
              
              maxnumbc_th = max(maxnumbc_th,int(disflagvalue(2)))
           endif


           IF(disflag(3).eq.1 )THEN ! mesh motion boundary condition
              DO i = 1,numprocpernd(id)
                 procid = procndlist(id,i)
                 numbc_meshmotion(procid) = numbc_meshmotion(procid) + 1
                 IF(bc_flag(3,id).EQ.0) numbc_flag(procid) = numbc_flag(procid) + 1 
              ENDDO
              bc_flag(3,id) = bc_flag(3,id) + 10000*int(disflagvalue(3))
              
              maxnumbc_mm = max(maxnumbc_mm,int(disflagvalue(3)))
           endif


        
        ELSE IF(itype .EQ. 7)THEN ! mesh motion boundary conditions
           id = id + accumnd
! 
!     CID = Coordinate frame ID
!     ICOMP = 6 displacement component flags (0 or 1)
!     note the flag is passed to node number one (i.e. fixed displacement)

           READ(iunit,'()')
           READ(iunit,*) value
           
!           ALLOCATE(BC_meshmotion_item)

!           BC_meshmotion_item%BC_nodeGlb = id
!           BC_meshmotion_item%BC_flagGlb = INT(value)
           
!           CALL add_BC(BC_meshmotion_item, BC_meshmotion_head, BC_meshmotion_tail)
           

! add to the number of boundary conditions per processor         

           DO i = 1,numprocpernd(id)
              procid = procndlist(id,i)
              numbc_meshmotion(procid) = numbc_meshmotion(procid) + 1
              IF(bc_flag(3,id).EQ.0) numbc_flag(procid) = numbc_flag(procid) + 1 
           ENDDO

           bc_flag(3,id) = bc_flag(3,id) + 10000*int(value)

           maxnumbc_mm = max(maxnumbc_mm,int(value))

        ELSE IF(itype .EQ. 5)THEN ! coordinate frames
           READ(iunit,'()')
           READ(iunit,'()')
           READ(iunit,'()')
           READ(iunit,'()')
        ELSE IF(itype .EQ. 10)THEN ! nodal temperature (corners)
           id = id + accumnd
           READ(iunit,*) value
           
!           ALLOCATE(BC_thermal_item)
           
!           BC_thermal_item%BC_nodeGlb = id
!           BC_thermal_item%BC_flagGlb = INT(value)
           
!           CALL add_BC(BC_thermal_item, BC_thermal_head, BC_thermal_tail)

! add to the number of boundary conditions per processor
  
           DO i = 1,numprocpernd(id)
              procid = procndlist(id,i)
              numbc_thermal(procid) = numbc_thermal(procid) + 1
              IF(bc_flag(2,id).EQ.0) numbc_flag(procid) = numbc_flag(procid) + 1 
           ENDDO

           bc_flag(2,id) = bc_flag(2,id) + 100*int(value)
           
           maxnumbc_th = max(maxnumbc_th,int(value))

        ELSE
           DO i = 1, kc
              READ(iunit,'()')
           ENDDO
        ENDIF
     ENDDO

     accumel =  sum(numel_glb(1:ibody))

     accumnd =  sum(numnp_glb(1:ibody))
  ENDDO

! if negative number then higher order element


  IF(tet4.AND..NOT.(tet10).AND..NOT.(hex8))THEN
     meshtype2d = 3
  ELSE IF (.NOT.(tet4).AND.tet10.AND..NOT.(hex8))THEN
     meshtype2d = 6
  ELSE IF (.NOT.(tet4).AND..NOT.(tet10).AND.hex8)THEN
     meshtype2d = 4
  ELSE
     meshtype2d = 5
  ENDIF

!  PRINT*,'Number of Surface elements, Interacting with burning =', NumEltet2D(1,:)
!  PRINT*,'Number of Surface elements, Interacting with Non-burning =', NumEltet2D(3,:)
!  PRINT*,'Number of Surface elements, Non-Interacting =', NumEltet2D(2,:)

  print*,'BOUNDARY CONDITIONS'
  print*,' Total of all types = ',numserbc
  print*,'  Number of structural = ',maxnumbc_str
  print*,'  Number of meshmotion = ',maxnumbc_mm
  print*,'  Number of thermal = ',maxnumbc_th

  IF(numvertx.EQ.4)THEN
     numvertx2d = 3
  ELSE IF(numvertx.EQ.10)THEN
     numvertx2d = 6
  ELSE IF(numvertx.EQ.8)THEN
     numvertx2d = 4
  ENDIF

  DO ibody = 1, numgeombodies
     iunit = 100 +  ibody 
     CLOSE(iunit)
  enddo
  
  RETURN
  
END SUBROUTINE read_patran