1 subroutine genbkb (ibksz) 2c 3c---------------------------------------------------------------------- 4c 5c This routine reads the boundary elements, reorders them and 6c generates traces for the gather/scatter operations. 7c 8c Zdenek Johan, Fall 1991. 9c---------------------------------------------------------------------- 10c 11 use dtnmod 12 use pointer_data 13 use phio 14 use iso_c_binding 15 include "common.h" 16 include "mpif.h" !Required to determine the max for itpblk 17 18 integer, target, allocatable :: ientp(:,:),iBCBtp(:,:) 19 real*8, target, allocatable :: BCBtp(:,:) 20 integer materb(ibksz) 21 integer, target :: intfromfile(50) ! integers read from headers 22 character*255 fname1 23 integer :: descriptor, descriptorG, GPID, color, nfiles, nfields 24 integer :: numparts, nppp, nprocs, writeLock 25 integer :: ierr_io, numprocs, itmp, itmp2 26 integer, target :: itpblktot,ierr 27 character*255 fname2 28 character(len=30) :: dataInt, dataDbl 29 dataInt = c_char_'integer'//c_null_char 30 dataDbl = c_char_'double'//c_null_char 31 32 nfiles = nsynciofiles 33 nfields = nsynciofieldsreadgeombc 34 numparts = numpe !This is the common settings. Beware if you try to compute several parts per process 35 nppp = numparts/numpe 36 ione=1 37 itwo=2 38 ieight=8 39 ieleven=11 40 itmp = int(log10(float(myrank+1)))+1 41 iel=1 42 itpblk=nelblb 43 44 ! Get the total number of different interior topologies in the whole domain. 45 ! Try to read from a field. If the field does not exist, scan the geombc file. 46 itpblktot=1 ! hardwired to monotopology for now 47 call phio_readheader(fhandle, 48 & c_char_'total number of boundary tpblocks' // char(0), 49 & c_loc(itpblktot), ione, dataInt, iotype) 50 51 if (itpblktot == -1) then 52 ! The field 'total number of different boundary tpblocks' was not found in the geombc file. 53 ! Scan all the geombc file for the 'connectivity interior' fields to get this information. 54 iblk=0 55 neltp=0 56 do while(neltp .ne. -1) 57 58 ! intfromfile is reinitialized to -1 every time. 59 ! If connectivity boundary@xxx is not found, then 60 ! readheader will return intfromfile unchanged 61 62 intfromfile(:)=-1 63 iblk = iblk+1 64 if(nfiles.gt.0)then 65 write (fname2,"('connectivity boundary',i1)") iblk 66 else 67 write (fname2,"('connectivity boundary linear tetrahedron')") 68 endif 69 70 call phio_readheader(fhandle, fname2 // char(0), 71 & c_loc(intfromfile), ieight, dataInt, iotype) 72 neltp = intfromfile(1) ! -1 if fname2 was not found, >=0 otherwise 73 end do 74 itpblktot = iblk-1 75 end if 76 77 if (myrank == 0) then 78 write(*,*) 'Number of boundary topologies: ',itpblktot 79 endif 80 81 nelblb=0 82 mattyp=0 83 ndofl = ndof 84 85 do iblk = 1, itpblktot 86 writeLock=0; 87 if(nfiles.gt.0)then 88 write (fname2,"('connectivity boundary',i1)") iblk 89 else 90 write (fname2,"('connectivity boundary linear tetrahedron')") 91 endif 92 93 ! Synchronization for performance monitoring, as some parts do not include some topologies 94 call MPI_Barrier(MPI_COMM_WORLD,ierr) 95 call phio_readheader(fhandle, fname2 // char(0), 96 & c_loc(intfromfile), ieight, dataInt, iotype) 97 neltp =intfromfile(1) 98 nenl =intfromfile(2) 99 ipordl=intfromfile(3) 100 nshl =intfromfile(4) 101 nshlb =intfromfile(5) 102 nenbl =intfromfile(6) 103 lcsyst=intfromfile(7) 104 numnbc=intfromfile(8) 105 106 allocate (ientp(neltp,nshl)) 107 allocate (iBCBtp(neltp,ndiBCB)) 108 allocate (BCBtp(neltp,ndBCB)) 109 iientpsiz=neltp*nshl 110 111 if (neltp==0) then 112 writeLock=1; 113 endif 114 115 call phio_readdatablock(fhandle, fname2 // char(0), 116 & c_loc(ientp),iientpsiz,dataInt,iotype) 117c 118c.... Read the boundary flux codes 119c 120 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 121 if(nfiles.gt.0)then 122 write (fname2,"('nbc codes',i1)") iblk 123 else 124 write (fname2,"('nbc codes linear tetrahedron')") 125 endif 126 127 call phio_readheader(fhandle, fname2 // char(0), 128 & c_loc(intfromfile), ieight, dataInt, iotype) 129 iiBCBtpsiz=neltp*ndiBCB 130 call phio_readdatablock(fhandle, fname2 // char(0), 131 & c_loc(iBCBtp),iiBCBtpsiz,dataInt,iotype) 132c 133c.... read the boundary condition data 134c 135 call MPI_BARRIER(MPI_COMM_WORLD, ierr) 136 if(nfiles.gt.0)then 137 write (fname2,"('nbc values',i1)") iblk 138 else 139 write (fname2,"('nbc values linear tetrahedron')") 140 endif 141 142 call phio_readheader(fhandle, fname2 // char(0), 143 & c_loc(intfromfile), ieight, dataInt, iotype) 144 BCBtp = zero 145 iBCBtpsiz=neltp*ndBCB 146 call phio_readdatablock(fhandle, fname2 // char(0), 147 & c_loc(BCBtp),iBCBtpsiz,dataDbl,iotype) 148c 149c This is a temporary fix until NSpre properly zeros this array where it 150c is not set. DEC has indigestion with these arrays though the 151c result is never used (never effects solution). 152c 153 if(writeLock==0) then 154 where(.not.btest(iBCBtp(:,1),0)) BCBtp(:,1)=zero 155 where(.not.btest(iBCBtp(:,1),1)) BCBtp(:,2)=zero 156 where(.not.btest(iBCBtp(:,1),3)) BCBtp(:,6)=zero 157 if(ndBCB.gt.6) then 158 do i=6,ndof 159 where(.not.btest(iBCBtp(:,1),i-1)) BCBtp(:,i+1)=zero 160 enddo 161 endif 162 where(.not.btest(iBCBtp(:,1),2)) 163 BCBtp(:,3)=zero 164 BCBtp(:,4)=zero 165 BCBtp(:,5)=zero 166 endwhere 167 168 do n=1,neltp,ibksz 169 nelblb=nelblb+1 170 npro= min(IBKSZ, neltp - n + 1) 171 lcblkb(1,nelblb) = iel 172 lcblkb(3,nelblb) = lcsyst 173 lcblkb(4,nelblb) = ipordl 174 lcblkb(5,nelblb) = nenl 175 lcblkb(6,nelblb) = nenbl 176 lcblkb(7,nelblb) = mattyp 177 lcblkb(8,nelblb) = ndofl 178 lcblkb(9,nelblb) = nshl 179 lcblkb(10,nelblb) = nshlb ! # of shape functions per elt 180c 181c.... save the element block 182c 183 n1=n 184 n2=n+npro-1 185 materb=1 ! all one material for now 186c 187c.... allocate memory for stack arrays 188c 189 allocate (mienb(nelblb)%p(npro,nshl)) 190 allocate (miBCB(nelblb)%p(npro,ndiBCB)) 191 allocate (mBCB(nelblb)%p(npro,nshlb,ndBCB)) 192 allocate (mmatb(nelblb)%p(npro)) 193c 194c.... save the boundary element block 195c 196 call gensvb (ientp(n1:n2,1:nshl), 197 & iBCBtp(n1:n2,:), BCBtp(n1:n2,:), 198 & materb, mienb(nelblb)%p, 199 & miBCB(nelblb)%p, mBCB(nelblb)%p, 200 & mmatb(nelblb)%p) 201 iel=iel+npro 202 enddo 203 endif 204 deallocate(ientp) 205 deallocate(iBCBtp) 206 deallocate(BCBtp) 207 208 enddo 209 lcblkb(1,nelblb+1) = iel 210 return 211c 212c.... end of file error handling 213c 214 911 call error ('genbcb ','end file',igeomBAK) 2151000 format(a80,//, 216 & ' B o u n d a r y E l e m e n t C o n n e c t i v i t y',//, 217 & ' Elem BC codes',/, 218 & ' Number C P V H ',5x,27('Node',i1,:,2x)) 2191100 format(2x,i5,2x,4i2,3x,27i7) 2202100 format(2x,i5,1p,1x,6e12.4) 221 end 222