| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | subroutine readat
|
|---|
| 3 | C
|
|---|
| 4 | C Read in data from preprocessor input file (.rea)
|
|---|
| 5 | C
|
|---|
| 6 | INCLUDE 'SIZE'
|
|---|
| 7 | INCLUDE 'INPUT'
|
|---|
| 8 | INCLUDE 'GEOM'
|
|---|
| 9 | INCLUDE 'PARALLEL'
|
|---|
| 10 | INCLUDE 'CTIMER'
|
|---|
| 11 |
|
|---|
| 12 | logical ifbswap,ifre2,parfound
|
|---|
| 13 | character*132 string
|
|---|
| 14 | integer idum(3*numsts+3)
|
|---|
| 15 |
|
|---|
| 16 | ierr = 0
|
|---|
| 17 | call flush_io
|
|---|
| 18 |
|
|---|
| 19 | ! check if new rea file version exists
|
|---|
| 20 | if(nid.eq.0) inquire(file=parfle, exist=parfound)
|
|---|
| 21 | call bcast(parfound,lsize)
|
|---|
| 22 | if (parfound) then
|
|---|
| 23 | if(nio.eq.0) write(6,'(A,A)') ' Reading ', parfle
|
|---|
| 24 | call readat_par
|
|---|
| 25 | goto 99
|
|---|
| 26 | endif
|
|---|
| 27 |
|
|---|
| 28 | etime0 = dnekclock_sync()
|
|---|
| 29 |
|
|---|
| 30 | if(nid.eq.0) then
|
|---|
| 31 | write(6,'(A,A)') ' Reading ', reafle
|
|---|
| 32 | open (unit=9,file=reafle,status='old', iostat=ierr)
|
|---|
| 33 | endif
|
|---|
| 34 |
|
|---|
| 35 | call bcast(ierr,isize)
|
|---|
| 36 | if (ierr .gt. 0) call exitti('Cannot open rea file!$',1)
|
|---|
| 37 |
|
|---|
| 38 | C Read parameters and logical flags
|
|---|
| 39 | call rdparam
|
|---|
| 40 | meshPartitioner=3 ! HYBRID (RSB+RCB)
|
|---|
| 41 |
|
|---|
| 42 | C Read Mesh Info
|
|---|
| 43 | if(nid.eq.0) then
|
|---|
| 44 | read(9,*) ! xfac,yfac,xzero,yzero
|
|---|
| 45 | read(9,*) ! dummy
|
|---|
| 46 | read(9,*) nelgs,ldimr,nelgv
|
|---|
| 47 | nelgt = abs(nelgs)
|
|---|
| 48 | endif
|
|---|
| 49 | call bcast(ldimr,ISIZE)
|
|---|
| 50 | call bcast(nelgs,ISIZE)
|
|---|
| 51 | call bcast(nelgv,ISIZE)
|
|---|
| 52 | call bcast(nelgt,ISIZE)
|
|---|
| 53 | ifre2 = .false.
|
|---|
| 54 | if (nelgs.lt.0) ifre2 = .true.
|
|---|
| 55 |
|
|---|
| 56 | call usrdat0
|
|---|
| 57 |
|
|---|
| 58 | if (nelgt.gt.350000 .and. .not.ifre2)
|
|---|
| 59 | $ call exitti('Problem size requires .re2!$',1)
|
|---|
| 60 |
|
|---|
| 61 | if (ifre2) call read_re2_hdr(ifbswap, .true.) ! rank0 will open and read
|
|---|
| 62 | call chk_nel ! make certain sufficient array sizes
|
|---|
| 63 |
|
|---|
| 64 | call mapelpr
|
|---|
| 65 |
|
|---|
| 66 | if (ifre2) then
|
|---|
| 67 | call read_re2_data(ifbswap, .true., .true., .true.)
|
|---|
| 68 | else
|
|---|
| 69 | maxrd = 32 ! max # procs to read at once
|
|---|
| 70 | mread = (np-1)/maxrd+1 ! mod param
|
|---|
| 71 | iread = 0 ! mod param
|
|---|
| 72 | x = 0
|
|---|
| 73 | do i=0,np-1,maxrd
|
|---|
| 74 | call nekgsync()
|
|---|
| 75 | if (mod(nid,mread).eq.iread) then
|
|---|
| 76 | if (nid.ne.0) then
|
|---|
| 77 | open(UNIT=9,FILE=REAFLE,STATUS='OLD')
|
|---|
| 78 | call cscan(string,'MESH DATA',9)
|
|---|
| 79 | read(9,*) string
|
|---|
| 80 | endif
|
|---|
| 81 | call rdmesh
|
|---|
| 82 | call rdcurve ! Curved side data
|
|---|
| 83 | call rdbdry ! Boundary Conditions
|
|---|
| 84 | if (nid.ne.0) close(unit=9)
|
|---|
| 85 | endif
|
|---|
| 86 | iread = iread + 1
|
|---|
| 87 | enddo
|
|---|
| 88 | endif
|
|---|
| 89 |
|
|---|
| 90 | C Read Restart options / Initial Conditions / Drive Force
|
|---|
| 91 | CALL RDICDF
|
|---|
| 92 | C Read materials property data
|
|---|
| 93 | CALL RDMATP
|
|---|
| 94 | C Read history data
|
|---|
| 95 | CALL RDHIST
|
|---|
| 96 | C Read output specs
|
|---|
| 97 | CALL RDOUT
|
|---|
| 98 | C Read objects
|
|---|
| 99 | CALL RDOBJ
|
|---|
| 100 |
|
|---|
| 101 | call nekgsync()
|
|---|
| 102 |
|
|---|
| 103 | C End of input data, close read file.
|
|---|
| 104 | if(nid.eq.0) then
|
|---|
| 105 | close(unit=9)
|
|---|
| 106 | call echopar
|
|---|
| 107 | write(6,'(A,g13.5,A,/)') ' done :: read .rea file ',
|
|---|
| 108 | $ dnekclock()-etime0,' sec'
|
|---|
| 109 | endif
|
|---|
| 110 |
|
|---|
| 111 | 99 call izero(boundaryID, size(boundaryID))
|
|---|
| 112 | call izero(boundaryIDt, size(boundaryIDt))
|
|---|
| 113 |
|
|---|
| 114 | ifld = 2
|
|---|
| 115 | if(ifflow) ifld = 1
|
|---|
| 116 | do iel = 1,nelv
|
|---|
| 117 | do ifc = 1,2*ndim
|
|---|
| 118 | boundaryID(ifc,iel) = bc(5,ifc,iel,ifld)
|
|---|
| 119 | enddo
|
|---|
| 120 | enddo
|
|---|
| 121 |
|
|---|
| 122 | ntmsh = 0
|
|---|
| 123 | do i=1,ldimt
|
|---|
| 124 | if(iftmsh(1+i)) ntmsh = ntmsh + 1
|
|---|
| 125 | enddo
|
|---|
| 126 |
|
|---|
| 127 | if (ntmsh.gt.0) then
|
|---|
| 128 | do iel = 1,nelt
|
|---|
| 129 | do ifc = 1,2*ndim
|
|---|
| 130 | boundaryIDt(ifc,iel) = bc(5,ifc,iel,2)
|
|---|
| 131 | enddo
|
|---|
| 132 | enddo
|
|---|
| 133 | endif
|
|---|
| 134 |
|
|---|
| 135 | return
|
|---|
| 136 | END
|
|---|
| 137 | c-----------------------------------------------------------------------
|
|---|
| 138 | subroutine vrdsmsh
|
|---|
| 139 | C
|
|---|
| 140 | C=====================================================================
|
|---|
| 141 | C Verify that mesh and dssum are properly defined by performing
|
|---|
| 142 | C a direct stiffness operation on the X,Y and Z coordinates.
|
|---|
| 143 | C Note that periodic faces are not checked here.
|
|---|
| 144 | C=====================================================================
|
|---|
| 145 | C
|
|---|
| 146 | INCLUDE 'SIZE'
|
|---|
| 147 | INCLUDE 'TOTAL'
|
|---|
| 148 | COMMON /SCRNS/ TA(LX1,LY1,LZ1,LELT),TB(LX1,LY1,LZ1,LELT)
|
|---|
| 149 | $ ,QMASK(LX1,LY1,LZ1,LELT),tmp(2)
|
|---|
| 150 | CHARACTER*3 CB
|
|---|
| 151 |
|
|---|
| 152 | c call vrdsmshx ! verify mesh topology
|
|---|
| 153 |
|
|---|
| 154 | if(nio.eq.0) write(*,*) 'verify mesh topology'
|
|---|
| 155 |
|
|---|
| 156 | IERR = 0
|
|---|
| 157 | EPS = 1.0e-04
|
|---|
| 158 | EPS = 1.0e-03
|
|---|
| 159 | IFIELD = 1
|
|---|
| 160 | if (nelgv.ne.nelgt .or. .not.ifflow) ifield = 2
|
|---|
| 161 | NXYZ1 = lx1*ly1*lz1
|
|---|
| 162 | NTOT = lx1*ly1*lz1*NELT
|
|---|
| 163 | NFACES = 2*ldim
|
|---|
| 164 |
|
|---|
| 165 | xmx = glmax(xm1,ntot)
|
|---|
| 166 | xmn = glmin(xm1,ntot)
|
|---|
| 167 | ymx = glmax(ym1,ntot)
|
|---|
| 168 | ymn = glmin(ym1,ntot)
|
|---|
| 169 | zmx = glmax(zm1,ntot)
|
|---|
| 170 | zmn = glmin(zm1,ntot)
|
|---|
| 171 | if (nio.eq.0) write(6,*) xmn,xmx,' Xrange'
|
|---|
| 172 | if (nio.eq.0) write(6,*) ymn,ymx,' Yrange'
|
|---|
| 173 | if (nio.eq.0) write(6,*) zmn,zmx,' Zrange'
|
|---|
| 174 | c return
|
|---|
| 175 | C
|
|---|
| 176 | C First check - use 1/Multiplicity
|
|---|
| 177 | C
|
|---|
| 178 | IF (IFHEAT) THEN
|
|---|
| 179 | CALL COPY(TA,TMULT,NTOT)
|
|---|
| 180 | ELSE
|
|---|
| 181 | CALL COPY(TA,VMULT,NTOT)
|
|---|
| 182 | ENDIF
|
|---|
| 183 | c
|
|---|
| 184 | c write(6,1)
|
|---|
| 185 | c $(nid,'tab4',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
|
|---|
| 186 | c 1 format(i3,a4,i3,16f5.2)
|
|---|
| 187 | c
|
|---|
| 188 | CALL DSSUM(TA,lx1,ly1,lz1)
|
|---|
| 189 | c
|
|---|
| 190 | c write(6,1)
|
|---|
| 191 | c $(nid,'taaf',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
|
|---|
| 192 | c
|
|---|
| 193 | CALL RONE (TB,NTOT)
|
|---|
| 194 | CALL SUB2 (TB,TA,NTOT)
|
|---|
| 195 | DO 1000 IE=1,NELT
|
|---|
| 196 | ieg=lglel(ie)
|
|---|
| 197 | DO 1000 IZ=1,lz1
|
|---|
| 198 | DO 1000 IY=1,ly1
|
|---|
| 199 | DO 1000 IX=1,lx1
|
|---|
| 200 | IF (ABS(TB(IX,IY,IZ,IE)).GT.EPS ) THEN
|
|---|
| 201 | WRITE(6,1005) IX,IY,IZ,IEG
|
|---|
| 202 | $ ,XM1(IX,IY,IZ,IE),YM1(IX,IY,IZ,IE),ZM1(IX,IY,IZ,IE)
|
|---|
| 203 | $ ,TA(IX,IY,IZ,IE),eps
|
|---|
| 204 | c WRITE(7,1005) IX,IY,IZ,IEG
|
|---|
| 205 | c $ ,XM1(IX,IY,IZ,IE),TB(IX,IY,IZ,IE),TA(IX,IY,IZ,IE)
|
|---|
| 206 | c $ ,QMASK(IX,IY,IZ,IE)
|
|---|
| 207 | 1005 FORMAT(2X,'WARNING: DSSUM problem at:',/
|
|---|
| 208 | $ ,1X,'I,J,K,IE:',3I5,i12,/
|
|---|
| 209 | $ ,2X,'Near X =',3G16.8,', d:',2G16.8)
|
|---|
| 210 | IERR=4
|
|---|
| 211 | ENDIF
|
|---|
| 212 | 1000 CONTINUE
|
|---|
| 213 | C
|
|---|
| 214 | C Set up QMASK quickly to annihilate checks on periodic bc's
|
|---|
| 215 | C
|
|---|
| 216 | CALL RONE(QMASK,NTOT)
|
|---|
| 217 | DO 100 IEL=1,NELT
|
|---|
| 218 | DO 100 IFACE=1,NFACES
|
|---|
| 219 | CB =CBC(IFACE,IEL,IFIELD)
|
|---|
| 220 | IF (CB.EQ.'P '.or.cb.eq.'p ')
|
|---|
| 221 | $ CALL FACEV(QMASK,IEL,IFACE,0.0,lx1,ly1,lz1)
|
|---|
| 222 | 100 CONTINUE
|
|---|
| 223 | CALL DSOP(QMASK,'MUL',lx1,ly1,lz1)
|
|---|
| 224 |
|
|---|
| 225 | c xxmin = glmin(xm1,ntot)
|
|---|
| 226 | c yymin = glmin(ym1,ntot)
|
|---|
| 227 | c zzmin = glmin(zm1,ntot)
|
|---|
| 228 | c xxmax = glmax(xm1,ntot)
|
|---|
| 229 | c yymax = glmax(ym1,ntot)
|
|---|
| 230 | c zzmax = glmax(zm1,ntot)
|
|---|
| 231 | c if (nio.eq.0) write(6,7) xxmin,yymin,zzmin,xxmax,yymax,zzmax
|
|---|
| 232 | c 7 format('xyz minmx2:',6g13.5)
|
|---|
| 233 |
|
|---|
| 234 |
|
|---|
| 235 |
|
|---|
| 236 | C
|
|---|
| 237 | C X-component
|
|---|
| 238 | C
|
|---|
| 239 | CALL COPY(TA,XM1,NTOT)
|
|---|
| 240 | CALL COPY(TB,XM1,NTOT)
|
|---|
| 241 | CALL DSOP(TA,'MIN',lx1,ly1,lz1)
|
|---|
| 242 | CALL DSOP(TB,'MAX',lx1,ly1,lz1)
|
|---|
| 243 | CALL SUB2(TA,XM1,NTOT)
|
|---|
| 244 | CALL SUB2(TB,XM1,NTOT)
|
|---|
| 245 | CALL COL2(TA,QMASK,NTOT)
|
|---|
| 246 | CALL COL2(TB,QMASK,NTOT)
|
|---|
| 247 | DO 1100 IE=1,NELT
|
|---|
| 248 | XSCMAX = VLMAX(XM1(1,1,1,IE),NXYZ1)
|
|---|
| 249 | XSCMIN = VLMIN(XM1(1,1,1,IE),NXYZ1)
|
|---|
| 250 | SCAL1=ABS(XSCMAX-XSCMIN)
|
|---|
| 251 | SCAL2=ABS(XSCMAX)
|
|---|
| 252 | SCAL3=ABS(XSCMIN)
|
|---|
| 253 | SCAL1=MAX(SCAL1,SCAL2)
|
|---|
| 254 | SCAL1=MAX(SCAL1,SCAL3)
|
|---|
| 255 | XSCALE = 1./SCAL1
|
|---|
| 256 | ieg=lglel(ie)
|
|---|
| 257 | DO 1100 IZ=1,lz1
|
|---|
| 258 | DO 1100 IY=1,ly1
|
|---|
| 259 | DO 1100 IX=1,lx1
|
|---|
| 260 | if (abs(ta(ix,iy,iz,ie)*xscale).gt.eps .or.
|
|---|
| 261 | $ abs(tb(ix,iy,iz,ie)*xscale).gt.eps ) then
|
|---|
| 262 | write(6,1105) ix,iy,iz,ieg
|
|---|
| 263 | $ ,xm1(ix,iy,iz,ie),ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
|
|---|
| 264 | $ ,tb(ix,iy,iz,ie),ta(ix,iy,iz,ie),XSCALE
|
|---|
| 265 | 1105 format(1x,'WARNING1 Element mesh mismatch at:',/
|
|---|
| 266 | $ ,1x,'i,j,k,ie:',3i5,I12,/
|
|---|
| 267 | $ ,1X,'Near X =',3G16.8,', d:',3G16.8)
|
|---|
| 268 | ierr=1
|
|---|
| 269 | endif
|
|---|
| 270 | 1100 CONTINUE
|
|---|
| 271 | C
|
|---|
| 272 | C Y-component
|
|---|
| 273 | C
|
|---|
| 274 | CALL COPY(TA,YM1,NTOT)
|
|---|
| 275 | CALL COPY(TB,YM1,NTOT)
|
|---|
| 276 | CALL DSOP(TA,'MIN',lx1,ly1,lz1)
|
|---|
| 277 | CALL DSOP(TB,'MAX',lx1,ly1,lz1)
|
|---|
| 278 | CALL SUB2(TA,YM1,NTOT)
|
|---|
| 279 | CALL SUB2(TB,YM1,NTOT)
|
|---|
| 280 | CALL COL2(TA,QMASK,NTOT)
|
|---|
| 281 | CALL COL2(TB,QMASK,NTOT)
|
|---|
| 282 | DO 1200 IE=1,NELT
|
|---|
| 283 | YSCMAX = VLMAX(YM1(1,1,1,IE),NXYZ1)
|
|---|
| 284 | YSCMIN = VLMIN(YM1(1,1,1,IE),NXYZ1)
|
|---|
| 285 | SCAL1=ABS(YSCMAX-YSCMIN)
|
|---|
| 286 | SCAL2=ABS(YSCMAX)
|
|---|
| 287 | SCAL3=ABS(YSCMIN)
|
|---|
| 288 | SCAL1=MAX(SCAL1,SCAL2)
|
|---|
| 289 | SCAL1=MAX(SCAL1,SCAL3)
|
|---|
| 290 | YSCALE = 1./SCAL1
|
|---|
| 291 | ieg=lglel(ie)
|
|---|
| 292 | DO 1200 IZ=1,lz1
|
|---|
| 293 | DO 1200 IY=1,ly1
|
|---|
| 294 | DO 1200 IX=1,lx1
|
|---|
| 295 | IF (ABS(TA(IX,IY,IZ,IE)*YSCALE).GT.EPS .OR.
|
|---|
| 296 | $ ABS(TB(IX,IY,IZ,IE)*YSCALE).GT.EPS ) THEN
|
|---|
| 297 | WRITE(6,1205) IX,IY,IZ,IEG
|
|---|
| 298 | $ ,XM1(IX,IY,IZ,IE),YM1(IX,IY,IZ,IE),ZM1(IX,IY,IZ,IE)
|
|---|
| 299 | $ ,TB(IX,IY,IZ,IE),TA(IX,IY,IZ,IE),yscale
|
|---|
| 300 | 1205 FORMAT(1X,'WARNING2 Element mesh mismatch at:',/
|
|---|
| 301 | $ ,1X,'I,J,K,IE:',3I5,i12,/
|
|---|
| 302 | $ ,1X,'Near Y =',3G16.8,', d:',3G16.8)
|
|---|
| 303 | IERR=2
|
|---|
| 304 | ENDIF
|
|---|
| 305 | 1200 CONTINUE
|
|---|
| 306 | C
|
|---|
| 307 | C Z-component
|
|---|
| 308 | C
|
|---|
| 309 | IF (IF3D) THEN
|
|---|
| 310 | CALL COPY(TA,ZM1,NTOT)
|
|---|
| 311 | CALL COPY(TB,ZM1,NTOT)
|
|---|
| 312 | CALL DSOP(TA,'MIN',lx1,ly1,lz1)
|
|---|
| 313 | CALL DSOP(TB,'MAX',lx1,ly1,lz1)
|
|---|
| 314 | CALL SUB2(TA,ZM1,NTOT)
|
|---|
| 315 | CALL SUB2(TB,ZM1,NTOT)
|
|---|
| 316 | CALL COL2(TA,QMASK,NTOT)
|
|---|
| 317 | CALL COL2(TB,QMASK,NTOT)
|
|---|
| 318 | DO 1300 IE=1,NELT
|
|---|
| 319 | ZSCMAX = VLMAX(ZM1(1,1,1,IE),NXYZ1)
|
|---|
| 320 | ZSCMIN = VLMIN(ZM1(1,1,1,IE),NXYZ1)
|
|---|
| 321 | SCAL1=ABS(ZSCMAX-ZSCMIN)
|
|---|
| 322 | SCAL2=ABS(ZSCMAX)
|
|---|
| 323 | SCAL3=ABS(ZSCMIN)
|
|---|
| 324 | SCAL1=MAX(SCAL1,SCAL2)
|
|---|
| 325 | SCAL1=MAX(SCAL1,SCAL3)
|
|---|
| 326 | ZSCALE = 1./SCAL1
|
|---|
| 327 | ieg=lglel(ie)
|
|---|
| 328 | DO 1300 IZ=1,lz1
|
|---|
| 329 | DO 1300 IY=1,ly1
|
|---|
| 330 | DO 1300 IX=1,lx1
|
|---|
| 331 | IF (ABS(TA(IX,IY,IZ,IE)*ZSCALE).GT.EPS .OR.
|
|---|
| 332 | $ ABS(TB(IX,IY,IZ,IE)*ZSCALE).GT.EPS ) THEN
|
|---|
| 333 | WRITE(6,1305) IX,IY,IZ,IEG
|
|---|
| 334 | $ ,XM1(IX,IY,IZ,IE),YM1(IX,IY,IZ,IE),ZM1(IX,IY,IZ,IE)
|
|---|
| 335 | $ ,TB(IX,IY,IZ,IE),TA(IX,IY,IZ,IE),zscale
|
|---|
| 336 | 1305 FORMAT(1X,'WARNING3 Element mesh mismatch at:',/
|
|---|
| 337 | $ ,1X,'I,J,K,IE:',3I5,i12,/
|
|---|
| 338 | $ ,1X,'Near Z =',3G16.8,', d:',3G16.8)
|
|---|
| 339 | IERR=3
|
|---|
| 340 | ENDIF
|
|---|
| 341 | 1300 CONTINUE
|
|---|
| 342 | ENDIF
|
|---|
| 343 |
|
|---|
| 344 | ierr = iglsum(ierr,1)
|
|---|
| 345 | IF (IERR.gt.0) THEN
|
|---|
| 346 | if(nid.eq.0) WRITE(6,1400)
|
|---|
| 347 | 1400 FORMAT
|
|---|
| 348 | $ (' Mesh consistency check failed. EXITING in VRDSMSH.')
|
|---|
| 349 | call exitt
|
|---|
| 350 | ENDIF
|
|---|
| 351 |
|
|---|
| 352 | tmp(1)=ierr
|
|---|
| 353 | CALL GOP(tmp,tmp(2),'M ',1)
|
|---|
| 354 | IF (tmp(1).ge.4.0) THEN
|
|---|
| 355 | WRITE(6,1400)
|
|---|
| 356 | $ (' Mesh consistency check failed. EXITING in VRDSMSH.')
|
|---|
| 357 | call exitt
|
|---|
| 358 | ENDIF
|
|---|
| 359 |
|
|---|
| 360 | if(nio.eq.0) then
|
|---|
| 361 | write(6,*) 'done :: verify mesh topology'
|
|---|
| 362 | write(6,*) ' '
|
|---|
| 363 | endif
|
|---|
| 364 |
|
|---|
| 365 | return
|
|---|
| 366 | end
|
|---|
| 367 | c-----------------------------------------------------------------------
|
|---|
| 368 | subroutine vrdsmshx ! verify mesh topology
|
|---|
| 369 | C
|
|---|
| 370 | C=====================================================================
|
|---|
| 371 | C Verify that mesh and dssum are properly defined by performing
|
|---|
| 372 | C a direct stiffness operation on the X,Y and Z coordinates.
|
|---|
| 373 | C Note that periodic faces are not checked here.
|
|---|
| 374 | C=====================================================================
|
|---|
| 375 | C
|
|---|
| 376 | INCLUDE 'SIZE'
|
|---|
| 377 | INCLUDE 'TOTAL'
|
|---|
| 378 | common /scrns/ tc(lx1,ly1,lz1,lelt),td(lx1,ly1,lz1,lelt)
|
|---|
| 379 | $ , ta(lx1,ly1,lz1,lelt),tb(lx1,ly1,lz1,lelt)
|
|---|
| 380 | $ , qmask(lx1,ly1,lz1,lelt)
|
|---|
| 381 | CHARACTER*3 CB
|
|---|
| 382 | C
|
|---|
| 383 | IERR = 0
|
|---|
| 384 | EPS = 1.0e-04
|
|---|
| 385 | EPS = 1.0e-03
|
|---|
| 386 | IFIELD = 1
|
|---|
| 387 | IF (IFHEAT) IFIELD = 2
|
|---|
| 388 | NXYZ1 = lx1*ly1*lz1
|
|---|
| 389 | NTOT = lx1*ly1*lz1*NELT
|
|---|
| 390 | NFACES = 2*ldim
|
|---|
| 391 |
|
|---|
| 392 | xmx = glmax(xm1,ntot)
|
|---|
| 393 | xmn = glmin(xm1,ntot)
|
|---|
| 394 | ymx = glmax(ym1,ntot)
|
|---|
| 395 | ymn = glmin(ym1,ntot)
|
|---|
| 396 | zmx = glmax(zm1,ntot)
|
|---|
| 397 | zmn = glmin(zm1,ntot)
|
|---|
| 398 | if (nio.eq.0) write(6,*) xmn,xmx,' Xrange'
|
|---|
| 399 | if (nio.eq.0) write(6,*) ymn,ymx,' Yrange'
|
|---|
| 400 | if (nio.eq.0) write(6,*) zmn,zmx,' Zrange'
|
|---|
| 401 | c return
|
|---|
| 402 | C
|
|---|
| 403 | C First check - use 1/Multiplicity
|
|---|
| 404 | C
|
|---|
| 405 | IF (IFHEAT) THEN
|
|---|
| 406 | CALL COPY(TA,TMULT,NTOT)
|
|---|
| 407 | ELSE
|
|---|
| 408 | CALL COPY(TA,VMULT,NTOT)
|
|---|
| 409 | ENDIF
|
|---|
| 410 | c
|
|---|
| 411 | c write(6,1)
|
|---|
| 412 | c $(nid,'tab4',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
|
|---|
| 413 | c 1 format(i3,a4,i3,16f5.2)
|
|---|
| 414 | c
|
|---|
| 415 | CALL DSSUM(TA,lx1,ly1,lz1)
|
|---|
| 416 | c
|
|---|
| 417 | c write(6,1)
|
|---|
| 418 | c $(nid,'taaf',lglel(ie),(ta(k,1,1,ie),k=1,lx1*ly1),ie=1,nelt)
|
|---|
| 419 | c
|
|---|
| 420 | CALL RONE (TB,NTOT)
|
|---|
| 421 | CALL SUB2 (TB,TA,NTOT)
|
|---|
| 422 | DO 1000 IE=1,NELT
|
|---|
| 423 | ieg=lglel(ie)
|
|---|
| 424 | DO 1000 IZ=1,lz1
|
|---|
| 425 | DO 1000 IY=1,ly1
|
|---|
| 426 | DO 1000 IX=1,lx1
|
|---|
| 427 | IF (ABS(TB(IX,IY,IZ,IE)).GT.EPS ) THEN
|
|---|
| 428 | WRITE(6,1005) IX,IY,IZ,IEG
|
|---|
| 429 | $ ,XM1(IX,IY,IZ,IE),YM1(IX,IY,IZ,IE),ZM1(IX,IY,IZ,IE)
|
|---|
| 430 | $ ,TA(IX,IY,IZ,IE),eps
|
|---|
| 431 | c WRITE(7,1005) IX,IY,IZ,IEG
|
|---|
| 432 | c $ ,XM1(IX,IY,IZ,IE),TB(IX,IY,IZ,IE),TA(IX,IY,IZ,IE)
|
|---|
| 433 | c $ ,QMASK(IX,IY,IZ,IE)
|
|---|
| 434 | 1005 FORMAT(2X,'WARNING: DSSUM problem at:',/
|
|---|
| 435 | $ ,1X,'I,J,K,IE:',3I5,i12,/
|
|---|
| 436 | $ ,2X,'Near X =',3G16.8,', d:',2G16.8)
|
|---|
| 437 | IERR=4
|
|---|
| 438 | ENDIF
|
|---|
| 439 | 1000 CONTINUE
|
|---|
| 440 | C
|
|---|
| 441 | C Set up QMASK quickly to annihilate checks on periodic bc's
|
|---|
| 442 | C
|
|---|
| 443 | CALL RONE(QMASK,NTOT)
|
|---|
| 444 | DO 100 IEL=1,NELT
|
|---|
| 445 | DO 100 IFACE=1,NFACES
|
|---|
| 446 | CB =CBC(IFACE,IEL,IFIELD)
|
|---|
| 447 | IF (CB.EQ.'P '.or.cb.eq.'p ')
|
|---|
| 448 | $ CALL FACEV(QMASK,IEL,IFACE,0.0,lx1,ly1,lz1)
|
|---|
| 449 | 100 CONTINUE
|
|---|
| 450 | call dsop(QMASK,'MUL',lx1,ly1,lz1)
|
|---|
| 451 |
|
|---|
| 452 | xxmin = glmin(xm1,ntot)
|
|---|
| 453 | yymin = glmin(ym1,ntot)
|
|---|
| 454 | zzmin = glmin(zm1,ntot)
|
|---|
| 455 | xxmax = glmax(xm1,ntot)
|
|---|
| 456 | yymax = glmax(ym1,ntot)
|
|---|
| 457 | zzmax = glmax(zm1,ntot)
|
|---|
| 458 | if (nio.eq.0) write(6,7) xxmin,yymin,zzmin,xxmax,yymax,zzmax
|
|---|
| 459 | 7 format('xyz minmx2:',6g13.5)
|
|---|
| 460 |
|
|---|
| 461 |
|
|---|
| 462 | C
|
|---|
| 463 | C X-component
|
|---|
| 464 | C
|
|---|
| 465 | call copy(ta,xm1,ntot)
|
|---|
| 466 | call copy(tb,xm1,ntot)
|
|---|
| 467 | call dsop(ta,'min',lx1,ly1,lz1)
|
|---|
| 468 | call dsop(tb,'max',lx1,ly1,lz1)
|
|---|
| 469 |
|
|---|
| 470 | call copy(tc,xm1,ntot)
|
|---|
| 471 | call copy(td,xm1,ntot)
|
|---|
| 472 | call dsop(tc,'min',lx1,ly1,lz1)
|
|---|
| 473 | call dsop(td,'max',lx1,ly1,lz1)
|
|---|
| 474 |
|
|---|
| 475 | xxmin = glmin(xm1,ntot)
|
|---|
| 476 | xxmax = glmax(xm1,ntot)
|
|---|
| 477 | yymax = glmax(ta ,ntot)
|
|---|
| 478 | yymin = glmin(ta ,ntot)
|
|---|
| 479 | zzmin = glmin(tb ,ntot)
|
|---|
| 480 | zzmax = glmax(tb ,ntot)
|
|---|
| 481 | if (nio.eq.0) write(6,9) xxmin,yymin,zzmin,xxmax,yymax,zzmax
|
|---|
| 482 | 9 format('xyz minmx3:',6g13.5)
|
|---|
| 483 |
|
|---|
| 484 | CALL SUB2(TA,XM1,NTOT)
|
|---|
| 485 | CALL SUB2(TB,XM1,NTOT)
|
|---|
| 486 |
|
|---|
| 487 | xxmin = glmin(qmask,ntot)
|
|---|
| 488 | xxmax = glmax(qmask,ntot)
|
|---|
| 489 | yymax = glmax(ta ,ntot)
|
|---|
| 490 | yymin = glmin(ta ,ntot)
|
|---|
| 491 | zzmin = glmin(tb ,ntot)
|
|---|
| 492 | zzmax = glmax(tb ,ntot)
|
|---|
| 493 | if (nio.eq.0) write(6,19) xxmin,yymin,zzmin,xxmax,yymax,zzmax
|
|---|
| 494 | 19 format('xyz minmx4:',6g13.5)
|
|---|
| 495 |
|
|---|
| 496 | CALL COL2(TA,QMASK,NTOT)
|
|---|
| 497 | CALL COL2(TB,QMASK,NTOT)
|
|---|
| 498 |
|
|---|
| 499 | xxmin = glmin(qmask,ntot)
|
|---|
| 500 | xxmax = glmax(qmask,ntot)
|
|---|
| 501 | yymax = glmax(ta ,ntot)
|
|---|
| 502 | yymin = glmin(ta ,ntot)
|
|---|
| 503 | zzmin = glmin(tb ,ntot)
|
|---|
| 504 | zzmax = glmax(tb ,ntot)
|
|---|
| 505 | if (nio.eq.0) write(6,29) xxmin,yymin,zzmin,xxmax,yymax,zzmax
|
|---|
| 506 | 29 format('xyz minmx5:',6g13.5)
|
|---|
| 507 |
|
|---|
| 508 | DO 1100 IE=1,NELT
|
|---|
| 509 | XSCMAX = VLMAX(XM1(1,1,1,IE),NXYZ1)
|
|---|
| 510 | XSCMIN = VLMIN(XM1(1,1,1,IE),NXYZ1)
|
|---|
| 511 | SCAL1=ABS(XSCMAX-XSCMIN)
|
|---|
| 512 | SCAL2=ABS(XSCMAX)
|
|---|
| 513 | SCAL3=ABS(XSCMIN)
|
|---|
| 514 | SCAL1=MAX(SCAL1,SCAL2)
|
|---|
| 515 | SCAL1=MAX(SCAL1,SCAL3)
|
|---|
| 516 | XSCALE = 1./SCAL1
|
|---|
| 517 | ieg=lglel(ie)
|
|---|
| 518 | DO 1100 IZ=1,lz1
|
|---|
| 519 | DO 1100 IY=1,ly1
|
|---|
| 520 | DO 1100 IX=1,lx1
|
|---|
| 521 | if (abs(ta(ix,iy,iz,ie)*xscale).gt.eps .or.
|
|---|
| 522 | $ abs(tb(ix,iy,iz,ie)*xscale).gt.eps ) then
|
|---|
| 523 | write(6,1105) nid,ix,iy,iz,ie,ieg
|
|---|
| 524 | $ ,xm1(ix,iy,iz,ie),tc(ix,iy,iz,ie),td(ix,iy,iz,ie)
|
|---|
| 525 | $ ,ym1(ix,iy,iz,ie),zm1(ix,iy,iz,ie)
|
|---|
| 526 | $ ,ta(ix,iy,iz,ie),tb(ix,iy,iz,ie),xscale
|
|---|
| 527 | $ ,qmask(ix,iy,iz,ie)
|
|---|
| 528 | 1105 format(i4.4,1x,'ie:',3i3,i10,i10,1p9e11.3)
|
|---|
| 529 | c1105 format(i4.4,1x,'ie:',3i3,i6,1p9e11.3)
|
|---|
| 530 | ierr=1
|
|---|
| 531 | goto 1101
|
|---|
| 532 | endif
|
|---|
| 533 | 1100 CONTINUE
|
|---|
| 534 | 1101 CONTINUE
|
|---|
| 535 |
|
|---|
| 536 | xxmin = glmin(xm1,ntot)
|
|---|
| 537 | xxmax = glmax(xm1,ntot)
|
|---|
| 538 | yymax = glmax(ta ,ntot)
|
|---|
| 539 | yymin = glmin(ta ,ntot)
|
|---|
| 540 | zzmin = glmin(tb ,ntot)
|
|---|
| 541 | zzmax = glmax(tb ,ntot)
|
|---|
| 542 | if (nio.eq.0) write(6,39) xxmin,yymin,zzmin,xxmax,yymax,zzmax
|
|---|
| 543 | 39 format('xyz minmx5:',6g13.5)
|
|---|
| 544 |
|
|---|
| 545 | c ifvo = .true.
|
|---|
| 546 | c ifpo = .false.
|
|---|
| 547 | c ifto = .true.
|
|---|
| 548 | c call outpost(xm1,ta,tb,pr,qmask,' ')
|
|---|
| 549 | c call exitt
|
|---|
| 550 |
|
|---|
| 551 | return
|
|---|
| 552 | end
|
|---|
| 553 | c-----------------------------------------------------------------------
|
|---|
| 554 | subroutine rotat2(xyz,angle,npts)
|
|---|
| 555 | C
|
|---|
| 556 | C Rotate NPTS through ANGLE (in two directions IF3D).
|
|---|
| 557 | C
|
|---|
| 558 | INCLUDE 'SIZE'
|
|---|
| 559 | INCLUDE 'INPUT'
|
|---|
| 560 | DIMENSION XYZ(3,1)
|
|---|
| 561 | COMMON /CTMP0/ RMTRX(3,3),RX(3,3),RZ(3,3),XYZN(3,10)
|
|---|
| 562 | C
|
|---|
| 563 | SINA=SIN(ANGLE)
|
|---|
| 564 | COSA=COS(ANGLE)
|
|---|
| 565 | CALL RZERO(RX,9)
|
|---|
| 566 | CALL RZERO(RZ,9)
|
|---|
| 567 | RX(1,1)=COSA
|
|---|
| 568 | RX(2,2)=COSA
|
|---|
| 569 | RX(1,2)=SINA
|
|---|
| 570 | RX(2,1)=-SINA
|
|---|
| 571 | RX(3,3)=1.0
|
|---|
| 572 | IF (IF3D) THEN
|
|---|
| 573 | RZ(1,1)=COSA
|
|---|
| 574 | RZ(3,3)=COSA
|
|---|
| 575 | RZ(1,3)=SINA
|
|---|
| 576 | RZ(3,1)=-SINA
|
|---|
| 577 | RZ(2,2)=1.0
|
|---|
| 578 | ELSE
|
|---|
| 579 | RZ(1,1)=1.0
|
|---|
| 580 | RZ(2,2)=1.0
|
|---|
| 581 | RZ(3,3)=1.0
|
|---|
| 582 | ENDIF
|
|---|
| 583 | CALL MXM(RX,3,RZ,3,RMTRX,3)
|
|---|
| 584 | C
|
|---|
| 585 | C Strip mine mxms in chunks of 10:
|
|---|
| 586 | DO 100 I=1,NPTS-10,10
|
|---|
| 587 | CALL MXM(RMTRX,3,XYZ(1,I),3,XYZN,10)
|
|---|
| 588 | CALL COPY(XYZ(1,I),XYZN,30)
|
|---|
| 589 | 100 CONTINUE
|
|---|
| 590 | N10=MOD1(NPTS,10)
|
|---|
| 591 | I=NPTS-N10+1
|
|---|
| 592 | CALL RZERO(XYZN,30)
|
|---|
| 593 | IF (N10.GT.0) THEN
|
|---|
| 594 | CALL MXM(RMTRX,3,XYZ(1,I),3,XYZN,N10)
|
|---|
| 595 | CALL COPY(XYZ(1,I),XYZN,3*N10)
|
|---|
| 596 | ENDIF
|
|---|
| 597 | C
|
|---|
| 598 | return
|
|---|
| 599 | END
|
|---|
| 600 | c-----------------------------------------------------------------------
|
|---|
| 601 | subroutine scale(xyzl,nl)
|
|---|
| 602 | C
|
|---|
| 603 | C Rescale XYZL such that the mean value of IXX=IYY=IZZ for each element.
|
|---|
| 604 | C
|
|---|
| 605 | INCLUDE 'SIZE'
|
|---|
| 606 | INCLUDE 'INPUT'
|
|---|
| 607 | DIMENSION XYZL(3,8,LELT)
|
|---|
| 608 | COMMON /CTMP0/ VO(LELT),XYZI(3,LELT),CG(3,LELT)
|
|---|
| 609 | $ ,TI(6),WORK(6)
|
|---|
| 610 | C
|
|---|
| 611 | C Compute volumes -
|
|---|
| 612 | C
|
|---|
| 613 | CALL VOLUME2(VO,XYZL,NL)
|
|---|
| 614 | VTOT=GLSUM (VO,NL)
|
|---|
| 615 | C
|
|---|
| 616 | C Compute (weighted) average inertia for each element.
|
|---|
| 617 | C
|
|---|
| 618 | NCRNR=2**ldim
|
|---|
| 619 | CALL RZERO(TI,6)
|
|---|
| 620 | DO 100 IL=1,NL
|
|---|
| 621 | VO0 = VO(IL)/VTOT
|
|---|
| 622 | CALL INRTIA(XYZI(1,IL),CG(1,IL),XYZL(1,1,IL),NCRNR,1)
|
|---|
| 623 | TI(1)=TI(1)+XYZI(1,IL)*VO0
|
|---|
| 624 | TI(2)=TI(2)+XYZI(2,IL)*VO0
|
|---|
| 625 | TI(3)=TI(3)+XYZI(3,IL)*VO0
|
|---|
| 626 | TI(4)=TI(4)+CG(1,IL) *VO0
|
|---|
| 627 | TI(5)=TI(5)+CG(2,IL) *VO0
|
|---|
| 628 | TI(6)=TI(6)+CG(3,IL) *VO0
|
|---|
| 629 | 100 CONTINUE
|
|---|
| 630 | CALL GOP(TI,WORK,'+ ',6)
|
|---|
| 631 | XI =SQRT(TI(1))
|
|---|
| 632 | YI =SQRT(TI(2))
|
|---|
| 633 | ZI =1.0
|
|---|
| 634 | IF (IF3D) ZI=SQRT(TI(3))
|
|---|
| 635 | C
|
|---|
| 636 | C Rescale ( & shift to a nearly mean zero )
|
|---|
| 637 | C
|
|---|
| 638 | DO 200 IL=1,NL
|
|---|
| 639 | DO 200 IC=1,NCRNR
|
|---|
| 640 | XYZL(1,IC,IL)=(XYZL(1,IC,IL)-TI(4))/XI
|
|---|
| 641 | XYZL(2,IC,IL)=(XYZL(2,IC,IL)-TI(5))/YI
|
|---|
| 642 | XYZL(3,IC,IL)=(XYZL(3,IC,IL)-TI(6))/ZI
|
|---|
| 643 | 200 CONTINUE
|
|---|
| 644 | C
|
|---|
| 645 | return
|
|---|
| 646 | END
|
|---|
| 647 | c-----------------------------------------------------------------------
|
|---|
| 648 | subroutine inrtia(xyzi,cg,xyzl,n,itype)
|
|---|
| 649 | C
|
|---|
| 650 | C Compute cg and inertia for a collection of unit point masses.
|
|---|
| 651 | C This is a global (multiprocessor) operation, only IF itype=2.
|
|---|
| 652 | C
|
|---|
| 653 | DIMENSION XYZI(3),CG(3),XYZL(3,1)
|
|---|
| 654 | DIMENSION TI(4),WORK(4)
|
|---|
| 655 | C
|
|---|
| 656 | TI(1)=0.0
|
|---|
| 657 | TI(2)=0.0
|
|---|
| 658 | TI(3)=0.0
|
|---|
| 659 | TI(4)=N
|
|---|
| 660 | DO 100 I=1,N
|
|---|
| 661 | TI(1)=TI(1)+XYZL(1,I)
|
|---|
| 662 | TI(2)=TI(2)+XYZL(2,I)
|
|---|
| 663 | TI(3)=TI(3)+XYZL(3,I)
|
|---|
| 664 | 100 CONTINUE
|
|---|
| 665 | IF (ITYPE.EQ.2) CALL GOP(TI,WORK,'+ ',4)
|
|---|
| 666 | IF (TI(4).EQ.0.0) TI(4)=1.0
|
|---|
| 667 | CG(1)=TI(1)/TI(4)
|
|---|
| 668 | CG(2)=TI(2)/TI(4)
|
|---|
| 669 | CG(3)=TI(3)/TI(4)
|
|---|
| 670 | C
|
|---|
| 671 | TI(1)=0.0
|
|---|
| 672 | TI(2)=0.0
|
|---|
| 673 | TI(3)=0.0
|
|---|
| 674 | DO 200 I=1,N
|
|---|
| 675 | TI(1)=TI(1)+( XYZL(1,I)-CG(1) )**2
|
|---|
| 676 | TI(2)=TI(2)+( XYZL(2,I)-CG(2) )**2
|
|---|
| 677 | TI(3)=TI(3)+( XYZL(3,I)-CG(3) )**2
|
|---|
| 678 | 200 CONTINUE
|
|---|
| 679 | IF (ITYPE.EQ.2) CALL GOP(TI,WORK,'+ ',3)
|
|---|
| 680 | TI(1)=TI(1)/TI(4)
|
|---|
| 681 | TI(2)=TI(2)/TI(4)
|
|---|
| 682 | TI(3)=TI(3)/TI(4)
|
|---|
| 683 | IF (ITYPE.EQ.2) THEN
|
|---|
| 684 | C std. def'n of inertia.
|
|---|
| 685 | XYZI(1)=TI(2)+TI(3)
|
|---|
| 686 | XYZI(2)=TI(3)+TI(1)
|
|---|
| 687 | XYZI(3)=TI(1)+TI(2)
|
|---|
| 688 | ELSE
|
|---|
| 689 | XYZI(1)=TI(1)
|
|---|
| 690 | XYZI(2)=TI(2)
|
|---|
| 691 | XYZI(3)=TI(3)
|
|---|
| 692 | ENDIF
|
|---|
| 693 | C
|
|---|
| 694 | return
|
|---|
| 695 | END
|
|---|
| 696 | c-----------------------------------------------------------------------
|
|---|
| 697 | subroutine volume2(vol,xyz,n)
|
|---|
| 698 | INCLUDE 'SIZE'
|
|---|
| 699 | INCLUDE 'INPUT'
|
|---|
| 700 | DIMENSION XYZ(3,2,2,2,1)
|
|---|
| 701 | DIMENSION VOL(1)
|
|---|
| 702 | C
|
|---|
| 703 | DO 1000 IE=1,N
|
|---|
| 704 | VOL(IE)=0.0
|
|---|
| 705 | IF (IF3D) THEN
|
|---|
| 706 | DO 20 K=1,2
|
|---|
| 707 | DO 20 J=1,2
|
|---|
| 708 | DO 20 I=1,2
|
|---|
| 709 | VOL1 = (XYZ(1,2,J,K,IE)-XYZ(1,1,J,K,IE))
|
|---|
| 710 | $ * (XYZ(2,I,2,K,IE)-XYZ(2,I,1,K,IE))
|
|---|
| 711 | $ * (XYZ(3,I,J,2,IE)-XYZ(3,I,J,1,IE))
|
|---|
| 712 | VOL2 = (XYZ(1,2,J,K,IE)-XYZ(1,1,J,K,IE))
|
|---|
| 713 | $ * (XYZ(2,I,J,2,IE)-XYZ(2,I,J,1,IE))
|
|---|
| 714 | $ * (XYZ(3,I,2,K,IE)-XYZ(3,I,1,K,IE))
|
|---|
| 715 | VOL3 = (XYZ(1,I,2,K,IE)-XYZ(1,I,1,K,IE))
|
|---|
| 716 | $ * (XYZ(2,2,J,K,IE)-XYZ(2,1,J,K,IE))
|
|---|
| 717 | $ * (XYZ(3,I,J,2,IE)-XYZ(3,I,J,1,IE))
|
|---|
| 718 | VOL4 = (XYZ(1,I,J,2,IE)-XYZ(1,I,J,1,IE))
|
|---|
| 719 | $ * (XYZ(2,I,2,K,IE)-XYZ(2,I,1,K,IE))
|
|---|
| 720 | $ * (XYZ(3,I,2,K,IE)-XYZ(3,I,1,K,IE))
|
|---|
| 721 | VOL5 = (XYZ(1,I,2,K,IE)-XYZ(1,I,1,K,IE))
|
|---|
| 722 | $ * (XYZ(2,I,J,2,IE)-XYZ(2,I,J,1,IE))
|
|---|
| 723 | $ * (XYZ(3,2,J,K,IE)-XYZ(3,1,J,K,IE))
|
|---|
| 724 | VOL6 = (XYZ(1,I,J,2,IE)-XYZ(1,I,J,1,IE))
|
|---|
| 725 | $ * (XYZ(2,I,2,K,IE)-XYZ(2,I,1,K,IE))
|
|---|
| 726 | $ * (XYZ(3,2,J,K,IE)-XYZ(3,1,J,K,IE))
|
|---|
| 727 | VOL(IE) = VOL(IE)+VOL1+VOL2+VOL3+VOL4+VOL5+VOL6
|
|---|
| 728 | 20 CONTINUE
|
|---|
| 729 | VOL(IE)=VOL(IE)/8.0
|
|---|
| 730 | ELSE
|
|---|
| 731 | C 2-D:
|
|---|
| 732 | DO 40 J=1,2
|
|---|
| 733 | DO 40 I=1,2
|
|---|
| 734 | VOL1 = (XYZ(1,2,J,1,IE)-XYZ(1,1,J,1,IE))
|
|---|
| 735 | $ * (XYZ(2,I,2,1,IE)-XYZ(2,I,1,1,IE))
|
|---|
| 736 | VOL3 = (XYZ(1,I,2,1,IE)-XYZ(1,I,1,1,IE))
|
|---|
| 737 | $ * (XYZ(2,2,J,1,IE)-XYZ(2,1,J,1,IE))
|
|---|
| 738 | VOL(IE)=VOL(IE)+VOL1+VOL3
|
|---|
| 739 | 40 CONTINUE
|
|---|
| 740 | VOL(IE)=VOL(IE)/4.0
|
|---|
| 741 | ENDIF
|
|---|
| 742 | VOL(IE)=ABS(VOL(IE))
|
|---|
| 743 | 1000 CONTINUE
|
|---|
| 744 | C
|
|---|
| 745 | return
|
|---|
| 746 | END
|
|---|
| 747 | c-----------------------------------------------------------------------
|
|---|
| 748 | subroutine findcg(cg,xyz,n)
|
|---|
| 749 | C
|
|---|
| 750 | C Compute cg for N elements.
|
|---|
| 751 | C
|
|---|
| 752 | INCLUDE 'SIZE'
|
|---|
| 753 | DIMENSION CG(3,1),XYZ(3,8,1)
|
|---|
| 754 | C
|
|---|
| 755 | NCRNR=2**ldim
|
|---|
| 756 | CALL RZERO(CG,3*N)
|
|---|
| 757 | DO 100 I =1,N
|
|---|
| 758 | DO 100 IC=1,NCRNR
|
|---|
| 759 | CG(1,I)=CG(1,I)+XYZ(1,IC,I)
|
|---|
| 760 | CG(2,I)=CG(2,I)+XYZ(2,IC,I)
|
|---|
| 761 | CG(3,I)=CG(3,I)+XYZ(3,IC,I)
|
|---|
| 762 | 100 CONTINUE
|
|---|
| 763 | TMP=1.0/(NCRNR)
|
|---|
| 764 | CALL CMULT(CG,TMP,3*N)
|
|---|
| 765 | return
|
|---|
| 766 | END
|
|---|
| 767 | c-----------------------------------------------------------------------
|
|---|
| 768 | subroutine divide(list1,list2,nl1,nl2,ifok,list,nl,xyzi,cg,WGT)
|
|---|
| 769 | C
|
|---|
| 770 | C Divide the elements associated with this subdomain according to
|
|---|
| 771 | C the direction having the smallest moment of inertia (the "long"
|
|---|
| 772 | C direction).
|
|---|
| 773 | C
|
|---|
| 774 | INCLUDE 'SIZE'
|
|---|
| 775 | INCLUDE 'INPUT'
|
|---|
| 776 | INCLUDE 'PARALLEL'
|
|---|
| 777 | INCLUDE 'TSTEP'
|
|---|
| 778 | C
|
|---|
| 779 | DIMENSION LIST(LELT),LIST1(LELT),LIST2(LELT)
|
|---|
| 780 | DIMENSION XYZI(3),CG(3,LELT),wgt(1)
|
|---|
| 781 | COMMON /CTMP0/ XCG(LELT),YCG(LELT),ZCG(LELT)
|
|---|
| 782 | REAL IXX,IYY,IZZ
|
|---|
| 783 | INTEGER WORK(2),WRK2(2)
|
|---|
| 784 | LOGICAL IFOK
|
|---|
| 785 | C
|
|---|
| 786 | C Choose "long" direction:
|
|---|
| 787 | C
|
|---|
| 788 | IXX=XYZI(1)
|
|---|
| 789 | IYY=XYZI(2)
|
|---|
| 790 | IZZ=XYZI(3)
|
|---|
| 791 | IF (IF3D) THEN
|
|---|
| 792 | IF (IXX.LE.IYY.AND.IXX.LE.IZZ) THEN
|
|---|
| 793 | DO 104 IE=1,NL
|
|---|
| 794 | XCG(IE)=CG(1,IE)
|
|---|
| 795 | YCG(IE)=CG(2,IE)
|
|---|
| 796 | ZCG(IE)=CG(3,IE)
|
|---|
| 797 | 104 CONTINUE
|
|---|
| 798 | ELSEIF (IYY.LE.IXX.AND.IYY.LE.IZZ) THEN
|
|---|
| 799 | DO 106 IE=1,NL
|
|---|
| 800 | XCG(IE)=CG(2,IE)
|
|---|
| 801 | YCG(IE)=CG(3,IE)
|
|---|
| 802 | ZCG(IE)=CG(1,IE)
|
|---|
| 803 | 106 CONTINUE
|
|---|
| 804 | ELSEIF (IZZ.LE.IXX.AND.IZZ.LE.IYY) THEN
|
|---|
| 805 | DO 108 IE=1,NL
|
|---|
| 806 | XCG(IE)=CG(3,IE)
|
|---|
| 807 | YCG(IE)=CG(1,IE)
|
|---|
| 808 | ZCG(IE)=CG(2,IE)
|
|---|
| 809 | 108 CONTINUE
|
|---|
| 810 | ENDIF
|
|---|
| 811 | ELSE
|
|---|
| 812 | C 2-D:
|
|---|
| 813 | IF (IXX.LE.IYY) THEN
|
|---|
| 814 | DO 114 IE=1,NL
|
|---|
| 815 | XCG(IE)=CG(1,IE)
|
|---|
| 816 | YCG(IE)=CG(2,IE)
|
|---|
| 817 | 114 CONTINUE
|
|---|
| 818 | ELSE
|
|---|
| 819 | DO 116 IE=1,NL
|
|---|
| 820 | XCG(IE)=CG(2,IE)
|
|---|
| 821 | YCG(IE)=CG(1,IE)
|
|---|
| 822 | 116 CONTINUE
|
|---|
| 823 | ENDIF
|
|---|
| 824 | ENDIF
|
|---|
| 825 | call col2(xcg,wgt,nl)
|
|---|
| 826 | call col2(ycg,wgt,nl)
|
|---|
| 827 | call col2(zcg,wgt,nl)
|
|---|
| 828 | C
|
|---|
| 829 | C Find median value of CG to determine dividing point:
|
|---|
| 830 | C
|
|---|
| 831 | XM=FMDIAN(XCG,NL,IFOK)
|
|---|
| 832 | YM=FMDIAN(YCG,NL,IFOK)
|
|---|
| 833 | ZM=0.0
|
|---|
| 834 | IF (IF3D) ZM=FMDIAN(ZCG,NL,IFOK)
|
|---|
| 835 | C
|
|---|
| 836 | C Diagnostics
|
|---|
| 837 | C
|
|---|
| 838 | IF (.NOT.IFOK) THEN
|
|---|
| 839 | WRITE(6,130) NID,NL,XM,YM,ZM
|
|---|
| 840 | DO 120 IL=1,NL
|
|---|
| 841 | WRITE(6,135) NID,IL,XCG(IL),YCG(IL),ZCG(IL)
|
|---|
| 842 | 120 CONTINUE
|
|---|
| 843 | 130 FORMAT(I10,'DIVIDE: NL,XM,YM,ZM',I3,3F12.5)
|
|---|
| 844 | 135 FORMAT(I10,'DIVIDE: NID,IL,XC,YC,ZCG',I4,3F12.5)
|
|---|
| 845 | ENDIF
|
|---|
| 846 | C
|
|---|
| 847 | C=============================================================
|
|---|
| 848 | C Divide LIST into LIST1 (XCG < XM) and LIST2 (XCG>XM).
|
|---|
| 849 | C=============================================================
|
|---|
| 850 | C
|
|---|
| 851 | NL1=0
|
|---|
| 852 | NL2=0
|
|---|
| 853 | DO 200 IE=1,NL
|
|---|
| 854 | IF (XCG(IE).LT.XM) THEN
|
|---|
| 855 | NL1=NL1+1
|
|---|
| 856 | LIST1(NL1)=LIST(IE)
|
|---|
| 857 | ENDIF
|
|---|
| 858 | IF (XCG(IE).GT.XM) THEN
|
|---|
| 859 | NL2=NL2+1
|
|---|
| 860 | LIST2(NL2)=LIST(IE)
|
|---|
| 861 | ENDIF
|
|---|
| 862 | IF (XCG(IE).EQ.XM) THEN
|
|---|
| 863 | C
|
|---|
| 864 | C We have to look at the other directions to arrive at
|
|---|
| 865 | C a unique subdivision algortithm.
|
|---|
| 866 | C
|
|---|
| 867 | C
|
|---|
| 868 | C More Diagnostics
|
|---|
| 869 | C
|
|---|
| 870 | IF (.NOT.IFOK) WRITE(6,201) NID,IE,XCG(IE),XM
|
|---|
| 871 | 201 FORMAT(I10,'DIVIDE: IE,XCG,XM:',I4,3F12.5)
|
|---|
| 872 | C
|
|---|
| 873 | IF (YCG(IE).LT.YM) THEN
|
|---|
| 874 | NL1=NL1+1
|
|---|
| 875 | LIST1(NL1)=LIST(IE)
|
|---|
| 876 | ENDIF
|
|---|
| 877 | IF (YCG(IE).GT.YM) THEN
|
|---|
| 878 | NL2=NL2+1
|
|---|
| 879 | LIST2(NL2)=LIST(IE)
|
|---|
| 880 | ENDIF
|
|---|
| 881 | IF (YCG(IE).EQ.YM) THEN
|
|---|
| 882 | C look at 3rd direction.
|
|---|
| 883 | IF (IF3D .AND. ZCG(IE).LT.ZM) THEN
|
|---|
| 884 | NL1=NL1+1
|
|---|
| 885 | LIST1(NL1)=LIST(IE)
|
|---|
| 886 | ELSE IF (IF3D .AND. ZCG(IE).GT.ZM) THEN
|
|---|
| 887 | NL2=NL2+1
|
|---|
| 888 | LIST2(NL2)=LIST(IE)
|
|---|
| 889 | ELSE
|
|---|
| 890 | C for 2- or 3-D intdeterminate case:
|
|---|
| 891 | NL1=NL1+1
|
|---|
| 892 | LIST1(NL1)=LIST(IE)
|
|---|
| 893 | ENDIF
|
|---|
| 894 | ENDIF
|
|---|
| 895 | C
|
|---|
| 896 | ENDIF
|
|---|
| 897 | 200 CONTINUE
|
|---|
| 898 | C
|
|---|
| 899 | C Check for an even distribution (i.e. - not different by
|
|---|
| 900 | C more than 1):
|
|---|
| 901 | C
|
|---|
| 902 | IFOK=.TRUE.
|
|---|
| 903 | WORK(1)=NL1
|
|---|
| 904 | WORK(2)=NL2
|
|---|
| 905 | CALL IGOP(WORK,WRK2,'+ ',2)
|
|---|
| 906 | IF (ABS(WORK(1)-WORK(2)).GT.1) IFOK=.FALSE.
|
|---|
| 907 | C
|
|---|
| 908 | return
|
|---|
| 909 | END
|
|---|
| 910 | c-----------------------------------------------------------------------
|
|---|
| 911 | subroutine bufchk(buf,n)
|
|---|
| 912 | integer n
|
|---|
| 913 | real buf(n)
|
|---|
| 914 | do i=1,n
|
|---|
| 915 | write(6,*) buf(i), ' whhhh'
|
|---|
| 916 | enddo
|
|---|
| 917 | return
|
|---|
| 918 | end
|
|---|
| 919 | c-----------------------------------------------------------------------
|
|---|
| 920 | subroutine chk_xyz
|
|---|
| 921 | include 'SIZE'
|
|---|
| 922 | include 'TOTAL'
|
|---|
| 923 | integer e,f,eg
|
|---|
| 924 |
|
|---|
| 925 | do e=1,nelt
|
|---|
| 926 | eg = lglel(e)
|
|---|
| 927 | write(6,1) nid,eg,e,(cbc(f,e,1),f=1,6)
|
|---|
| 928 | enddo
|
|---|
| 929 | 1 format(3i12,6(1x,a3),' cbc')
|
|---|
| 930 |
|
|---|
| 931 | return
|
|---|
| 932 | end
|
|---|
| 933 | c-----------------------------------------------------------------------
|
|---|
| 934 | subroutine chk_nel
|
|---|
| 935 | include 'SIZE'
|
|---|
| 936 | include 'TOTAL'
|
|---|
| 937 |
|
|---|
| 938 | neltmx=np*lelt
|
|---|
| 939 | nelvmx=np*lelv
|
|---|
| 940 |
|
|---|
| 941 | neltmx=min(neltmx,lelg)
|
|---|
| 942 | nelvmx=min(nelvmx,lelg)
|
|---|
| 943 |
|
|---|
| 944 | nelgt = iglmax(nelgt,1)
|
|---|
| 945 | nelgv = iglmax(nelgv,1)
|
|---|
| 946 |
|
|---|
| 947 | c write(6,*) nid,' inside chk_nel',nelgt,neltmx,nelvmx
|
|---|
| 948 |
|
|---|
| 949 | if (nelgt.gt.neltmx.or.nelgv.gt.nelvmx) then
|
|---|
| 950 | if (nid.eq.0) then
|
|---|
| 951 | lelt_needed = nelgt/np
|
|---|
| 952 | if (mod(nelgt,np).ne.0) lelt_needed = lelt_needed + 1
|
|---|
| 953 | write(6,12) lelt,lelg,lelt_needed,np,nelgt
|
|---|
| 954 | 12 format(//,2X,'ABORT: Problem size too large!'
|
|---|
| 955 | $ ,/,2X
|
|---|
| 956 | $ ,/,2X,'This solver has been compiled for:'
|
|---|
| 957 | $ ,/,2X,' number of elements/proc (lelt):',i12
|
|---|
| 958 | $ ,/,2X,' total number of elements (lelg):',i12
|
|---|
| 959 | $ ,/,2X
|
|---|
| 960 | $ ,/,2X,'Recompile with the following SIZE parameters:'
|
|---|
| 961 | $ ,/,2X,' lelt >= ',i12,' for np = ',i12
|
|---|
| 962 | $ ,/,2X,' lelg >= ',i12,/)
|
|---|
| 963 | c write(6,*)'help:',lp,np,nelvmx,nelgv,neltmx,nelgt
|
|---|
| 964 | c write(6,*)'help:',lelt,lelv,lelgv
|
|---|
| 965 | endif
|
|---|
| 966 | call exitt
|
|---|
| 967 | endif
|
|---|
| 968 |
|
|---|
| 969 | if(nelgt.gt.nelgt_max) then
|
|---|
| 970 | if(nid.eq.0) write(6,*)
|
|---|
| 971 | $ 'ABORT: Total number of elements too large!',
|
|---|
| 972 | $ ' nel_max = ', nelgt_max
|
|---|
| 973 | call exitt
|
|---|
| 974 | endif
|
|---|
| 975 |
|
|---|
| 976 | if (nelt.gt.lelt) then
|
|---|
| 977 | write(6,'(A,3I12)') 'ABORT: nelt>lelt!', nid, nelt, lelt
|
|---|
| 978 | call exitt
|
|---|
| 979 | endif
|
|---|
| 980 |
|
|---|
| 981 | return
|
|---|
| 982 | end
|
|---|
| 983 | c-----------------------------------------------------------------------
|
|---|
| 984 | subroutine cscan(sout,key,nk)
|
|---|
| 985 |
|
|---|
| 986 | character*132 sout,key
|
|---|
| 987 | character*132 string
|
|---|
| 988 | character*1 string1(132)
|
|---|
| 989 | equivalence (string1,string)
|
|---|
| 990 | c
|
|---|
| 991 | do i=1,100000000
|
|---|
| 992 | call blank(string,132)
|
|---|
| 993 | read (nk,80,end=100,err=100) string
|
|---|
| 994 | call chcopy(sout,string,132)
|
|---|
| 995 | c write (6,*) string
|
|---|
| 996 | if (indx1(string,key,nk).ne.0) return
|
|---|
| 997 | enddo
|
|---|
| 998 | 100 continue
|
|---|
| 999 | c
|
|---|
| 1000 | 80 format(a132)
|
|---|
| 1001 | return
|
|---|
| 1002 |
|
|---|
| 1003 | end
|
|---|
| 1004 | c-----------------------------------------------------------------------
|
|---|
