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