| 1 | c-----------------------------------------------------------------------
|
|---|
| 2 | c
|
|---|
| 3 | subroutine set_vert(glo_num,ngv,nx,nel,vertex,ifcenter)
|
|---|
| 4 | c
|
|---|
| 5 | c Given global array, vertex, pointing to hex vertices, set up
|
|---|
| 6 | c a new array of global pointers for an nx^ldim set of elements.
|
|---|
| 7 | c
|
|---|
| 8 | include 'SIZE'
|
|---|
| 9 | include 'INPUT'
|
|---|
| 10 | c
|
|---|
| 11 | integer*8 glo_num(1),ngv
|
|---|
| 12 | integer vertex(1),nx
|
|---|
| 13 | logical ifcenter
|
|---|
| 14 |
|
|---|
| 15 | if (if3d) then
|
|---|
| 16 | call setvert3d(glo_num,ngv,nx,nel,vertex,ifcenter)
|
|---|
| 17 | else
|
|---|
| 18 | call setvert2d(glo_num,ngv,nx,nel,vertex,ifcenter)
|
|---|
| 19 | endif
|
|---|
| 20 |
|
|---|
| 21 | c Check for single-element periodicity 'p' bc
|
|---|
| 22 | nz = 1
|
|---|
| 23 | if (if3d) nz = nx
|
|---|
| 24 | call check_p_bc(glo_num,nx,nx,nz,nel)
|
|---|
| 25 |
|
|---|
| 26 | if(nio.eq.0 .and. loglevel.gt.2)
|
|---|
| 27 | $ write(6,*) 'call usrsetvert'
|
|---|
| 28 | call usrsetvert(glo_num,nel,nx,nx,nz)
|
|---|
| 29 | if(nio.eq.0 .and. loglevel.gt.2)
|
|---|
| 30 | $ write(6,'(A,/)') ' done :: usrsetvert'
|
|---|
| 31 |
|
|---|
| 32 | return
|
|---|
| 33 | end
|
|---|
| 34 | c
|
|---|
| 35 | c-----------------------------------------------------------------------
|
|---|
| 36 | c
|
|---|
| 37 | subroutine crs_solve_l2(uf,vf)
|
|---|
| 38 | c
|
|---|
| 39 | c Given an input vector v, this generates the H1 coarse-grid solution
|
|---|
| 40 | c
|
|---|
| 41 | include 'SIZE'
|
|---|
| 42 | include 'DOMAIN'
|
|---|
| 43 | include 'ESOLV'
|
|---|
| 44 | include 'GEOM'
|
|---|
| 45 | include 'PARALLEL'
|
|---|
| 46 | include 'SOLN'
|
|---|
| 47 | include 'INPUT'
|
|---|
| 48 | include 'TSTEP'
|
|---|
| 49 | real uf(1),vf(1)
|
|---|
| 50 | common /scrpre/ uc(lcr*lelt),w(2*lx1*ly1*lz1)
|
|---|
| 51 |
|
|---|
| 52 | call map_f_to_c_l2_bilin(uf,vf,w)
|
|---|
| 53 | call fgslib_crs_solve(xxth(ifield),uc,uf)
|
|---|
| 54 | call map_c_to_f_l2_bilin(uf,uc,w)
|
|---|
| 55 |
|
|---|
| 56 | return
|
|---|
| 57 | end
|
|---|
| 58 | c
|
|---|
| 59 | c-----------------------------------------------------------------------
|
|---|
| 60 | c subroutine test_h1_crs
|
|---|
| 61 | c include 'SIZE'
|
|---|
| 62 | c include 'DOMAIN'
|
|---|
| 63 | c common /scrxxt/ x(lcr*lelv),b(lcr*lelv)
|
|---|
| 64 | c common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
|
|---|
| 65 | c real x,b
|
|---|
| 66 | c ntot=nelv*nxyz_c
|
|---|
| 67 | c do i=1,12
|
|---|
| 68 | c call rzero(b,ntot)
|
|---|
| 69 | c if(mp.eq.1) then
|
|---|
| 70 | c b(i)=1
|
|---|
| 71 | c else if(mid.eq.0) then
|
|---|
| 72 | c if(i.gt.8) b(i-8)=1
|
|---|
| 73 | c else
|
|---|
| 74 | c if(i.le.8) b(i)=1
|
|---|
| 75 | c endif
|
|---|
| 76 | c call hsmg_coarse_solve(x,b)
|
|---|
| 77 | c print *, 'Column ',i,':',(x(j),j=1,ntot)
|
|---|
| 78 | c enddo
|
|---|
| 79 | c return
|
|---|
| 80 | c end
|
|---|
| 81 | c-----------------------------------------------------------------------
|
|---|
| 82 | c
|
|---|
| 83 | subroutine set_up_h1_crs
|
|---|
| 84 |
|
|---|
| 85 | include 'SIZE'
|
|---|
| 86 | include 'GEOM'
|
|---|
| 87 | include 'DOMAIN'
|
|---|
| 88 | include 'INPUT'
|
|---|
| 89 | include 'PARALLEL'
|
|---|
| 90 | include 'TSTEP'
|
|---|
| 91 | common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
|
|---|
| 92 |
|
|---|
| 93 | common /ivrtx/ vertex ((2**ldim)*lelt)
|
|---|
| 94 | integer vertex
|
|---|
| 95 |
|
|---|
| 96 | integer gs_handle
|
|---|
| 97 | integer null_space,e
|
|---|
| 98 |
|
|---|
| 99 | character*3 cb
|
|---|
| 100 | common /scrxxt/ cmlt(lcr,lelv),mask(lcr,lelv)
|
|---|
| 101 | common /scrxxti/ ia(lcr,lcr,lelv), ja(lcr,lcr,lelv)
|
|---|
| 102 | real mask
|
|---|
| 103 | integer ia,ja
|
|---|
| 104 | real z
|
|---|
| 105 |
|
|---|
| 106 | integer key(2),aa(2)
|
|---|
| 107 | common /scrch/ iwork(2,lx1*ly1*lz1*lelv)
|
|---|
| 108 | common /scrns/ w(7*lx1*ly1*lz1*lelv)
|
|---|
| 109 | common /vptsol/ a(27*lx1*ly1*lz1*lelv)
|
|---|
| 110 | integer w
|
|---|
| 111 | real wr(1)
|
|---|
| 112 | equivalence (wr,w)
|
|---|
| 113 |
|
|---|
| 114 | common /scrvhx/ h1(lx1*ly1*lz1*lelv),h2(lx1*ly1*lz1*lelv)
|
|---|
| 115 | common /scrmgx/ w1(lx1*ly1*lz1*lelv),w2(lx1*ly1*lz1*lelv)
|
|---|
| 116 |
|
|---|
| 117 | integer*8 ngv
|
|---|
| 118 | character*132 amgfile_c
|
|---|
| 119 | character*1 fname1(132)
|
|---|
| 120 | equivalence (fname1,amgfile_c)
|
|---|
| 121 | integer nnamg
|
|---|
| 122 |
|
|---|
| 123 | t0 = dnekclock()
|
|---|
| 124 |
|
|---|
| 125 | c nxc is order of coarse grid space + 1, nxc=2, linear, 3=quad,etc.
|
|---|
| 126 | c nxc=param(82)
|
|---|
| 127 | c if (nxc.gt.lxc) then
|
|---|
| 128 | c nxc=lxc
|
|---|
| 129 | c write(6,*) 'WARNING :: coarse grid space too large',nxc,lxc
|
|---|
| 130 | c endif
|
|---|
| 131 | c if (nxc.lt.2) nxc=2
|
|---|
| 132 |
|
|---|
| 133 | nxc = 2
|
|---|
| 134 | nx_crs = nxc
|
|---|
| 135 |
|
|---|
| 136 | if(nio.eq.0) write(6,*) 'setup h1 coarse grid, nx_crs=', nx_crs
|
|---|
| 137 |
|
|---|
| 138 | ncr = nxc**ldim
|
|---|
| 139 | nxyz_c = ncr
|
|---|
| 140 | c
|
|---|
| 141 | c Set SEM_to_GLOB
|
|---|
| 142 | c
|
|---|
| 143 | call get_vertex
|
|---|
| 144 | call set_vert(se_to_gcrs,ngv,nxc,nelv,vertex,.true.)
|
|---|
| 145 |
|
|---|
| 146 | c Set mask
|
|---|
| 147 | z=0
|
|---|
| 148 | ntot=nelv*nxyz_c
|
|---|
| 149 | nzc=1
|
|---|
| 150 | if (if3d) nzc=nxc
|
|---|
| 151 | call rone(mask,ntot)
|
|---|
| 152 | call rone(cmlt,ntot)
|
|---|
| 153 | nfaces=2*ldim
|
|---|
| 154 | c ifield=1 !c? avo: set in set_overlap through 'TSTEP'?
|
|---|
| 155 |
|
|---|
| 156 | if (ifield.eq.1) then
|
|---|
| 157 | do ie=1,nelv
|
|---|
| 158 | do iface=1,nfaces
|
|---|
| 159 | cb=cbc(iface,ie,ifield)
|
|---|
| 160 | if (cb.eq.'o ' .or. cb.eq.'on ' .or.
|
|---|
| 161 | $ cb.eq.'O ' .or. cb.eq.'ON ' .or. cb.eq.'MM ' .or.
|
|---|
| 162 | $ cb.eq.'mm ' .or. cb.eq.'ms ' .or. cb.eq.'MS ')
|
|---|
| 163 | $ call facev(mask,ie,iface,z,nxc,nxc,nzc) ! 'S* ' & 's* ' ?avo?
|
|---|
| 164 | enddo
|
|---|
| 165 | enddo
|
|---|
| 166 | elseif (ifield.eq.ifldmhd) then ! no ifmhd ?avo?
|
|---|
| 167 | do ie=1,nelv
|
|---|
| 168 | do iface=1,nfaces
|
|---|
| 169 | cb=cbc(iface,ie,ifield)
|
|---|
| 170 | if (cb.eq.'ndd' .or. cb.eq.'dnd' .or. cb.eq.'ddn')
|
|---|
| 171 | $ call facev(mask,ie,iface,z,nxc,nxc,nzc)
|
|---|
| 172 | enddo
|
|---|
| 173 | enddo
|
|---|
| 174 | endif
|
|---|
| 175 |
|
|---|
| 176 | c Set global index of dirichlet nodes to zero; xxt will ignore them
|
|---|
| 177 |
|
|---|
| 178 | call fgslib_gs_setup(gs_handle,se_to_gcrs,ntot,nekcomm,mp)
|
|---|
| 179 | call fgslib_gs_op (gs_handle,mask,1,2,0) ! "*"
|
|---|
| 180 | call fgslib_gs_op (gs_handle,cmlt,1,1,0) ! "+"
|
|---|
| 181 | call fgslib_gs_free (gs_handle)
|
|---|
| 182 | call set_jl_crs_mask(ntot,mask,se_to_gcrs)
|
|---|
| 183 |
|
|---|
| 184 | call invcol1(cmlt,ntot)
|
|---|
| 185 |
|
|---|
| 186 | c Setup local SEM-based Neumann operators (for now, just full...)
|
|---|
| 187 |
|
|---|
| 188 | c if (param(51).eq.1) then ! old coarse grid
|
|---|
| 189 | c nxyz1=lx1*ly1*lz1
|
|---|
| 190 | c lda = 27*nxyz1*lelt
|
|---|
| 191 | c ldw = 7*nxyz1*lelt
|
|---|
| 192 | c call get_local_crs(a,lda,nxc,h1,h2,w,ldw)
|
|---|
| 193 | c else
|
|---|
| 194 | c NOTE: a(),h1,...,w2() must all be large enough
|
|---|
| 195 | n = lx1*ly1*lz1*nelv
|
|---|
| 196 | call rone (h1,n)
|
|---|
| 197 | call rzero(h2,n)
|
|---|
| 198 | call get_local_crs_galerkin(a,ncr,nxc,h1,h2,w1,w2)
|
|---|
| 199 | c endif
|
|---|
| 200 |
|
|---|
| 201 | call set_mat_ij(ia,ja,ncr,nelv)
|
|---|
| 202 | null_space=0
|
|---|
| 203 | if (ifield.eq.1) then
|
|---|
| 204 | if (ifvcor) null_space=1
|
|---|
| 205 | elseif (ifield.eq.ifldmhd) then
|
|---|
| 206 | if (ifbcor) null_space=1
|
|---|
| 207 | endif
|
|---|
| 208 |
|
|---|
| 209 | nz=ncr*ncr*nelv
|
|---|
| 210 | isolver = param(40)
|
|---|
| 211 |
|
|---|
| 212 | call blank(fname1,132)
|
|---|
| 213 | lamgn = ltrunc(amgfile,len(amgfile))
|
|---|
| 214 | call chcopy(fname1,amgfile,lamgn)
|
|---|
| 215 | call chcopy(fname1(lamgn+1),char(0),1)
|
|---|
| 216 |
|
|---|
| 217 | ierr = 0
|
|---|
| 218 | call fgslib_crs_setup(xxth(ifield),isolver,nekcomm,mp,ntot,
|
|---|
| 219 | $ se_to_gcrs,nz,ia,ja,a, null_space, crs_param,
|
|---|
| 220 | $ amgfile_c,ierr)
|
|---|
| 221 | ierr = iglmax(ierr,1)
|
|---|
| 222 | if (ifneknek) ierr = iglmax_ms(ierr,1)
|
|---|
| 223 | if (ierr.eq.1) then
|
|---|
| 224 | call exitt
|
|---|
| 225 | endif
|
|---|
| 226 |
|
|---|
| 227 | t0 = dnekclock()-t0
|
|---|
| 228 | if (nio.eq.0) then
|
|---|
| 229 | write(6,*) 'done :: setup h1 coarse grid ',t0, ' sec'
|
|---|
| 230 | write(6,*) ' '
|
|---|
| 231 | endif
|
|---|
| 232 |
|
|---|
| 233 | return
|
|---|
| 234 | end
|
|---|
| 235 | c
|
|---|
| 236 | c-----------------------------------------------------------------------
|
|---|
| 237 | subroutine set_jl_crs_mask(n, mask, se_to_gcrs)
|
|---|
| 238 | real mask(1)
|
|---|
| 239 | integer*8 se_to_gcrs(1)
|
|---|
| 240 | do i=1,n
|
|---|
| 241 | if(mask(i).lt.0.1) se_to_gcrs(i)=0
|
|---|
| 242 | enddo
|
|---|
| 243 | return
|
|---|
| 244 | end
|
|---|
| 245 | c-----------------------------------------------------------------------
|
|---|
| 246 | subroutine set_mat_ij(ia,ja,n,ne)
|
|---|
| 247 | integer n,ne
|
|---|
| 248 | integer ia(n,n,ne), ja(n,n,ne)
|
|---|
| 249 | c
|
|---|
| 250 | integer i,j,ie
|
|---|
| 251 | do ie=1,ne
|
|---|
| 252 | do j=1,n
|
|---|
| 253 | do i=1,n
|
|---|
| 254 | ia(i,j,ie)=(ie-1)*n+i-1
|
|---|
| 255 | ja(i,j,ie)=(ie-1)*n+j-1
|
|---|
| 256 | enddo
|
|---|
| 257 | enddo
|
|---|
| 258 | enddo
|
|---|
| 259 | return
|
|---|
| 260 | end
|
|---|
| 261 | c-----------------------------------------------------------------------
|
|---|
| 262 | subroutine irank_vec(ind,nn,a,m,n,key,nkey,aa)
|
|---|
| 263 | c
|
|---|
| 264 | c Compute rank of each unique entry a(1,i)
|
|---|
| 265 | c
|
|---|
| 266 | c Output: ind(i) i=1,...,n (global) rank of entry a(*,i)
|
|---|
| 267 | c nn = max(rank)
|
|---|
| 268 | c a(j,i) is destroyed
|
|---|
| 269 | c
|
|---|
| 270 | c Input: a(j,i) j=1,...,m; i=1,...,n
|
|---|
| 271 | c m : leading dim. of v (ldv must be .ge. m)
|
|---|
| 272 | c key : sort key
|
|---|
| 273 | c nkey :
|
|---|
| 274 | c
|
|---|
| 275 | c Although not mandatory, this ranking procedure is probably
|
|---|
| 276 | c most effectively employed when the keys are pre-sorted. Thus,
|
|---|
| 277 | c the option is provided to sort vi() prior to the ranking.
|
|---|
| 278 | c
|
|---|
| 279 | c
|
|---|
| 280 | integer ind(n),a(m,n)
|
|---|
| 281 | integer key(nkey),aa(m)
|
|---|
| 282 | logical iftuple_ianeb,a_ne_b
|
|---|
| 283 | c
|
|---|
| 284 | if (m.eq.1) then
|
|---|
| 285 | c
|
|---|
| 286 | write(6,*)
|
|---|
| 287 | $ 'WARNING: For single key, not clear that rank is unique!'
|
|---|
| 288 | call irank(a,ind,n)
|
|---|
| 289 | return
|
|---|
| 290 | endif
|
|---|
| 291 | c
|
|---|
| 292 | c
|
|---|
| 293 | nk = min(nkey,m)
|
|---|
| 294 | call ituple_sort(a,m,n,key,nk,ind,aa)
|
|---|
| 295 | c
|
|---|
| 296 | c Find unique a's
|
|---|
| 297 | c
|
|---|
| 298 | nn=1
|
|---|
| 299 | c
|
|---|
| 300 | call icopy(aa,a,m)
|
|---|
| 301 | a(1,1) = nn
|
|---|
| 302 | a(2,1)=ind(1)
|
|---|
| 303 | c
|
|---|
| 304 | do i=2,n
|
|---|
| 305 | a_ne_b = iftuple_ianeb(aa,a(1,i),key,nk)
|
|---|
| 306 | if (a_ne_b) then
|
|---|
| 307 | call icopy(aa,a(1,i),m)
|
|---|
| 308 | nn = nn+1
|
|---|
| 309 | endif
|
|---|
| 310 | a(1,i) = nn
|
|---|
| 311 | a(2,i) = ind(i)
|
|---|
| 312 | enddo
|
|---|
| 313 | c
|
|---|
| 314 | c Set ind() to rank
|
|---|
| 315 | c
|
|---|
| 316 | do i=1,n
|
|---|
| 317 | iold=a(2,i)
|
|---|
| 318 | ind(iold) = a(1,i)
|
|---|
| 319 | enddo
|
|---|
| 320 | c
|
|---|
| 321 | return
|
|---|
| 322 | end
|
|---|
| 323 | c
|
|---|
| 324 | c-----------------------------------------------------------------------
|
|---|
| 325 | c
|
|---|
| 326 | subroutine ituple_sort(a,lda,n,key,nkey,ind,aa)
|
|---|
| 327 | C
|
|---|
| 328 | C Use Heap Sort (p 231 Num. Rec., 1st Ed.)
|
|---|
| 329 | C
|
|---|
| 330 | integer a(lda,n),aa(lda)
|
|---|
| 331 | integer ind(1),key(nkey)
|
|---|
| 332 | logical iftuple_ialtb
|
|---|
| 333 | C
|
|---|
| 334 | dO 10 j=1,n
|
|---|
| 335 | ind(j)=j
|
|---|
| 336 | 10 continue
|
|---|
| 337 | C
|
|---|
| 338 | if (n.le.1) return
|
|---|
| 339 | L=n/2+1
|
|---|
| 340 | ir=n
|
|---|
| 341 | 100 continue
|
|---|
| 342 | if (l.gt.1) then
|
|---|
| 343 | l=l-1
|
|---|
| 344 | c aa = a (l)
|
|---|
| 345 | call icopy(aa,a(1,l),lda)
|
|---|
| 346 | ii = ind(l)
|
|---|
| 347 | else
|
|---|
| 348 | c aa = a(ir)
|
|---|
| 349 | call icopy(aa,a(1,ir),lda)
|
|---|
| 350 | ii = ind(ir)
|
|---|
| 351 | c a(ir) = a( 1)
|
|---|
| 352 | call icopy(a(1,ir),a(1,1),lda)
|
|---|
| 353 | ind(ir) = ind( 1)
|
|---|
| 354 | ir=ir-1
|
|---|
| 355 | if (ir.eq.1) then
|
|---|
| 356 | c a(1) = aa
|
|---|
| 357 | call icopy(a(1,1),aa,lda)
|
|---|
| 358 | ind(1) = ii
|
|---|
| 359 | return
|
|---|
| 360 | endif
|
|---|
| 361 | endif
|
|---|
| 362 | i=l
|
|---|
| 363 | j=l+l
|
|---|
| 364 | 200 continue
|
|---|
| 365 | if (j.le.ir) then
|
|---|
| 366 | if (j.lt.ir) then
|
|---|
| 367 | if (iftuple_ialtb(a(1,j),a(1,j+1),key,nkey)) j=j+1
|
|---|
| 368 | endif
|
|---|
| 369 | if (iftuple_ialtb(aa,a(1,j),key,nkey)) then
|
|---|
| 370 | c a(i) = a(j)
|
|---|
| 371 | call icopy(a(1,i),a(1,j),lda)
|
|---|
| 372 | ind(i) = ind(j)
|
|---|
| 373 | i=j
|
|---|
| 374 | j=j+j
|
|---|
| 375 | else
|
|---|
| 376 | j=ir+1
|
|---|
| 377 | endif
|
|---|
| 378 | GOTO 200
|
|---|
| 379 | endif
|
|---|
| 380 | c a(i) = aa
|
|---|
| 381 | call icopy(a(1,i),aa,lda)
|
|---|
| 382 | ind(i) = ii
|
|---|
| 383 | GOTO 100
|
|---|
| 384 | end
|
|---|
| 385 | c
|
|---|
| 386 | c-----------------------------------------------------------------------
|
|---|
| 387 | c
|
|---|
| 388 | subroutine tuple_sort(a,lda,n,key,nkey,ind,aa)
|
|---|
| 389 | C
|
|---|
| 390 | C Use Heap Sort (p 231 Num. Rec., 1st Ed.)
|
|---|
| 391 | C
|
|---|
| 392 | real a(lda,n),aa(lda)
|
|---|
| 393 | integer ind(1),key(nkey)
|
|---|
| 394 | logical iftuple_altb
|
|---|
| 395 | C
|
|---|
| 396 | dO 10 j=1,n
|
|---|
| 397 | ind(j)=j
|
|---|
| 398 | 10 continue
|
|---|
| 399 | C
|
|---|
| 400 | if (n.le.1) return
|
|---|
| 401 | L=n/2+1
|
|---|
| 402 | ir=n
|
|---|
| 403 | 100 continue
|
|---|
| 404 | if (l.gt.1) then
|
|---|
| 405 | l=l-1
|
|---|
| 406 | c aa = a (l)
|
|---|
| 407 | call copy(aa,a(1,l),lda)
|
|---|
| 408 | ii = ind(l)
|
|---|
| 409 | else
|
|---|
| 410 | c aa = a(ir)
|
|---|
| 411 | call copy(aa,a(1,ir),lda)
|
|---|
| 412 | ii = ind(ir)
|
|---|
| 413 | c a(ir) = a( 1)
|
|---|
| 414 | call copy(a(1,ir),a(1,1),lda)
|
|---|
| 415 | ind(ir) = ind( 1)
|
|---|
| 416 | ir=ir-1
|
|---|
| 417 | if (ir.eq.1) then
|
|---|
| 418 | c a(1) = aa
|
|---|
| 419 | call copy(a(1,1),aa,lda)
|
|---|
| 420 | ind(1) = ii
|
|---|
| 421 | return
|
|---|
| 422 | endif
|
|---|
| 423 | endif
|
|---|
| 424 | i=l
|
|---|
| 425 | j=l+l
|
|---|
| 426 | 200 continue
|
|---|
| 427 | if (j.le.ir) then
|
|---|
| 428 | if (j.lt.ir) then
|
|---|
| 429 | c if ( a(j).lt.a(j+1) ) j=j+1
|
|---|
| 430 | if (iftuple_altb(a(1,j),a(1,j+1),key,nkey)) j=j+1
|
|---|
| 431 | endif
|
|---|
| 432 | c if (aa.lt.a(j)) then
|
|---|
| 433 | if (iftuple_altb(aa,a(1,j),key,nkey)) then
|
|---|
| 434 | c a(i) = a(j)
|
|---|
| 435 | call copy(a(1,i),a(1,j),lda)
|
|---|
| 436 | ind(i) = ind(j)
|
|---|
| 437 | i=j
|
|---|
| 438 | j=j+j
|
|---|
| 439 | else
|
|---|
| 440 | j=ir+1
|
|---|
| 441 | endif
|
|---|
| 442 | GOTO 200
|
|---|
| 443 | endif
|
|---|
| 444 | c a(i) = aa
|
|---|
| 445 | call copy(a(1,i),aa,lda)
|
|---|
| 446 | ind(i) = ii
|
|---|
| 447 | GOTO 100
|
|---|
| 448 | end
|
|---|
| 449 | c
|
|---|
| 450 | c-----------------------------------------------------------------------
|
|---|
| 451 | c
|
|---|
| 452 | logical function iftuple_ialtb(a,b,key,nkey)
|
|---|
| 453 | integer a(1),b(1)
|
|---|
| 454 | integer key(nkey)
|
|---|
| 455 | c
|
|---|
| 456 | do i=1,nkey
|
|---|
| 457 | k=key(i)
|
|---|
| 458 | if (a(k).lt.b(k)) then
|
|---|
| 459 | iftuple_ialtb = .true.
|
|---|
| 460 | return
|
|---|
| 461 | elseif (a(k).gt.b(k)) then
|
|---|
| 462 | iftuple_ialtb = .false.
|
|---|
| 463 | return
|
|---|
| 464 | endif
|
|---|
| 465 | enddo
|
|---|
| 466 | iftuple_ialtb = .false.
|
|---|
| 467 | return
|
|---|
| 468 | end
|
|---|
| 469 | c
|
|---|
| 470 | c-----------------------------------------------------------------------
|
|---|
| 471 | c
|
|---|
| 472 | logical function iftuple_altb(a,b,key,nkey)
|
|---|
| 473 | real a(1),b(1)
|
|---|
| 474 | integer key(nkey)
|
|---|
| 475 | c
|
|---|
| 476 | do i=1,nkey
|
|---|
| 477 | k=key(i)
|
|---|
| 478 | if (a(k).lt.b(k)) then
|
|---|
| 479 | iftuple_altb = .true.
|
|---|
| 480 | return
|
|---|
| 481 | elseif (a(k).gt.b(k)) then
|
|---|
| 482 | iftuple_altb = .false.
|
|---|
| 483 | return
|
|---|
| 484 | endif
|
|---|
| 485 | enddo
|
|---|
| 486 | iftuple_altb = .false.
|
|---|
| 487 | return
|
|---|
| 488 | end
|
|---|
| 489 | c
|
|---|
| 490 | c-----------------------------------------------------------------------
|
|---|
| 491 | c
|
|---|
| 492 | logical function iftuple_ianeb(a,b,key,nkey)
|
|---|
| 493 | integer a(1),b(1)
|
|---|
| 494 | integer key(nkey)
|
|---|
| 495 | c
|
|---|
| 496 | do i=1,nkey
|
|---|
| 497 | k=key(i)
|
|---|
| 498 | if (a(k).ne.b(k)) then
|
|---|
| 499 | iftuple_ianeb = .true.
|
|---|
| 500 | return
|
|---|
| 501 | endif
|
|---|
| 502 | enddo
|
|---|
| 503 | iftuple_ianeb = .false.
|
|---|
| 504 | return
|
|---|
| 505 | end
|
|---|
| 506 | c
|
|---|
| 507 | c-----------------------------------------------------------------------
|
|---|
| 508 | c
|
|---|
| 509 | subroutine get_local_crs(a,lda,nxc,h1,h2,w,ldw)
|
|---|
| 510 | c
|
|---|
| 511 | c This routine generates Nelv submatrices of order nxc^ldim.
|
|---|
| 512 | c
|
|---|
| 513 | include 'SIZE'
|
|---|
| 514 | include 'GEOM'
|
|---|
| 515 | include 'INPUT'
|
|---|
| 516 | include 'TSTEP'
|
|---|
| 517 | include 'DOMAIN'
|
|---|
| 518 | include 'PARALLEL'
|
|---|
| 519 | c
|
|---|
| 520 | c
|
|---|
| 521 | c Generate local triangular matrix
|
|---|
| 522 | c
|
|---|
| 523 | real a(1),h1(1),h2(1),w(ldw)
|
|---|
| 524 | c
|
|---|
| 525 | parameter (lcrd=lx1**ldim)
|
|---|
| 526 | common /ctmp1/ x(lcrd),y(lcrd),z(lcrd)
|
|---|
| 527 | c
|
|---|
| 528 | c
|
|---|
| 529 | ncrs_loc = nxc**ldim
|
|---|
| 530 | n2 = ncrs_loc*ncrs_loc
|
|---|
| 531 | c
|
|---|
| 532 | c Required storage for a:
|
|---|
| 533 | nda = n2*nelv
|
|---|
| 534 | if (nda.gt.lda) then
|
|---|
| 535 | write(6,*)nid,'ERROR: increase storage get_local_crs:',nda,lda
|
|---|
| 536 | call exitt
|
|---|
| 537 | endif
|
|---|
| 538 | c
|
|---|
| 539 | c
|
|---|
| 540 | l = 1
|
|---|
| 541 | do ie=1,nelv
|
|---|
| 542 | c
|
|---|
| 543 | call map_m_to_n(x,nxc,xm1(1,1,1,ie),lx1,if3d,w,ldw)
|
|---|
| 544 | call map_m_to_n(y,nxc,ym1(1,1,1,ie),lx1,if3d,w,ldw)
|
|---|
| 545 | if (if3d) call map_m_to_n(z,nxc,zm1(1,1,1,ie),lx1,if3d,w,ldw)
|
|---|
| 546 | c.later. call map_m_to_n(hl1,nxc,h1(1,1,1,ie),lx1,if3d,w,ldw)
|
|---|
| 547 | c.later. call map_m_to_n(hl2,nxc,h2(1,1,1,ie),lx1,if3d,w,ldw)
|
|---|
| 548 | c
|
|---|
| 549 | call a_crs_enriched(a(l),h1,h2,x,y,z,nxc,if3d,ie)
|
|---|
| 550 | l=l+n2
|
|---|
| 551 | c
|
|---|
| 552 | enddo
|
|---|
| 553 | c
|
|---|
| 554 | return
|
|---|
| 555 | end
|
|---|
| 556 | c
|
|---|
| 557 | c-----------------------------------------------------------------------
|
|---|
| 558 | c
|
|---|
| 559 | subroutine a_crs_enriched(a,h1,h2,x1,y1,z1,nxc,if3d,ie)
|
|---|
| 560 | c
|
|---|
| 561 | c This sets up a matrix for a single array of tensor-product
|
|---|
| 562 | c gridpoints (e.g., an array defined by SEM-GLL vertices)
|
|---|
| 563 | c
|
|---|
| 564 | c For example, suppose ldim=3.
|
|---|
| 565 | c
|
|---|
| 566 | c Then, there would be ncrs_loc := nxc^3 dofs for this matrix,
|
|---|
| 567 | c
|
|---|
| 568 | c and the matrix size would be (ncrs_loc x ncrs_loc).
|
|---|
| 569 | c
|
|---|
| 570 | c
|
|---|
| 571 | c
|
|---|
| 572 | include 'SIZE'
|
|---|
| 573 | c
|
|---|
| 574 | real a(1),h1(1),h2(1)
|
|---|
| 575 | real x1(nxc,nxc,1),y1(nxc,nxc,1),z1(nxc,nxc,1)
|
|---|
| 576 | logical if3d
|
|---|
| 577 | c
|
|---|
| 578 | parameter (ldm2=2**ldim)
|
|---|
| 579 | real a_loc(ldm2,ldm2)
|
|---|
| 580 | real x(8),y(8),z(8)
|
|---|
| 581 | c
|
|---|
| 582 | ncrs_loc = nxc**ldim
|
|---|
| 583 | n2 = ncrs_loc*ncrs_loc
|
|---|
| 584 | call rzero(a,n2)
|
|---|
| 585 | c
|
|---|
| 586 | nyc=nxc
|
|---|
| 587 | nzc=2
|
|---|
| 588 | if (if3d) nzc=nxc
|
|---|
| 589 | nz =0
|
|---|
| 590 | if (if3d) nz=1
|
|---|
| 591 | c
|
|---|
| 592 | c Here, we march across sub-cubes
|
|---|
| 593 | c
|
|---|
| 594 | do kz=1,nzc-1
|
|---|
| 595 | do ky=1,nyc-1
|
|---|
| 596 | do kx=1,nxc-1
|
|---|
| 597 | k = 0
|
|---|
| 598 | do iz=0,nz
|
|---|
| 599 | do iy=0,1
|
|---|
| 600 | do ix=0,1
|
|---|
| 601 | k = k+1
|
|---|
| 602 | x(k) = x1(kx+ix,ky+iy,kz+iz)
|
|---|
| 603 | y(k) = y1(kx+ix,ky+iy,kz+iz)
|
|---|
| 604 | z(k) = z1(kx+ix,ky+iy,kz+iz)
|
|---|
| 605 | enddo
|
|---|
| 606 | enddo
|
|---|
| 607 | enddo
|
|---|
| 608 | if (if3d) then
|
|---|
| 609 | call a_crs_3d(a_loc,h1,h2,x,y,z,ie)
|
|---|
| 610 | else
|
|---|
| 611 | call a_crs_2d(a_loc,h1,h2,x,y,ie)
|
|---|
| 612 | endif
|
|---|
| 613 | c call outmat(a_loc,ldm2,ldm2,'A_loc ',ie)
|
|---|
| 614 | c
|
|---|
| 615 | c Assemble:
|
|---|
| 616 | c
|
|---|
| 617 | j = 0
|
|---|
| 618 | do jz=0,nz
|
|---|
| 619 | do jy=0,1
|
|---|
| 620 | do jx=0,1
|
|---|
| 621 | j = j+1
|
|---|
| 622 | ja = (kx+jx) + nxc*(ky+jy-1) + nxc*nyc*(kz+jz-1)
|
|---|
| 623 | c
|
|---|
| 624 | i = 0
|
|---|
| 625 | do iz=0,nz
|
|---|
| 626 | do iy=0,1
|
|---|
| 627 | do ix=0,1
|
|---|
| 628 | i = i+1
|
|---|
| 629 | ia = (kx+ix) + nxc*(ky+iy-1) + nxc*nyc*(kz+iz-1)
|
|---|
| 630 | c
|
|---|
| 631 | ija = ia + ncrs_loc*(ja-1)
|
|---|
| 632 | a(ija) = a(ija) + a_loc(i,j)
|
|---|
| 633 | c
|
|---|
| 634 | enddo
|
|---|
| 635 | enddo
|
|---|
| 636 | enddo
|
|---|
| 637 | c
|
|---|
| 638 | enddo
|
|---|
| 639 | enddo
|
|---|
| 640 | enddo
|
|---|
| 641 | enddo
|
|---|
| 642 | enddo
|
|---|
| 643 | enddo
|
|---|
| 644 | c
|
|---|
| 645 | return
|
|---|
| 646 | end
|
|---|
| 647 | c
|
|---|
| 648 | c-----------------------------------------------------------------------
|
|---|
| 649 | c
|
|---|
| 650 | subroutine a_crs_3d(a,h1,h2,xc,yc,zc,ie)
|
|---|
| 651 | c
|
|---|
| 652 | c Generate stiffness matrix for 3D coarse grid problem.
|
|---|
| 653 | c
|
|---|
| 654 | c This is done by using two tetrahedrizations of each
|
|---|
| 655 | c hexahedral subdomain (element) such that each of the
|
|---|
| 656 | c 6 panels (faces) on the sides of an element has a big X.
|
|---|
| 657 | c
|
|---|
| 658 | c
|
|---|
| 659 | real a(0:7,0:7),h1(0:7),h2(0:7)
|
|---|
| 660 | real xc(0:7),yc(0:7),zc(0:7)
|
|---|
| 661 | c
|
|---|
| 662 | real a_loc(4,4)
|
|---|
| 663 | real xt(4),yt(4),zt(4)
|
|---|
| 664 | c
|
|---|
| 665 | integer vertex(4,5,2)
|
|---|
| 666 | save vertex
|
|---|
| 667 | data vertex / 000 , 001 , 010 , 100
|
|---|
| 668 | $ , 000 , 001 , 011 , 101
|
|---|
| 669 | $ , 011 , 010 , 000 , 110
|
|---|
| 670 | $ , 011 , 010 , 001 , 111
|
|---|
| 671 | $ , 000 , 110 , 101 , 011
|
|---|
| 672 | c
|
|---|
| 673 | $ , 101 , 100 , 110 , 000
|
|---|
| 674 | $ , 101 , 100 , 111 , 001
|
|---|
| 675 | $ , 110 , 111 , 100 , 010
|
|---|
| 676 | $ , 110 , 111 , 101 , 011
|
|---|
| 677 | $ , 111 , 001 , 100 , 010 /
|
|---|
| 678 | c
|
|---|
| 679 | integer icalld
|
|---|
| 680 | save icalld
|
|---|
| 681 | data icalld/0/
|
|---|
| 682 | c
|
|---|
| 683 | if (icalld.eq.0) then
|
|---|
| 684 | do i=1,40
|
|---|
| 685 | call bindec(vertex(i,1,1))
|
|---|
| 686 | enddo
|
|---|
| 687 | endif
|
|---|
| 688 | icalld=icalld+1
|
|---|
| 689 | c
|
|---|
| 690 | call rzero(a,64)
|
|---|
| 691 | do k=1,10
|
|---|
| 692 | do iv=1,4
|
|---|
| 693 | xt(iv) = xc(vertex(iv,k,1))
|
|---|
| 694 | yt(iv) = yc(vertex(iv,k,1))
|
|---|
| 695 | zt(iv) = zc(vertex(iv,k,1))
|
|---|
| 696 | enddo
|
|---|
| 697 | call get_local_A_tet(a_loc,xt,yt,zt,k,ie)
|
|---|
| 698 | do j=1,4
|
|---|
| 699 | jj = vertex(j,k,1)
|
|---|
| 700 | do i=1,4
|
|---|
| 701 | ii = vertex(i,k,1)
|
|---|
| 702 | a(ii,jj) = a(ii,jj) + 0.5*a_loc(i,j)
|
|---|
| 703 | enddo
|
|---|
| 704 | enddo
|
|---|
| 705 | enddo
|
|---|
| 706 | return
|
|---|
| 707 | end
|
|---|
| 708 | c
|
|---|
| 709 | c-----------------------------------------------------------------------
|
|---|
| 710 | c
|
|---|
| 711 | subroutine bindec(bin_in)
|
|---|
| 712 | integer bin_in,d,b,b2
|
|---|
| 713 | c
|
|---|
| 714 | keep = bin_in
|
|---|
| 715 | d = bin_in
|
|---|
| 716 | b2 = 1
|
|---|
| 717 | b = 0
|
|---|
| 718 | do l=1,12
|
|---|
| 719 | b = b + b2*mod(d,10)
|
|---|
| 720 | d = d/10
|
|---|
| 721 | b2 = b2*2
|
|---|
| 722 | if (d.eq.0) goto 1
|
|---|
| 723 | enddo
|
|---|
| 724 | 1 continue
|
|---|
| 725 | bin_in = b
|
|---|
| 726 | return
|
|---|
| 727 | end
|
|---|
| 728 | c
|
|---|
| 729 | c-----------------------------------------------------------------------
|
|---|
| 730 | subroutine get_local_A_tet(a,x,y,z,kt,ie)
|
|---|
| 731 | c
|
|---|
| 732 | c Generate local tetrahedral matrix
|
|---|
| 733 | c
|
|---|
| 734 | c
|
|---|
| 735 | real a(4,4), g(4,4)
|
|---|
| 736 | real x(4),y(4),z(4)
|
|---|
| 737 | c
|
|---|
| 738 | 11 continue
|
|---|
| 739 | x23 = x(2) - x(3)
|
|---|
| 740 | y23 = y(2) - y(3)
|
|---|
| 741 | z23 = z(2) - z(3)
|
|---|
| 742 | x34 = x(3) - x(4)
|
|---|
| 743 | y34 = y(3) - y(4)
|
|---|
| 744 | z34 = z(3) - z(4)
|
|---|
| 745 | x41 = x(4) - x(1)
|
|---|
| 746 | y41 = y(4) - y(1)
|
|---|
| 747 | z41 = z(4) - z(1)
|
|---|
| 748 | x12 = x(1) - x(2)
|
|---|
| 749 | y12 = y(1) - y(2)
|
|---|
| 750 | z12 = z(1) - z(2)
|
|---|
| 751 | c
|
|---|
| 752 | xy234 = x34*y23 - x23*y34
|
|---|
| 753 | xy341 = x34*y41 - x41*y34
|
|---|
| 754 | xy412 = x12*y41 - x41*y12
|
|---|
| 755 | xy123 = x12*y23 - x23*y12
|
|---|
| 756 | xz234 = x23*z34 - x34*z23
|
|---|
| 757 | xz341 = x41*z34 - x34*z41
|
|---|
| 758 | xz412 = x41*z12 - x12*z41
|
|---|
| 759 | xz123 = x23*z12 - x12*z23
|
|---|
| 760 | yz234 = y34*z23 - y23*z34
|
|---|
| 761 | yz341 = y34*z41 - y41*z34
|
|---|
| 762 | yz412 = y12*z41 - y41*z12
|
|---|
| 763 | yz123 = y12*z23 - y23*z12
|
|---|
| 764 | c
|
|---|
| 765 | g(1,1) = -(x(2)*yz234 + y(2)*xz234 + z(2)*xy234)
|
|---|
| 766 | g(2,1) = -(x(3)*yz341 + y(3)*xz341 + z(3)*xy341)
|
|---|
| 767 | g(3,1) = -(x(4)*yz412 + y(4)*xz412 + z(4)*xy412)
|
|---|
| 768 | g(4,1) = -(x(1)*yz123 + y(1)*xz123 + z(1)*xy123)
|
|---|
| 769 | g(1,2) = yz234
|
|---|
| 770 | g(2,2) = yz341
|
|---|
| 771 | g(3,2) = yz412
|
|---|
| 772 | g(4,2) = yz123
|
|---|
| 773 | g(1,3) = xz234
|
|---|
| 774 | g(2,3) = xz341
|
|---|
| 775 | g(3,3) = xz412
|
|---|
| 776 | g(4,3) = xz123
|
|---|
| 777 | g(1,4) = xy234
|
|---|
| 778 | g(2,4) = xy341
|
|---|
| 779 | g(3,4) = xy412
|
|---|
| 780 | g(4,4) = xy123
|
|---|
| 781 | c
|
|---|
| 782 | c vol36 = 1/(36*volume) = 1/(6*determinant)
|
|---|
| 783 | c
|
|---|
| 784 | det = x(1)*yz234 + x(2)*yz341 + x(3)*yz412 + x(4)*yz123
|
|---|
| 785 | vol36 = 1.0/(6.0*det)
|
|---|
| 786 | if (vol36.lt.0) then
|
|---|
| 787 | write(6,*) 'Error: tetrahedron not right-handed',ie
|
|---|
| 788 | write(6,1) 'x',(x(k),k=1,4)
|
|---|
| 789 | write(6,1) 'y',(y(k),k=1,4)
|
|---|
| 790 | write(6,1) 'z',(z(k),k=1,4)
|
|---|
| 791 | 1 format(a1,1p4e15.5)
|
|---|
| 792 |
|
|---|
| 793 | c call exitt ! Option 1
|
|---|
| 794 |
|
|---|
| 795 | xx = x(1) ! Option 2
|
|---|
| 796 | x(1) = x(2) ! -- this is the option that
|
|---|
| 797 | x(2) = xx ! actually works. 11/25/07
|
|---|
| 798 |
|
|---|
| 799 | xx = y(1)
|
|---|
| 800 | y(1) = y(2)
|
|---|
| 801 | y(2) = xx
|
|---|
| 802 |
|
|---|
| 803 | xx = z(1)
|
|---|
| 804 | z(1) = z(2)
|
|---|
| 805 | z(2) = xx
|
|---|
| 806 |
|
|---|
| 807 | goto 11
|
|---|
| 808 |
|
|---|
| 809 | c call rzero(a,16) ! Option 3
|
|---|
| 810 | c return
|
|---|
| 811 |
|
|---|
| 812 | c vol36 = abs(vol36) ! Option 4
|
|---|
| 813 |
|
|---|
| 814 | endif
|
|---|
| 815 | c
|
|---|
| 816 | do j=1,4
|
|---|
| 817 | do i=1,4
|
|---|
| 818 | a(i,j)=vol36*(g(i,2)*g(j,2)+g(i,3)*g(j,3)+g(i,4)*g(j,4))
|
|---|
| 819 | enddo
|
|---|
| 820 | enddo
|
|---|
| 821 | c
|
|---|
| 822 | return
|
|---|
| 823 | end
|
|---|
| 824 | c
|
|---|
| 825 | c-----------------------------------------------------------------------
|
|---|
| 826 | c
|
|---|
| 827 | subroutine a_crs_2d(a,h1,h2,x,y,ie)
|
|---|
| 828 | c
|
|---|
| 829 | include 'SIZE'
|
|---|
| 830 | include 'GEOM'
|
|---|
| 831 | include 'INPUT'
|
|---|
| 832 | include 'TSTEP'
|
|---|
| 833 | include 'DOMAIN'
|
|---|
| 834 | include 'PARALLEL'
|
|---|
| 835 | c
|
|---|
| 836 | c Generate local triangle-based stiffnes matrix for quad
|
|---|
| 837 | c
|
|---|
| 838 | real a(4,4),h1(1),h2(1)
|
|---|
| 839 | real x(1),y(1)
|
|---|
| 840 | c
|
|---|
| 841 | c Triangle to Square pointers
|
|---|
| 842 | c
|
|---|
| 843 | integer elem(3,2)
|
|---|
| 844 | save elem
|
|---|
| 845 | data elem / 1,2,4 , 1,4,3 /
|
|---|
| 846 | c
|
|---|
| 847 | real a_loc(3,3)
|
|---|
| 848 | c
|
|---|
| 849 | c
|
|---|
| 850 | call rzero(a,16)
|
|---|
| 851 | c
|
|---|
| 852 | do i=1,2
|
|---|
| 853 | j1 = elem(1,i)
|
|---|
| 854 | j2 = elem(2,i)
|
|---|
| 855 | j3 = elem(3,i)
|
|---|
| 856 | x1=x(j1)
|
|---|
| 857 | y1=y(j1)
|
|---|
| 858 | x2=x(j2)
|
|---|
| 859 | y2=y(j2)
|
|---|
| 860 | x3=x(j3)
|
|---|
| 861 | y3=y(j3)
|
|---|
| 862 | c
|
|---|
| 863 | y23=y2-y3
|
|---|
| 864 | y31=y3-y1
|
|---|
| 865 | y12=y1-y2
|
|---|
| 866 | c
|
|---|
| 867 | x32=x3-x2
|
|---|
| 868 | x13=x1-x3
|
|---|
| 869 | x21=x2-x1
|
|---|
| 870 | c
|
|---|
| 871 | c area4 = 1/(4*area)
|
|---|
| 872 | area4 = 0.50/(x21*y31 - y12*x13)
|
|---|
| 873 | c
|
|---|
| 874 | a_loc(1, 1) = area4*( y23*y23+x32*x32 )
|
|---|
| 875 | a_loc(1, 2) = area4*( y23*y31+x32*x13 )
|
|---|
| 876 | a_loc(1, 3) = area4*( y23*y12+x32*x21 )
|
|---|
| 877 | c
|
|---|
| 878 | a_loc(2, 1) = area4*( y31*y23+x13*x32 )
|
|---|
| 879 | a_loc(2, 2) = area4*( y31*y31+x13*x13 )
|
|---|
| 880 | a_loc(2, 3) = area4*( y31*y12+x13*x21 )
|
|---|
| 881 | c
|
|---|
| 882 | a_loc(3, 1) = area4*( y12*y23+x21*x32 )
|
|---|
| 883 | a_loc(3, 2) = area4*( y12*y31+x21*x13 )
|
|---|
| 884 | a_loc(3, 3) = area4*( y12*y12+x21*x21 )
|
|---|
| 885 | c
|
|---|
| 886 | c Store in "4 x 4" format
|
|---|
| 887 | c
|
|---|
| 888 | do il=1,3
|
|---|
| 889 | iv = elem(il,i)
|
|---|
| 890 | do jl=1,3
|
|---|
| 891 | jv = elem(jl,i)
|
|---|
| 892 | a(iv,jv) = a(iv,jv) + a_loc(il,jl)
|
|---|
| 893 | enddo
|
|---|
| 894 | enddo
|
|---|
| 895 | enddo
|
|---|
| 896 | c
|
|---|
| 897 | return
|
|---|
| 898 | end
|
|---|
| 899 | c
|
|---|
| 900 | c-----------------------------------------------------------------------
|
|---|
| 901 | c
|
|---|
| 902 | subroutine map_m_to_n(a,na,b,nb,if3d,w,ldw)
|
|---|
| 903 | c
|
|---|
| 904 | c Input: b
|
|---|
| 905 | c Output: a
|
|---|
| 906 | c
|
|---|
| 907 | real a(1),b(1),w(1)
|
|---|
| 908 | logical if3d
|
|---|
| 909 | c
|
|---|
| 910 | parameter(lx=50)
|
|---|
| 911 | real za(lx),zb(lx)
|
|---|
| 912 | c
|
|---|
| 913 | real iba(lx*lx),ibat(lx*lx)
|
|---|
| 914 | save iba,ibat
|
|---|
| 915 | c
|
|---|
| 916 | integer nao,nbo
|
|---|
| 917 | save nao,nbo
|
|---|
| 918 | data nao,nbo / -9, -9/
|
|---|
| 919 | c
|
|---|
| 920 | if (na.gt.lx.or.nb.gt.lx) then
|
|---|
| 921 | write(6,*)'ERROR: increase lx in map_m_to_n to max:',na,nb
|
|---|
| 922 | call exitt
|
|---|
| 923 | endif
|
|---|
| 924 | c
|
|---|
| 925 | if (na.ne.nao .or. nb.ne.nbo) then
|
|---|
| 926 | nao = na
|
|---|
| 927 | nbo = nb
|
|---|
| 928 | call zwgll(za,w,na)
|
|---|
| 929 | call zwgll(zb,w,nb)
|
|---|
| 930 | call igllm(iba,ibat,zb,za,nb,na,nb,na)
|
|---|
| 931 | endif
|
|---|
| 932 | c
|
|---|
| 933 | call specmpn(a,na,b,nb,iba,ibat,if3d,w,ldw)
|
|---|
| 934 | c
|
|---|
| 935 | return
|
|---|
| 936 | end
|
|---|
| 937 | c
|
|---|
| 938 | c-----------------------------------------------------------------------
|
|---|
| 939 | c
|
|---|
| 940 | subroutine specmpn(b,nb,a,na,ba,ab,if3d,w,ldw)
|
|---|
| 941 | C
|
|---|
| 942 | C - Spectral interpolation from A to B via tensor products
|
|---|
| 943 | C - scratch arrays: w(na*na*nb + nb*nb*na)
|
|---|
| 944 | C
|
|---|
| 945 | C 5/3/00 -- this routine replaces specmp in navier1.f, which
|
|---|
| 946 | c has a potential memory problem
|
|---|
| 947 | C
|
|---|
| 948 | C
|
|---|
| 949 | logical if3d
|
|---|
| 950 | c
|
|---|
| 951 | real b(nb,nb,nb),a(na,na,na)
|
|---|
| 952 | real w(ldw)
|
|---|
| 953 | c
|
|---|
| 954 | ltest = na*nb
|
|---|
| 955 | if (if3d) ltest = na*na*nb + nb*na*na
|
|---|
| 956 | if (ldw.lt.ltest) then
|
|---|
| 957 | write(6,*) 'ERROR specmp:',ldw,ltest,if3d
|
|---|
| 958 | call exitt
|
|---|
| 959 | endif
|
|---|
| 960 | c
|
|---|
| 961 | if (if3d) then
|
|---|
| 962 | nab = na*nb
|
|---|
| 963 | nbb = nb*nb
|
|---|
| 964 | call mxm(ba,nb,a,na,w,na*na)
|
|---|
| 965 | k=1
|
|---|
| 966 | l=na*na*nb + 1
|
|---|
| 967 | do iz=1,na
|
|---|
| 968 | call mxm(w(k),nb,ab,na,w(l),nb)
|
|---|
| 969 | k=k+nab
|
|---|
| 970 | l=l+nbb
|
|---|
| 971 | enddo
|
|---|
| 972 | l=na*na*nb + 1
|
|---|
| 973 | call mxm(w(l),nbb,ab,na,b,nb)
|
|---|
| 974 | else
|
|---|
| 975 | call mxm(ba,nb,a,na,w,na)
|
|---|
| 976 | call mxm(w,nb,ab,na,b,nb)
|
|---|
| 977 | endif
|
|---|
| 978 | return
|
|---|
| 979 | end
|
|---|
| 980 | c
|
|---|
| 981 | c-----------------------------------------------------------------------
|
|---|
| 982 | c
|
|---|
| 983 | subroutine irank(A,IND,N)
|
|---|
| 984 | C
|
|---|
| 985 | C Use Heap Sort (p 233 Num. Rec.), 5/26/93 pff.
|
|---|
| 986 | C
|
|---|
| 987 | integer A(1),IND(1)
|
|---|
| 988 | C
|
|---|
| 989 | if (n.le.1) return
|
|---|
| 990 | DO 10 J=1,N
|
|---|
| 991 | IND(j)=j
|
|---|
| 992 | 10 continue
|
|---|
| 993 | C
|
|---|
| 994 | if (n.eq.1) return
|
|---|
| 995 | L=n/2+1
|
|---|
| 996 | ir=n
|
|---|
| 997 | 100 continue
|
|---|
| 998 | IF (l.gt.1) THEN
|
|---|
| 999 | l=l-1
|
|---|
| 1000 | indx=ind(l)
|
|---|
| 1001 | q=a(indx)
|
|---|
| 1002 | ELSE
|
|---|
| 1003 | indx=ind(ir)
|
|---|
| 1004 | q=a(indx)
|
|---|
| 1005 | ind(ir)=ind(1)
|
|---|
| 1006 | ir=ir-1
|
|---|
| 1007 | if (ir.eq.1) then
|
|---|
| 1008 | ind(1)=indx
|
|---|
| 1009 | return
|
|---|
| 1010 | endif
|
|---|
| 1011 | ENDIF
|
|---|
| 1012 | i=l
|
|---|
| 1013 | j=l+l
|
|---|
| 1014 | 200 continue
|
|---|
| 1015 | IF (J.le.IR) THEN
|
|---|
| 1016 | IF (J.lt.IR) THEN
|
|---|
| 1017 | IF ( A(IND(j)).lt.A(IND(j+1)) ) j=j+1
|
|---|
| 1018 | ENDIF
|
|---|
| 1019 | IF (q.lt.A(IND(j))) THEN
|
|---|
| 1020 | IND(I)=IND(J)
|
|---|
| 1021 | I=J
|
|---|
| 1022 | J=J+J
|
|---|
| 1023 | ELSE
|
|---|
| 1024 | J=IR+1
|
|---|
| 1025 | ENDIF
|
|---|
| 1026 | GOTO 200
|
|---|
| 1027 | ENDIF
|
|---|
| 1028 | IND(I)=INDX
|
|---|
| 1029 | GOTO 100
|
|---|
| 1030 | END
|
|---|
| 1031 | c-----------------------------------------------------------------------
|
|---|
| 1032 | subroutine iranku(r,input,n,w,ind)
|
|---|
| 1033 | c
|
|---|
| 1034 | c Return the rank of each input value, and the maximum rank.
|
|---|
| 1035 | c
|
|---|
| 1036 | c OUTPUT: r(k) = rank of each entry, k=1,..,n
|
|---|
| 1037 | c maxr = max( r )
|
|---|
| 1038 | c w(i) = sorted & compressed list of input values
|
|---|
| 1039 | c
|
|---|
| 1040 | integer r(1),input(1),ind(1),w(1)
|
|---|
| 1041 | c
|
|---|
| 1042 | call icopy(r,input,n)
|
|---|
| 1043 | call isort(r,ind,n)
|
|---|
| 1044 | c
|
|---|
| 1045 | maxr = 1
|
|---|
| 1046 | rlast = r(1)
|
|---|
| 1047 | do i=1,n
|
|---|
| 1048 | c Bump rank only when r_i changes
|
|---|
| 1049 | if (r(i).ne.rlast) then
|
|---|
| 1050 | rlast = r(i)
|
|---|
| 1051 | maxr = maxr + 1
|
|---|
| 1052 | endif
|
|---|
| 1053 | r(i) = maxr
|
|---|
| 1054 | enddo
|
|---|
| 1055 | call iunswap(r,ind,n,w)
|
|---|
| 1056 | c
|
|---|
| 1057 | return
|
|---|
| 1058 | end
|
|---|
| 1059 | c
|
|---|
| 1060 | c-----------------------------------------------------------------------
|
|---|
| 1061 | c
|
|---|
| 1062 | subroutine ifacev_redef(a,ie,iface,val,nx,ny,nz)
|
|---|
| 1063 | C
|
|---|
| 1064 | C Assign the value VAL to face(IFACE,IE) of array A.
|
|---|
| 1065 | C IFACE is the input in the pre-processor ordering scheme.
|
|---|
| 1066 | C
|
|---|
| 1067 | include 'SIZE'
|
|---|
| 1068 | integer a(nx,ny,nz,lelt),val
|
|---|
| 1069 | call facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
|
|---|
| 1070 | do 100 iz=kz1,kz2
|
|---|
| 1071 | do 100 iy=ky1,ky2
|
|---|
| 1072 | do 100 ix=kx1,kx2
|
|---|
| 1073 | a(ix,iy,iz,ie)=val
|
|---|
| 1074 | 100 continue
|
|---|
| 1075 | return
|
|---|
| 1076 | end
|
|---|
| 1077 | c
|
|---|
| 1078 | c-----------------------------------------------------------------------
|
|---|
| 1079 | subroutine map_c_to_f_l2_bilin(uf,uc,w)
|
|---|
| 1080 | c
|
|---|
| 1081 | c H1 Iterpolation operator: linear --> spectral GLL mesh
|
|---|
| 1082 | c
|
|---|
| 1083 | include 'SIZE'
|
|---|
| 1084 | include 'DOMAIN'
|
|---|
| 1085 | include 'INPUT'
|
|---|
| 1086 | c
|
|---|
| 1087 | parameter (lxyz = lx2*ly2*lz2)
|
|---|
| 1088 | real uc(nxyz_c,lelt),uf(lxyz,lelt),w(1)
|
|---|
| 1089 |
|
|---|
| 1090 | ltot22 = 2*lx2*ly2*lz2
|
|---|
| 1091 | nx_crs = 2 ! bilinear only
|
|---|
| 1092 |
|
|---|
| 1093 | do ie=1,nelv
|
|---|
| 1094 | call maph1_to_l2(uf(1,ie),lx2,uc(1,ie),nx_crs,if3d,w,ltot22)
|
|---|
| 1095 | enddo
|
|---|
| 1096 | c
|
|---|
| 1097 | return
|
|---|
| 1098 | end
|
|---|
| 1099 | c
|
|---|
| 1100 | c-----------------------------------------------------------------------
|
|---|
| 1101 | c
|
|---|
| 1102 | subroutine map_f_to_c_l2_bilin(uc,uf,w)
|
|---|
| 1103 |
|
|---|
| 1104 | c TRANSPOSE of L2 Iterpolation operator: T
|
|---|
| 1105 | c (linear --> spectral GLL mesh)
|
|---|
| 1106 |
|
|---|
| 1107 | include 'SIZE'
|
|---|
| 1108 | include 'DOMAIN'
|
|---|
| 1109 | include 'INPUT'
|
|---|
| 1110 |
|
|---|
| 1111 | parameter (lxyz = lx2*ly2*lz2)
|
|---|
| 1112 | real uc(nxyz_c,lelt),uf(lxyz,lelt),w(1)
|
|---|
| 1113 |
|
|---|
| 1114 | ltot22 = 2*lx2*ly2*lz2
|
|---|
| 1115 | nx_crs = 2 ! bilinear only
|
|---|
| 1116 |
|
|---|
| 1117 | do ie=1,nelv
|
|---|
| 1118 | call maph1_to_l2t(uc(1,ie),nx_crs,uf(1,ie),lx2,if3d,w,ltot22)
|
|---|
| 1119 | enddo
|
|---|
| 1120 | c
|
|---|
| 1121 | return
|
|---|
| 1122 | end
|
|---|
| 1123 | c
|
|---|
| 1124 | c-----------------------------------------------------------------------
|
|---|
| 1125 | c
|
|---|
| 1126 | subroutine maph1_to_l2(a,na,b,nb,if3d,w,ldw)
|
|---|
| 1127 | c
|
|---|
| 1128 | c Input: b
|
|---|
| 1129 | c Output: a
|
|---|
| 1130 | c
|
|---|
| 1131 | real a(1),b(1),w(1)
|
|---|
| 1132 | logical if3d
|
|---|
| 1133 | c
|
|---|
| 1134 | parameter(lx=50)
|
|---|
| 1135 | real za(lx),zb(lx)
|
|---|
| 1136 | c
|
|---|
| 1137 | real iba(lx*lx),ibat(lx*lx)
|
|---|
| 1138 | save iba,ibat
|
|---|
| 1139 | c
|
|---|
| 1140 | integer nao,nbo
|
|---|
| 1141 | save nao,nbo
|
|---|
| 1142 | data nao,nbo / -9, -9/
|
|---|
| 1143 | c
|
|---|
| 1144 | c
|
|---|
| 1145 | if (na.gt.lx.or.nb.gt.lx) then
|
|---|
| 1146 | write(6,*)'ERROR: increase lx in maph1_to_l2 to max:',na,nb
|
|---|
| 1147 | call exitt
|
|---|
| 1148 | endif
|
|---|
| 1149 | c
|
|---|
| 1150 | if (na.ne.nao .or. nb.ne.nbo) then
|
|---|
| 1151 | nao = na
|
|---|
| 1152 | nbo = nb
|
|---|
| 1153 | call zwgl (za,w,na)
|
|---|
| 1154 | call zwgll(zb,w,nb)
|
|---|
| 1155 | call igllm(iba,ibat,zb,za,nb,na,nb,na)
|
|---|
| 1156 | endif
|
|---|
| 1157 | c
|
|---|
| 1158 | call specmpn(a,na,b,nb,iba,ibat,if3d,w,ldw)
|
|---|
| 1159 | c
|
|---|
| 1160 | return
|
|---|
| 1161 | end
|
|---|
| 1162 | c
|
|---|
| 1163 | c-----------------------------------------------------------------------
|
|---|
| 1164 | c
|
|---|
| 1165 | subroutine maph1_to_l2t(b,nb,a,na,if3d,w,ldw)
|
|---|
| 1166 | c
|
|---|
| 1167 | c Input: a
|
|---|
| 1168 | c Output: b
|
|---|
| 1169 | c
|
|---|
| 1170 | real a(1),b(1),w(1)
|
|---|
| 1171 | logical if3d
|
|---|
| 1172 | c
|
|---|
| 1173 | parameter(lx=50)
|
|---|
| 1174 | real za(lx),zb(lx)
|
|---|
| 1175 | c
|
|---|
| 1176 | real iba(lx*lx),ibat(lx*lx)
|
|---|
| 1177 | save iba,ibat
|
|---|
| 1178 | c
|
|---|
| 1179 | integer nao,nbo
|
|---|
| 1180 | save nao,nbo
|
|---|
| 1181 | data nao,nbo / -9, -9/
|
|---|
| 1182 | c
|
|---|
| 1183 | c
|
|---|
| 1184 | if (na.gt.lx.or.nb.gt.lx) then
|
|---|
| 1185 | write(6,*)'ERROR: increase lx in maph1_to_l2 to max:',na,nb
|
|---|
| 1186 | call exitt
|
|---|
| 1187 | endif
|
|---|
| 1188 | c
|
|---|
| 1189 | if (na.ne.nao .or. nb.ne.nbo) then
|
|---|
| 1190 | nao = na
|
|---|
| 1191 | nbo = nb
|
|---|
| 1192 | call zwgl (za,w,na)
|
|---|
| 1193 | call zwgll(zb,w,nb)
|
|---|
| 1194 | call igllm(iba,ibat,zb,za,nb,na,nb,na)
|
|---|
| 1195 | endif
|
|---|
| 1196 | c
|
|---|
| 1197 | call specmpn(b,nb,a,na,ibat,iba,if3d,w,ldw)
|
|---|
| 1198 | c
|
|---|
| 1199 | return
|
|---|
| 1200 | end
|
|---|
| 1201 | c
|
|---|
| 1202 | c-----------------------------------------------------------------------
|
|---|
| 1203 | subroutine irank_vec_tally(ind,nn,a,m,n,key,nkey,key2,aa)
|
|---|
| 1204 | c
|
|---|
| 1205 | c Compute rank of each unique entry a(1,i)
|
|---|
| 1206 | c
|
|---|
| 1207 | c Output: ind(i) i=1,...,n (global) rank of entry a(*,i)
|
|---|
| 1208 | c nn = max(rank)
|
|---|
| 1209 | c a(j,i) is destroyed
|
|---|
| 1210 | c a(1,i) tally of preceding structure values
|
|---|
| 1211 | c
|
|---|
| 1212 | c Input: a(j,i) j=1,...,m; i=1,...,n
|
|---|
| 1213 | c m : leading dim. of v (ldv must be .ge. m)
|
|---|
| 1214 | c key : sort key
|
|---|
| 1215 | c nkey :
|
|---|
| 1216 | c
|
|---|
| 1217 | c Although not mandatory, this ranking procedure is probably
|
|---|
| 1218 | c most effectively employed when the keys are pre-sorted. Thus,
|
|---|
| 1219 | c the option is provided to sort vi() prior to the ranking.
|
|---|
| 1220 | c
|
|---|
| 1221 | c
|
|---|
| 1222 | integer ind(n),a(m,n)
|
|---|
| 1223 | integer key(nkey),key2(0:3),aa(m)
|
|---|
| 1224 | logical iftuple_ianeb,a_ne_b
|
|---|
| 1225 | c
|
|---|
| 1226 | c
|
|---|
| 1227 | nk = min(nkey,m)
|
|---|
| 1228 | call ituple_sort(a,m,n,key,nk,ind,aa)
|
|---|
| 1229 | c do i=1,n
|
|---|
| 1230 | c write(6,*) i,' sort:',(a(k,i),k=1,3)
|
|---|
| 1231 | c enddo
|
|---|
| 1232 | c
|
|---|
| 1233 | c
|
|---|
| 1234 | c Find unique a's
|
|---|
| 1235 | c
|
|---|
| 1236 | call icopy(aa,a,m)
|
|---|
| 1237 | nn=1
|
|---|
| 1238 | mm=0
|
|---|
| 1239 | c
|
|---|
| 1240 | a(1,1) = nn
|
|---|
| 1241 | a(2,1)=ind(1)
|
|---|
| 1242 | a(3,1)=mm
|
|---|
| 1243 | c
|
|---|
| 1244 | do i=2,n
|
|---|
| 1245 | a_ne_b = iftuple_ianeb(aa,a(1,i),key,nk)
|
|---|
| 1246 | if (a_ne_b) then ! new structure
|
|---|
| 1247 | ms = aa(3) ! structure type
|
|---|
| 1248 | if (aa(2).eq.0) ms = aa(2) ! structure type
|
|---|
| 1249 | mm = mm+key2(ms) ! n dofs
|
|---|
| 1250 | call icopy(aa,a(1,i),m)
|
|---|
| 1251 | nn = nn+1
|
|---|
| 1252 | endif
|
|---|
| 1253 | a(1,i) = nn
|
|---|
| 1254 | a(2,i) = ind(i)
|
|---|
| 1255 | a(3,i) = mm
|
|---|
| 1256 | enddo
|
|---|
| 1257 | ms = aa(3)
|
|---|
| 1258 | if (aa(2).eq.0) ms = aa(2) ! structure type
|
|---|
| 1259 | nn = mm+key2(ms)
|
|---|
| 1260 | c
|
|---|
| 1261 | c Set ind() to rank
|
|---|
| 1262 | c
|
|---|
| 1263 | do i=1,n
|
|---|
| 1264 | iold=a(2,i)
|
|---|
| 1265 | ind(iold) = a(1,i)
|
|---|
| 1266 | enddo
|
|---|
| 1267 | c
|
|---|
| 1268 | c Set a1() to number of preceding dofs
|
|---|
| 1269 | c
|
|---|
| 1270 | do i=1,n
|
|---|
| 1271 | iold=a(2,i)
|
|---|
| 1272 | a(1,iold) = a(3,i)
|
|---|
| 1273 | enddo
|
|---|
| 1274 | c
|
|---|
| 1275 | return
|
|---|
| 1276 | end
|
|---|
| 1277 | c
|
|---|
| 1278 | c-----------------------------------------------------------------------
|
|---|
| 1279 | c
|
|---|
| 1280 | subroutine out_se1(se2crs,nx,name)
|
|---|
| 1281 | c
|
|---|
| 1282 | include 'SIZE'
|
|---|
| 1283 | integer se2crs(nx,nx,1)
|
|---|
| 1284 | character*4 name
|
|---|
| 1285 | c
|
|---|
| 1286 | write(6,*)
|
|---|
| 1287 | write(6,*) 'out_se',nx,name
|
|---|
| 1288 | do ie=nelv-1,1,-2
|
|---|
| 1289 | write(6,*)
|
|---|
| 1290 | do j=nx,1,-1
|
|---|
| 1291 | if(nx.eq.4) then
|
|---|
| 1292 | write(6,4) name,((se2crs(i,j,k+ie),i=1,nx),k=0,1)
|
|---|
| 1293 | elseif(nx.eq.3) then
|
|---|
| 1294 | write(6,3) name,((se2crs(i,j,k+ie),i=1,nx),k=0,1)
|
|---|
| 1295 | else
|
|---|
| 1296 | write(6,2) name,((se2crs(i,j,k+ie),i=1,nx),k=0,1)
|
|---|
| 1297 | endif
|
|---|
| 1298 | enddo
|
|---|
| 1299 | enddo
|
|---|
| 1300 | c
|
|---|
| 1301 | 4 format(a4,5x,2(4i5,3x))
|
|---|
| 1302 | 3 format(a4,5x,2(3i5,3x))
|
|---|
| 1303 | 2 format(a4,5x,2(2i5,3x))
|
|---|
| 1304 | c
|
|---|
| 1305 | return
|
|---|
| 1306 | end
|
|---|
| 1307 | c
|
|---|
| 1308 | c-----------------------------------------------------------------------
|
|---|
| 1309 | c
|
|---|
| 1310 | subroutine out_se0(se2crs,nx,nel,name)
|
|---|
| 1311 | c
|
|---|
| 1312 | include 'SIZE'
|
|---|
| 1313 | integer se2crs(nx,nx,1)
|
|---|
| 1314 | character*4 name
|
|---|
| 1315 | c
|
|---|
| 1316 | write(6,*)
|
|---|
| 1317 | write(6,*) 'out_se',nx,name,nel
|
|---|
| 1318 | do ie=nel-3,1,-4
|
|---|
| 1319 | write(6,*)
|
|---|
| 1320 | do j=nx,1,-1
|
|---|
| 1321 | if(nx.eq.4) then
|
|---|
| 1322 | write(6,4) name,((se2crs(i,j,k+ie),i=1,nx),k=0,3)
|
|---|
| 1323 | elseif(nx.eq.3) then
|
|---|
| 1324 | write(6,3) name,((se2crs(i,j,k+ie),i=1,nx),k=0,3)
|
|---|
| 1325 | else
|
|---|
| 1326 | write(6,2) name,((se2crs(i,j,k+ie),i=1,nx),k=0,3)
|
|---|
| 1327 | endif
|
|---|
| 1328 | enddo
|
|---|
| 1329 | enddo
|
|---|
| 1330 | c
|
|---|
| 1331 | 4 format(a4,5x,4(4i5,3x))
|
|---|
| 1332 | 3 format(a4,5x,4(3i5,3x))
|
|---|
| 1333 | 2 format(a4,5x,4(2i5,3x))
|
|---|
| 1334 | c
|
|---|
| 1335 | return
|
|---|
| 1336 | end
|
|---|
| 1337 | c
|
|---|
| 1338 | c-----------------------------------------------------------------------
|
|---|
| 1339 | c
|
|---|
| 1340 | subroutine crs_solve_h1(uf,vf)
|
|---|
| 1341 | c
|
|---|
| 1342 | c Given an input vector v, this generates the H1 coarse-grid solution
|
|---|
| 1343 | c
|
|---|
| 1344 | include 'SIZE'
|
|---|
| 1345 | include 'DOMAIN'
|
|---|
| 1346 | include 'INPUT'
|
|---|
| 1347 | include 'GEOM'
|
|---|
| 1348 | include 'SOLN'
|
|---|
| 1349 | include 'PARALLEL'
|
|---|
| 1350 | include 'CTIMER'
|
|---|
| 1351 | include 'TSTEP'
|
|---|
| 1352 |
|
|---|
| 1353 | real uf(1),vf(1)
|
|---|
| 1354 | common /scrpre/ uc(lcr*lelt)
|
|---|
| 1355 | common /scrpr2/ vc(lcr*lelt)
|
|---|
| 1356 | common /scrxxt/ cmlt(lcr,lelv),mask(lcr,lelv)
|
|---|
| 1357 |
|
|---|
| 1358 | integer icalld1
|
|---|
| 1359 | save icalld1
|
|---|
| 1360 | data icalld1 /0/
|
|---|
| 1361 |
|
|---|
| 1362 |
|
|---|
| 1363 | if (icalld1.eq.0) then ! timer info
|
|---|
| 1364 | ncrsl=0
|
|---|
| 1365 | tcrsl=0.0
|
|---|
| 1366 | icalld1=1
|
|---|
| 1367 | endif
|
|---|
| 1368 | ncrsl = ncrsl + 1
|
|---|
| 1369 |
|
|---|
| 1370 | ntot = nelv*lx1*ly1*lz1
|
|---|
| 1371 | call col3(uf,vf,vmult,ntot)
|
|---|
| 1372 |
|
|---|
| 1373 | call map_f_to_c_h1_bilin(vc,uf) ! additive Schwarz
|
|---|
| 1374 |
|
|---|
| 1375 | #ifdef TIMER
|
|---|
| 1376 | etime1=dnekclock()
|
|---|
| 1377 | #endif
|
|---|
| 1378 | call fgslib_crs_solve(xxth(ifield),uc,vc)
|
|---|
| 1379 | #ifdef TIMER
|
|---|
| 1380 | tcrsl=tcrsl+dnekclock()-etime1
|
|---|
| 1381 | #endif
|
|---|
| 1382 |
|
|---|
| 1383 | call map_c_to_f_h1_bilin(uf,uc)
|
|---|
| 1384 |
|
|---|
| 1385 |
|
|---|
| 1386 | return
|
|---|
| 1387 | end
|
|---|
| 1388 | c-----------------------------------------------------------------------
|
|---|
| 1389 | subroutine set_h1_basis_bilin
|
|---|
| 1390 | c
|
|---|
| 1391 | include 'SIZE'
|
|---|
| 1392 | include 'DOMAIN'
|
|---|
| 1393 | include 'WZ'
|
|---|
| 1394 | c
|
|---|
| 1395 | do ix=1,lx1
|
|---|
| 1396 | h1_basis(ix) = 0.5*(1.0-zgm1(ix,1))
|
|---|
| 1397 | h1_basis(ix+lx1) = 0.5*(1.0+zgm1(ix,1))
|
|---|
| 1398 | enddo
|
|---|
| 1399 | call transpose(h1_basist,2,h1_basis,lx1)
|
|---|
| 1400 | c
|
|---|
| 1401 | return
|
|---|
| 1402 | end
|
|---|
| 1403 | c
|
|---|
| 1404 | c-----------------------------------------------------------------------
|
|---|
| 1405 | c
|
|---|
| 1406 | subroutine map_c_to_f_h1_bilin(uf,uc)
|
|---|
| 1407 | c
|
|---|
| 1408 | c H1 Iterpolation operator: linear --> spectral GLL mesh
|
|---|
| 1409 | c
|
|---|
| 1410 | include 'SIZE'
|
|---|
| 1411 | include 'INPUT'
|
|---|
| 1412 | include 'DOMAIN'
|
|---|
| 1413 | c
|
|---|
| 1414 | parameter (lxyz = lx1*ly1*lz1)
|
|---|
| 1415 | real uc(2,2,ldim-1,lelt),uf(lxyz,lelt)
|
|---|
| 1416 | parameter (l2 = ldim-1)
|
|---|
| 1417 | common /ctmp0/ w(lx1,lx1,2),v(lx1,2,l2,lelt)
|
|---|
| 1418 | c
|
|---|
| 1419 | integer icalld
|
|---|
| 1420 | save icalld
|
|---|
| 1421 | data icalld/0/
|
|---|
| 1422 | if (icalld.eq.0) then
|
|---|
| 1423 | icalld=icalld+1
|
|---|
| 1424 | call set_h1_basis_bilin
|
|---|
| 1425 | endif
|
|---|
| 1426 | c
|
|---|
| 1427 | c
|
|---|
| 1428 | n2 = 2
|
|---|
| 1429 | if (if3d) then
|
|---|
| 1430 | c
|
|---|
| 1431 | n31 = n2*n2*nelv
|
|---|
| 1432 | n13 = lx1*lx1
|
|---|
| 1433 | c
|
|---|
| 1434 | call mxm(h1_basis,lx1,uc,n2,v,n31)
|
|---|
| 1435 | do ie=1,nelv
|
|---|
| 1436 | do iz=1,n2
|
|---|
| 1437 | call mxm(v(1,1,iz,ie),lx1,h1_basist,n2,w(1,1,iz),lx1)
|
|---|
| 1438 | enddo
|
|---|
| 1439 | call mxm(w,n13,h1_basist,n2,uf(1,ie),lx1)
|
|---|
| 1440 | enddo
|
|---|
| 1441 | c
|
|---|
| 1442 | else
|
|---|
| 1443 | c
|
|---|
| 1444 | n31 = 2*nelv
|
|---|
| 1445 | call mxm(h1_basis,lx1,uc,n2,v,n31)
|
|---|
| 1446 | do ie=1,nelv
|
|---|
| 1447 | call mxm(v(1,1,1,ie),lx1,h1_basist,n2,uf(1,ie),lx1)
|
|---|
| 1448 | enddo
|
|---|
| 1449 | endif
|
|---|
| 1450 | return
|
|---|
| 1451 | end
|
|---|
| 1452 | c
|
|---|
| 1453 | c-----------------------------------------------------------------------
|
|---|
| 1454 | c
|
|---|
| 1455 | subroutine map_f_to_c_h1_bilin(uc,uf)
|
|---|
| 1456 | c
|
|---|
| 1457 | c TRANSPOSE of H1 Iterpolation operator: T
|
|---|
| 1458 | c (linear --> spectral GLL mesh)
|
|---|
| 1459 | c
|
|---|
| 1460 | include 'SIZE'
|
|---|
| 1461 | include 'DOMAIN'
|
|---|
| 1462 | include 'INPUT'
|
|---|
| 1463 | c
|
|---|
| 1464 | parameter (lxyz = lx1*ly1*lz1)
|
|---|
| 1465 | real uc(lcr,lelt),uf(lx1,ly1,lz1,lelt)
|
|---|
| 1466 | common /ctmp0/ w(2,2,lx1),v(2,ly1,lz1,lelt)
|
|---|
| 1467 | c
|
|---|
| 1468 | integer icalld
|
|---|
| 1469 | save icalld
|
|---|
| 1470 | data icalld/0/
|
|---|
| 1471 | if (icalld.eq.0) then
|
|---|
| 1472 | icalld=icalld+1
|
|---|
| 1473 | call set_h1_basis_bilin
|
|---|
| 1474 | endif
|
|---|
| 1475 | c
|
|---|
| 1476 | n2 = 2
|
|---|
| 1477 | if (if3d) then
|
|---|
| 1478 | n31 = ly1*lz1*nelv
|
|---|
| 1479 | n13 = n2*n2
|
|---|
| 1480 | call mxm(h1_basist,n2,uf,lx1,v,n31)
|
|---|
| 1481 | do ie=1,nelv
|
|---|
| 1482 | do iz=1,lz1
|
|---|
| 1483 | call mxm(v(1,1,iz,ie),n2,h1_basis,lx1,w(1,1,iz),n2)
|
|---|
| 1484 | enddo
|
|---|
| 1485 | call mxm(w,n13,h1_basis,lx1,uc(1,ie),n2)
|
|---|
| 1486 | enddo
|
|---|
| 1487 | else
|
|---|
| 1488 | n31 = ly1*nelv
|
|---|
| 1489 | call mxm(h1_basist,n2,uf,lx1,v,n31)
|
|---|
| 1490 | do ie=1,nelv
|
|---|
| 1491 | call mxm(v(1,1,1,ie),n2,h1_basis,lx1,uc(1,ie),n2)
|
|---|
| 1492 | enddo
|
|---|
| 1493 | endif
|
|---|
| 1494 |
|
|---|
| 1495 | return
|
|---|
| 1496 | end
|
|---|
| 1497 | c-----------------------------------------------------------------------
|
|---|
| 1498 | subroutine get_local_crs_galerkin(a,ncl,nxc,h1,h2,w1,w2)
|
|---|
| 1499 |
|
|---|
| 1500 | c This routine generates Nelv submatrices of order ncl using
|
|---|
| 1501 | c Galerkin projection
|
|---|
| 1502 |
|
|---|
| 1503 | include 'SIZE'
|
|---|
| 1504 |
|
|---|
| 1505 | real a(ncl,ncl,1),h1(1),h2(1)
|
|---|
| 1506 | real w1(lx1*ly1*lz1,nelv),w2(lx1*ly1*lz1,nelv)
|
|---|
| 1507 |
|
|---|
| 1508 | parameter (lcrd=lx1**ldim)
|
|---|
| 1509 | common /ctmp1z/ b(lcrd,8)
|
|---|
| 1510 |
|
|---|
| 1511 | integer e
|
|---|
| 1512 |
|
|---|
| 1513 | do j=1,ncl
|
|---|
| 1514 | call gen_crs_basis(b(1,j),j) ! bi- or tri-linear interpolant
|
|---|
| 1515 | enddo
|
|---|
| 1516 |
|
|---|
| 1517 | isd = 1
|
|---|
| 1518 | imsh = 1
|
|---|
| 1519 |
|
|---|
| 1520 | nxyz = lx1*ly1*lz1
|
|---|
| 1521 | do j = 1,ncl
|
|---|
| 1522 | do e = 1,nelv
|
|---|
| 1523 | call copy(w1(1,e),b(1,j),nxyz)
|
|---|
| 1524 | enddo
|
|---|
| 1525 |
|
|---|
| 1526 | call axhelm (w2,w1,h1,h2,imsh,isd) ! A^e * bj
|
|---|
| 1527 |
|
|---|
| 1528 | do e = 1,nelv
|
|---|
| 1529 | do i = 1,ncl
|
|---|
| 1530 | a(i,j,e) = vlsc2(b(1,i),w2(1,e),nxyz) ! bi^T * A^e * bj
|
|---|
| 1531 | enddo
|
|---|
| 1532 | enddo
|
|---|
| 1533 |
|
|---|
| 1534 | enddo
|
|---|
| 1535 |
|
|---|
| 1536 | return
|
|---|
| 1537 | end
|
|---|
| 1538 | c-----------------------------------------------------------------------
|
|---|
| 1539 | subroutine gen_crs_basis(b,j) ! bi- tri-linear
|
|---|
| 1540 |
|
|---|
| 1541 | include 'SIZE'
|
|---|
| 1542 | real b(lx1,ly1,lz1)
|
|---|
| 1543 |
|
|---|
| 1544 | real z0(lx1),z1(lx1)
|
|---|
| 1545 | real zr(lx1),zs(lx1),zt(lx1)
|
|---|
| 1546 |
|
|---|
| 1547 | integer p,q,r
|
|---|
| 1548 |
|
|---|
| 1549 | call zwgll(zr,zs,lx1)
|
|---|
| 1550 |
|
|---|
| 1551 | do i=1,lx1
|
|---|
| 1552 | z0(i) = .5*(1-zr(i)) ! 1-->0
|
|---|
| 1553 | z1(i) = .5*(1+zr(i)) ! 0-->1
|
|---|
| 1554 | enddo
|
|---|
| 1555 |
|
|---|
| 1556 | call copy(zr,z0,lx1)
|
|---|
| 1557 | call copy(zs,z0,lx1)
|
|---|
| 1558 | call copy(zt,z0,lx1)
|
|---|
| 1559 |
|
|---|
| 1560 | if (mod(j,2).eq.0) call copy(zr,z1,lx1)
|
|---|
| 1561 | if (j.eq.3.or.j.eq.4.or.j.eq.7.or.j.eq.8) call copy(zs,z1,lx1)
|
|---|
| 1562 | if (j.gt.4) call copy(zt,z1,lx1)
|
|---|
| 1563 |
|
|---|
| 1564 | if (ldim.eq.3) then
|
|---|
| 1565 | do r=1,lx1
|
|---|
| 1566 | do q=1,lx1
|
|---|
| 1567 | do p=1,lx1
|
|---|
| 1568 | b(p,q,r) = zr(p)*zs(q)*zt(r)
|
|---|
| 1569 | enddo
|
|---|
| 1570 | enddo
|
|---|
| 1571 | enddo
|
|---|
| 1572 | else
|
|---|
| 1573 | do q=1,lx1
|
|---|
| 1574 | do p=1,lx1
|
|---|
| 1575 | b(p,q,1) = zr(p)*zs(q)
|
|---|
| 1576 | enddo
|
|---|
| 1577 | enddo
|
|---|
| 1578 | endif
|
|---|
| 1579 |
|
|---|
| 1580 | return
|
|---|
| 1581 | end
|
|---|
| 1582 | c-----------------------------------------------------------------------
|
|---|
| 1583 | subroutine gen_crs_basis2(b,j) ! bi- tri-quadratic
|
|---|
| 1584 |
|
|---|
| 1585 | include 'SIZE'
|
|---|
| 1586 | real b(lx1,ly1,lz1)
|
|---|
| 1587 |
|
|---|
| 1588 | real z0(lx1),z1(lx1),z2(lx1)
|
|---|
| 1589 | real zr(lx1),zs(lx1),zt(lx1)
|
|---|
| 1590 |
|
|---|
| 1591 | integer p,q,r
|
|---|
| 1592 |
|
|---|
| 1593 | call zwgll(zr,zs,lx1)
|
|---|
| 1594 |
|
|---|
| 1595 | do i=1,lx1
|
|---|
| 1596 | z0(i) = .5*(zr(i)-1)*zr(i) ! 1-->0 ! Lagrangian, ordered
|
|---|
| 1597 | z1(i) = 4.*(1+zr(i))*(1-zr(i)) ! lexicographically
|
|---|
| 1598 | z2(i) = .5*(zr(i)+1)*zr(i) ! 0-->1 !
|
|---|
| 1599 | enddo
|
|---|
| 1600 |
|
|---|
| 1601 | call copy(zr,z0,lx1)
|
|---|
| 1602 | call copy(zs,z0,lx1)
|
|---|
| 1603 | call copy(zt,z0,lx1)
|
|---|
| 1604 |
|
|---|
| 1605 | if (mod(j,2).eq.0) call copy(zr,z1,lx1)
|
|---|
| 1606 | if (j.eq.3.or.j.eq.4.or.j.eq.7.or.j.eq.8) call copy(zs,z1,lx1)
|
|---|
| 1607 | if (j.gt.4) call copy(zt,z1,lx1)
|
|---|
| 1608 |
|
|---|
| 1609 | if (ldim.eq.3) then
|
|---|
| 1610 | do r=1,lx1
|
|---|
| 1611 | do q=1,lx1
|
|---|
| 1612 | do p=1,lx1
|
|---|
| 1613 | b(p,q,r) = zr(p)*zs(q)*zt(r)
|
|---|
| 1614 | enddo
|
|---|
| 1615 | enddo
|
|---|
| 1616 | enddo
|
|---|
| 1617 | else
|
|---|
| 1618 | do q=1,lx1
|
|---|
| 1619 | do p=1,lx1
|
|---|
| 1620 | b(p,q,1) = zr(p)*zs(q)
|
|---|
| 1621 | enddo
|
|---|
| 1622 | enddo
|
|---|
| 1623 | endif
|
|---|
| 1624 |
|
|---|
| 1625 | return
|
|---|
| 1626 | end
|
|---|
| 1627 | c-----------------------------------------------------------------------
|
|---|
| 1628 | subroutine get_vertex
|
|---|
| 1629 | include 'SIZE'
|
|---|
| 1630 | include 'TOTAL'
|
|---|
| 1631 |
|
|---|
| 1632 | common /ivrtx/ vertex ((2**ldim)*lelt)
|
|---|
| 1633 | integer vertex
|
|---|
| 1634 |
|
|---|
| 1635 | call get_vert
|
|---|
| 1636 |
|
|---|
| 1637 | return
|
|---|
| 1638 | end
|
|---|
| 1639 | c-----------------------------------------------------------------------
|
|---|
| 1640 | c-----------------------------------------------------------------------
|
|---|
| 1641 | subroutine irank_vecn(ind,nn,a,m,n,key,nkey,aa)
|
|---|
| 1642 | c
|
|---|
| 1643 | c Compute rank of each unique entry a(1,i)
|
|---|
| 1644 | c
|
|---|
| 1645 | c Output: ind(i) i=1,...,n (global) rank of entry a(*,i)
|
|---|
| 1646 | c nn = max(rank)
|
|---|
| 1647 | c a(j,i) is permuted
|
|---|
| 1648 | c
|
|---|
| 1649 | c Input: a(j,i) j=1,...,m; i=1,...,n
|
|---|
| 1650 | c m : leading dim. of v (ldv must be .ge. m)
|
|---|
| 1651 | c key : sort key
|
|---|
| 1652 | c nkey :
|
|---|
| 1653 | c
|
|---|
| 1654 | c Although not mandatory, this ranking procedure is probably
|
|---|
| 1655 | c most effectively employed when the keys are pre-sorted. Thus,
|
|---|
| 1656 | c the option is provided to sort vi() prior to the ranking.
|
|---|
| 1657 | c
|
|---|
| 1658 | c
|
|---|
| 1659 | integer ind(n),a(m,n)
|
|---|
| 1660 | integer key(nkey),aa(m)
|
|---|
| 1661 | logical iftuple_ianeb,a_ne_b
|
|---|
| 1662 |
|
|---|
| 1663 | nk = min(nkey,m)
|
|---|
| 1664 | call ituple_sort(a,m,n,key,nk,ind,aa)
|
|---|
| 1665 |
|
|---|
| 1666 | c Find unique a's
|
|---|
| 1667 | call icopy(aa,a,m)
|
|---|
| 1668 | nn = 1
|
|---|
| 1669 | ind(1) = nn
|
|---|
| 1670 | c
|
|---|
| 1671 | do i=2,n
|
|---|
| 1672 | a_ne_b = iftuple_ianeb(aa,a(1,i),key,nk)
|
|---|
| 1673 | if (a_ne_b) then
|
|---|
| 1674 | call icopy(aa,a(1,i),m)
|
|---|
| 1675 | nn = nn+1
|
|---|
| 1676 | endif
|
|---|
| 1677 | ind(i) = nn ! set ind() to rank
|
|---|
| 1678 | enddo
|
|---|
| 1679 |
|
|---|
| 1680 | return
|
|---|
| 1681 | end
|
|---|
| 1682 | c-----------------------------------------------------------------------
|
|---|
| 1683 | subroutine gbtuple_rank(tuple,m,n,nmax,cr_h,nid,np,ind)
|
|---|
| 1684 | c
|
|---|
| 1685 | c Return a unique rank for each matched tuple set. Global. Balanced.
|
|---|
| 1686 | c
|
|---|
| 1687 | c tuple is destroyed.
|
|---|
| 1688 | c
|
|---|
| 1689 | c By "balanced" we mean that none of the tuple entries is likely to
|
|---|
| 1690 | c be much more uniquely populated than any other, so that any of
|
|---|
| 1691 | c the tuples can serve as an initial (parallel) sort key
|
|---|
| 1692 | c
|
|---|
| 1693 | c First two slots in tuple(:,i) assumed empty
|
|---|
| 1694 | c
|
|---|
| 1695 | integer ind(nmax),tuple(m,nmax),cr_h
|
|---|
| 1696 |
|
|---|
| 1697 | parameter (mmax=40)
|
|---|
| 1698 | integer key(mmax),wtuple(mmax)
|
|---|
| 1699 |
|
|---|
| 1700 | if (m.gt.mmax) then
|
|---|
| 1701 | write(6,*) nid,m,mmax,' gbtuple_rank fail'
|
|---|
| 1702 | call exitt
|
|---|
| 1703 | endif
|
|---|
| 1704 |
|
|---|
| 1705 | do i=1,n
|
|---|
| 1706 | tuple(1,i) = mod(tuple(3,i),np) ! destination processor
|
|---|
| 1707 | tuple(2,i) = i ! return location
|
|---|
| 1708 | enddo
|
|---|
| 1709 |
|
|---|
| 1710 | ni= n
|
|---|
| 1711 | ky=1 ! Assumes crystal_new already called
|
|---|
| 1712 | call fgslib_crystal_ituple_transfer(cr_h, tuple,m,ni,nmax, ky)
|
|---|
| 1713 |
|
|---|
| 1714 | nimx = iglmax(ni,1)
|
|---|
| 1715 | if (ni.gt.nmax) write(6,*) ni,nmax,n,'cr_xfer problem, A'
|
|---|
| 1716 | if (nimx.gt.nmax) call exitt
|
|---|
| 1717 |
|
|---|
| 1718 | nkey = m-2
|
|---|
| 1719 | do k=1,nkey
|
|---|
| 1720 | key(k) = k+2
|
|---|
| 1721 | enddo
|
|---|
| 1722 |
|
|---|
| 1723 | call irank_vecn(ind,nu,tuple,m,ni,key,nkey,wtuple)! tuple re-ordered,
|
|---|
| 1724 | ! but contents same
|
|---|
| 1725 |
|
|---|
| 1726 | nu_tot = igl_running_sum(nu) ! running sum over P processors
|
|---|
| 1727 | nu_prior = nu_tot - nu
|
|---|
| 1728 |
|
|---|
| 1729 | do i=1,ni
|
|---|
| 1730 | tuple(3,i) = ind(i) + nu_prior ! global ranking
|
|---|
| 1731 | enddo
|
|---|
| 1732 |
|
|---|
| 1733 | call fgslib_crystal_ituple_transfer(cr_h, tuple,m,ni,nmax, ky)
|
|---|
| 1734 |
|
|---|
| 1735 | nk = 1 ! restore to original order, local rank: 2; global: 3
|
|---|
| 1736 | ky = 2
|
|---|
| 1737 | call ituple_sort(tuple,m,n,ky,nk,ind,wtuple)
|
|---|
| 1738 |
|
|---|
| 1739 |
|
|---|
| 1740 | return
|
|---|
| 1741 | end
|
|---|
| 1742 | c-----------------------------------------------------------------------
|
|---|
| 1743 | subroutine setvert3d(glo_num,ngv,nx,nel,vertex,ifcenter)
|
|---|
| 1744 | c
|
|---|
| 1745 | c setup unique ids for dssum
|
|---|
| 1746 | c note:
|
|---|
| 1747 | c total number of unique vertices, edges and faces has to be smaller
|
|---|
| 1748 | c than 2**31 (integer-4 limit).
|
|---|
| 1749 | c if nelgt < 2**31/12 we're ok for sure (independent of N)!
|
|---|
| 1750 | c
|
|---|
| 1751 | include 'SIZE'
|
|---|
| 1752 | include 'CTIMER'
|
|---|
| 1753 | include 'PARALLEL'
|
|---|
| 1754 | include 'TOPOL'
|
|---|
| 1755 | include 'GEOM'
|
|---|
| 1756 |
|
|---|
| 1757 | integer*8 glo_num(1),ngv
|
|---|
| 1758 | integer vertex(0:1,0:1,0:1,1),nx
|
|---|
| 1759 | logical ifcenter
|
|---|
| 1760 |
|
|---|
| 1761 | integer edge(0:1,0:1,0:1,3,lelt),enum(12,lelt),fnum(6,lelt)
|
|---|
| 1762 | common /scrmg/ edge,enum,fnum
|
|---|
| 1763 |
|
|---|
| 1764 | parameter (nsafe=8) ! OFTEN, nsafe=2 suffices
|
|---|
| 1765 | integer etuple(4,12*lelt*nsafe),ftuple(5,6,lelt*nsafe)
|
|---|
| 1766 | integer ind(12*lelt*nsafe)
|
|---|
| 1767 | common /scrns/ ind,etuple
|
|---|
| 1768 | equivalence (etuple,ftuple)
|
|---|
| 1769 |
|
|---|
| 1770 | integer gvf(4),facet(4),aa(3),key(3),e
|
|---|
| 1771 | logical ifij
|
|---|
| 1772 |
|
|---|
| 1773 | integer*8 igv,ig0
|
|---|
| 1774 | integer*8 ngvv,ngve,ngvs,ngvi,ngvm
|
|---|
| 1775 | integer*8 n_on_edge,n_on_face,n_in_interior
|
|---|
| 1776 | integer*8 i8glmax
|
|---|
| 1777 | c
|
|---|
| 1778 | ny = nx
|
|---|
| 1779 | nz = nx
|
|---|
| 1780 | nxyz = nx*ny*nz
|
|---|
| 1781 | c
|
|---|
| 1782 | key(1)=1
|
|---|
| 1783 | key(2)=2
|
|---|
| 1784 | key(3)=3
|
|---|
| 1785 | c
|
|---|
| 1786 | c Assign hypercube ordering of vertices
|
|---|
| 1787 | c -------------------------------------
|
|---|
| 1788 | c
|
|---|
| 1789 | c Count number of unique vertices
|
|---|
| 1790 | nlv = 2**ldim
|
|---|
| 1791 | ngvv = iglmax(vertex,nlv*nel)
|
|---|
| 1792 | c
|
|---|
| 1793 | do e=1,nel
|
|---|
| 1794 | do k=0,1
|
|---|
| 1795 | do j=0,1
|
|---|
| 1796 | do i=0,1
|
|---|
| 1797 | c Local to global node number (vertex)
|
|---|
| 1798 | il = 1 + (nx-1)*i + nx*(nx-1)*j + nx*nx*(nx-1)*k
|
|---|
| 1799 | ile = il + nx*ny*nz*(e-1)
|
|---|
| 1800 | glo_num(ile) = vertex(i,j,k,e)
|
|---|
| 1801 | enddo
|
|---|
| 1802 | enddo
|
|---|
| 1803 | enddo
|
|---|
| 1804 | enddo
|
|---|
| 1805 | ngv = ngvv
|
|---|
| 1806 | c
|
|---|
| 1807 | if (nx.eq.2) return
|
|---|
| 1808 | c
|
|---|
| 1809 | c Assign global vertex numbers to SEM nodes on each edge
|
|---|
| 1810 | c ------------------------------------------------------
|
|---|
| 1811 | c
|
|---|
| 1812 | c Assign edge labels by bounding vertices.
|
|---|
| 1813 | do e=1,nel
|
|---|
| 1814 | do k=0,1
|
|---|
| 1815 | do j=0,1
|
|---|
| 1816 | do i=0,1
|
|---|
| 1817 | edge(i,j,k,1,e) = vertex(i,j,k,e) ! r-edge
|
|---|
| 1818 | edge(j,i,k,2,e) = vertex(i,j,k,e) ! s-edge
|
|---|
| 1819 | edge(k,i,j,3,e) = vertex(i,j,k,e) ! t-edge
|
|---|
| 1820 | enddo
|
|---|
| 1821 | enddo
|
|---|
| 1822 | enddo
|
|---|
| 1823 | enddo
|
|---|
| 1824 | c
|
|---|
| 1825 | c Sort edges by bounding vertices.
|
|---|
| 1826 | do i=0,12*nel-1
|
|---|
| 1827 | if (edge(0,i,0,1,1).gt.edge(1,i,0,1,1)) then
|
|---|
| 1828 | kswap = edge(0,i,0,1,1)
|
|---|
| 1829 | edge(0,i,0,1,1) = edge(1,i,0,1,1)
|
|---|
| 1830 | edge(1,i,0,1,1) = kswap
|
|---|
| 1831 | endif
|
|---|
| 1832 | etuple(3,i+1) = edge(0,i,0,1,1)
|
|---|
| 1833 | etuple(4,i+1) = edge(1,i,0,1,1)
|
|---|
| 1834 | enddo
|
|---|
| 1835 | c
|
|---|
| 1836 | c Assign a number (rank) to each unique edge
|
|---|
| 1837 | m = 4
|
|---|
| 1838 | n = 12*nel
|
|---|
| 1839 | nmax = 12*lelt*nsafe ! nsafe for crystal router factor of safety
|
|---|
| 1840 | call gbtuple_rank(etuple,m,n,nmax,cr_h,nid,np,ind)
|
|---|
| 1841 | do i=1,12*nel
|
|---|
| 1842 | enum(i,1) = etuple(3,i)
|
|---|
| 1843 | enddo
|
|---|
| 1844 | n_unique_edges = iglmax(enum,12*nel)
|
|---|
| 1845 | c
|
|---|
| 1846 | n_on_edge = nx-2
|
|---|
| 1847 | ngve = n_unique_edges*n_on_edge
|
|---|
| 1848 | do e=1,nel
|
|---|
| 1849 | iedg_loc = 0
|
|---|
| 1850 | c
|
|---|
| 1851 | c Edges 1-4
|
|---|
| 1852 | do k=0,1
|
|---|
| 1853 | do j=0,1
|
|---|
| 1854 | igv = ngv + n_on_edge*(enum(iedg_loc+1,e)-1)
|
|---|
| 1855 | i0 = nx*(nx-1)*j + nx*nx*(nx-1)*k
|
|---|
| 1856 | i0e = i0 + nxyz*(e-1)
|
|---|
| 1857 | if (glo_num(i0e+1).lt.glo_num(i0e+nx)) then
|
|---|
| 1858 | do i=2,nx-1 ! std forward case
|
|---|
| 1859 | glo_num(i0e+i) = igv + i-1
|
|---|
| 1860 | enddo
|
|---|
| 1861 | else
|
|---|
| 1862 | do i=2,nx-1 ! backward case
|
|---|
| 1863 | glo_num(i0e+i) = igv + 1 + n_on_edge-(i-1)
|
|---|
| 1864 | enddo
|
|---|
| 1865 | endif
|
|---|
| 1866 | iedg_loc = iedg_loc + 1
|
|---|
| 1867 | enddo
|
|---|
| 1868 | enddo
|
|---|
| 1869 | c
|
|---|
| 1870 | c Edges 5-8
|
|---|
| 1871 | do k=0,1
|
|---|
| 1872 | do i=0,1
|
|---|
| 1873 | igv = ngv + n_on_edge*(enum(iedg_loc+1,e)-1)
|
|---|
| 1874 | i0 = 1+(nx-1)*i + nx*nx*(nx-1)*k
|
|---|
| 1875 | i0e = i0 + nxyz*(e-1)
|
|---|
| 1876 | if (glo_num(i0e).lt.glo_num(i0e+nx*(nx-1))) then
|
|---|
| 1877 | do j=2,nx-1 ! std forward case
|
|---|
| 1878 | glo_num(i0e+(j-1)*nx) = igv + j-1
|
|---|
| 1879 | enddo
|
|---|
| 1880 | else
|
|---|
| 1881 | do j=2,nx-1 ! backward case
|
|---|
| 1882 | glo_num(i0e+(j-1)*nx) = igv + 1 + n_on_edge-(j-1)
|
|---|
| 1883 | enddo
|
|---|
| 1884 | endif
|
|---|
| 1885 | iedg_loc = iedg_loc + 1
|
|---|
| 1886 | enddo
|
|---|
| 1887 | enddo
|
|---|
| 1888 | c
|
|---|
| 1889 | c Edges 9-12
|
|---|
| 1890 | do j=0,1
|
|---|
| 1891 | do i=0,1
|
|---|
| 1892 | igv = ngv + n_on_edge*(enum(iedg_loc+1,e)-1)
|
|---|
| 1893 | i0 = 1 + (nx-1)*i + nx*(nx-1)*j
|
|---|
| 1894 | i0e = i0 + nxyz*(e-1)
|
|---|
| 1895 | if (glo_num(i0e).lt.glo_num(i0e+nx*nx*(nx-1))) then
|
|---|
| 1896 | do k=2,nx-1 ! std forward case
|
|---|
| 1897 | glo_num(i0e+(k-1)*nx*nx) = igv + k-1
|
|---|
| 1898 | enddo
|
|---|
| 1899 | else
|
|---|
| 1900 | do k=2,nx-1 ! backward case
|
|---|
| 1901 | glo_num(i0e+(k-1)*nx*nx) = igv + 1 + n_on_edge-(k-1)
|
|---|
| 1902 | enddo
|
|---|
| 1903 | endif
|
|---|
| 1904 | iedg_loc = iedg_loc + 1
|
|---|
| 1905 | enddo
|
|---|
| 1906 | enddo
|
|---|
| 1907 | enddo
|
|---|
| 1908 | ngv = ngv + ngve
|
|---|
| 1909 | c
|
|---|
| 1910 | c Asign global node numbers on the interior of each face
|
|---|
| 1911 | c ------------------------------------------------------
|
|---|
| 1912 | c
|
|---|
| 1913 | c Assign faces by 3-tuples
|
|---|
| 1914 | c
|
|---|
| 1915 | c (The following variables all take the symmetric
|
|---|
| 1916 | c notation of IFACE as arguments:)
|
|---|
| 1917 | c
|
|---|
| 1918 | c ICFACE(i,IFACE) - Gives the 4 vertices which reside on face IFACE
|
|---|
| 1919 | c as depicted below, e.g. ICFACE(i,2)=2,4,6,8.
|
|---|
| 1920 | c
|
|---|
| 1921 | c 3+-----+4 ^ Y
|
|---|
| 1922 | c / 2 /| |
|
|---|
| 1923 | c Edge 1 extends / / | |
|
|---|
| 1924 | c from vertex 7+-----+8 +2 +----> X
|
|---|
| 1925 | c 1 to 2. | 4 | / /
|
|---|
| 1926 | c | |/ /
|
|---|
| 1927 | c 5+-----+6 Z
|
|---|
| 1928 | c 3
|
|---|
| 1929 | c
|
|---|
| 1930 | nfaces=ldim*2
|
|---|
| 1931 | ncrnr =2**(ldim-1)
|
|---|
| 1932 | do e=1,nel
|
|---|
| 1933 | do ifac=1,nfaces
|
|---|
| 1934 | do icrn=1,ncrnr
|
|---|
| 1935 | i = icface(icrn,ifac)-1
|
|---|
| 1936 | facet(icrn) = vertex(i,0,0,e)
|
|---|
| 1937 | enddo
|
|---|
| 1938 | call isort(facet,ind,ncrnr)
|
|---|
| 1939 | call icopy(ftuple(3,ifac,e),facet,ncrnr-1)
|
|---|
| 1940 | enddo
|
|---|
| 1941 | enddo
|
|---|
| 1942 |
|
|---|
| 1943 | c Assign a number (rank) to each unique face
|
|---|
| 1944 | m = 5
|
|---|
| 1945 | n = 6*nel
|
|---|
| 1946 | nmax = 6*lelt*nsafe ! nsafe for crystal router factor of safety
|
|---|
| 1947 | call gbtuple_rank(ftuple,m,n,nmax,cr_h,nid,np,ind)
|
|---|
| 1948 | do i=1,6*nel
|
|---|
| 1949 | fnum(i,1) = ftuple(3,i,1)
|
|---|
| 1950 | enddo
|
|---|
| 1951 | n_unique_faces = iglmax(fnum,6*nel)
|
|---|
| 1952 | c
|
|---|
| 1953 | call dsset (nx,ny,nz)
|
|---|
| 1954 | do e=1,nel
|
|---|
| 1955 | do iface=1,nfaces
|
|---|
| 1956 | i0 = skpdat(1,iface)
|
|---|
| 1957 | i1 = skpdat(2,iface)
|
|---|
| 1958 | is = skpdat(3,iface)
|
|---|
| 1959 | j0 = skpdat(4,iface)
|
|---|
| 1960 | j1 = skpdat(5,iface)
|
|---|
| 1961 | js = skpdat(6,iface)
|
|---|
| 1962 | c
|
|---|
| 1963 | c On each face, count from minimum global vertex number,
|
|---|
| 1964 | c towards smallest adjacent vertex number. e.g., suppose
|
|---|
| 1965 | c the face is defined by the following global vertex numbers:
|
|---|
| 1966 | c
|
|---|
| 1967 | c
|
|---|
| 1968 | c 11+--------+81
|
|---|
| 1969 | c |c d|
|
|---|
| 1970 | c | |
|
|---|
| 1971 | c | |
|
|---|
| 1972 | c |a b|
|
|---|
| 1973 | c 15+--------+62
|
|---|
| 1974 | c
|
|---|
| 1975 | c We would count from c-->a, then towards d.
|
|---|
| 1976 | c
|
|---|
| 1977 | gvf(1) = glo_num(i0+nx*(j0-1)+nxyz*(e-1))
|
|---|
| 1978 | gvf(2) = glo_num(i1+nx*(j0-1)+nxyz*(e-1))
|
|---|
| 1979 | gvf(3) = glo_num(i0+nx*(j1-1)+nxyz*(e-1))
|
|---|
| 1980 | gvf(4) = glo_num(i1+nx*(j1-1)+nxyz*(e-1))
|
|---|
| 1981 | c
|
|---|
| 1982 | call irank(gvf,ind,4)
|
|---|
| 1983 | c
|
|---|
| 1984 | c ind(1) tells which element of gvf() is smallest.
|
|---|
| 1985 | c
|
|---|
| 1986 | ifij = .false.
|
|---|
| 1987 | if (ind(1).eq.1) then
|
|---|
| 1988 | idir = 1
|
|---|
| 1989 | jdir = 1
|
|---|
| 1990 | if (gvf(2).lt.gvf(3)) ifij = .true.
|
|---|
| 1991 | elseif (ind(1).eq.2) then
|
|---|
| 1992 | idir = -1
|
|---|
| 1993 | jdir = 1
|
|---|
| 1994 | if (gvf(1).lt.gvf(4)) ifij = .true.
|
|---|
| 1995 | elseif (ind(1).eq.3) then
|
|---|
| 1996 | idir = 1
|
|---|
| 1997 | jdir = -1
|
|---|
| 1998 | if (gvf(4).lt.gvf(1)) ifij = .true.
|
|---|
| 1999 | elseif (ind(1).eq.4) then
|
|---|
| 2000 | idir = -1
|
|---|
| 2001 | jdir = -1
|
|---|
| 2002 | if (gvf(3).lt.gvf(2)) ifij = .true.
|
|---|
| 2003 | endif
|
|---|
| 2004 | c
|
|---|
| 2005 | if (idir.lt.0) then
|
|---|
| 2006 | it=i0
|
|---|
| 2007 | i0=i1
|
|---|
| 2008 | i1=it
|
|---|
| 2009 | is=-is
|
|---|
| 2010 | endif
|
|---|
| 2011 | c
|
|---|
| 2012 | if (jdir.lt.0) then
|
|---|
| 2013 | jt=j0
|
|---|
| 2014 | j0=j1
|
|---|
| 2015 | j1=jt
|
|---|
| 2016 | js=-js
|
|---|
| 2017 | endif
|
|---|
| 2018 | c
|
|---|
| 2019 | nxx = nx*nx
|
|---|
| 2020 | n_on_face = (nx-2)*(ny-2)
|
|---|
| 2021 | ngvs = n_unique_faces*n_on_face
|
|---|
| 2022 | ig0 = ngv + n_on_face*(fnum(iface,e)-1)
|
|---|
| 2023 | if (ifij) then
|
|---|
| 2024 | k=0
|
|---|
| 2025 | l=0
|
|---|
| 2026 | do j=j0,j1,js
|
|---|
| 2027 | do i=i0,i1,is
|
|---|
| 2028 | k=k+1
|
|---|
| 2029 | c this is a serious kludge to stay on the face interior
|
|---|
| 2030 | if (k.gt.nx.and.k.lt.nxx-nx .and.
|
|---|
| 2031 | $ mod(k,nx).ne.1.and.mod(k,nx).ne.0) then
|
|---|
| 2032 | c interior
|
|---|
| 2033 | l = l+1
|
|---|
| 2034 | glo_num(i+nx*(j-1)+nxyz*(e-1)) = l + ig0
|
|---|
| 2035 | endif
|
|---|
| 2036 | enddo
|
|---|
| 2037 | enddo
|
|---|
| 2038 | else
|
|---|
| 2039 | k=0
|
|---|
| 2040 | l=0
|
|---|
| 2041 | do i=i0,i1,is
|
|---|
| 2042 | do j=j0,j1,js
|
|---|
| 2043 | k=k+1
|
|---|
| 2044 | c this is a serious kludge to stay on the face interior
|
|---|
| 2045 | if (k.gt.nx.and.k.lt.nxx-nx .and.
|
|---|
| 2046 | $ mod(k,nx).ne.1.and.mod(k,nx).ne.0) then
|
|---|
| 2047 | c interior
|
|---|
| 2048 | l = l+1
|
|---|
| 2049 | glo_num(i+nx*(j-1)+nxyz*(e-1)) = l + ig0
|
|---|
| 2050 | endif
|
|---|
| 2051 | enddo
|
|---|
| 2052 | enddo
|
|---|
| 2053 | endif
|
|---|
| 2054 | enddo
|
|---|
| 2055 | enddo
|
|---|
| 2056 | ngv = ngv + ngvs
|
|---|
| 2057 | c
|
|---|
| 2058 | c Finally, number interiors (only ifcenter=.true.)
|
|---|
| 2059 | c -------------------------------------------------
|
|---|
| 2060 | c
|
|---|
| 2061 | n_in_interior = (nx-2)*(ny-2)*(nz-2)
|
|---|
| 2062 | ngvi = n_in_interior*nelgt
|
|---|
| 2063 | if (ifcenter) then
|
|---|
| 2064 | do e=1,nel
|
|---|
| 2065 | ig0 = ngv + n_in_interior*(lglel(e)-1)
|
|---|
| 2066 | l = 0
|
|---|
| 2067 | do k=2,nz-1
|
|---|
| 2068 | do j=2,ny-1
|
|---|
| 2069 | do i=2,nx-1
|
|---|
| 2070 | l = l+1
|
|---|
| 2071 | glo_num(i+nx*(j-1)+nx*ny*(k-1)+nxyz*(e-1)) = ig0+l
|
|---|
| 2072 | enddo
|
|---|
| 2073 | enddo
|
|---|
| 2074 | enddo
|
|---|
| 2075 | enddo
|
|---|
| 2076 | ngv = ngv + ngvi
|
|---|
| 2077 | else
|
|---|
| 2078 | do e=1,nel
|
|---|
| 2079 | l = 0
|
|---|
| 2080 | do k=2,nz-1
|
|---|
| 2081 | do j=2,ny-1
|
|---|
| 2082 | do i=2,nx-1
|
|---|
| 2083 | l = l+1
|
|---|
| 2084 | glo_num(i+nx*(j-1)+nx*ny*(k-1)+nxyz*(e-1)) = 0
|
|---|
| 2085 | enddo
|
|---|
| 2086 | enddo
|
|---|
| 2087 | enddo
|
|---|
| 2088 | enddo
|
|---|
| 2089 | endif
|
|---|
| 2090 | c
|
|---|
| 2091 | c Quick check on maximum #dofs:
|
|---|
| 2092 | m = nxyz*nelt
|
|---|
| 2093 | ngvm = i8glmax(glo_num,m)
|
|---|
| 2094 | ngvv = ngvv + ngve + ngvs ! number of unique ids w/o interior
|
|---|
| 2095 | ngvi = ngvi + ngvv ! total number of unique ids
|
|---|
| 2096 | if (nio.eq.0) write(6,1) nx,ngvv,ngvi,ngv,ngvm
|
|---|
| 2097 | 1 format(' setvert3d:',i4,4i12)
|
|---|
| 2098 | c
|
|---|
| 2099 | return
|
|---|
| 2100 | end
|
|---|
| 2101 | c-----------------------------------------------------------------------
|
|---|
| 2102 | subroutine setvert2d(glo_num,ngv,nx,nel,vertex,ifcenter)
|
|---|
| 2103 | c
|
|---|
| 2104 | c setup unique ids for dssum
|
|---|
| 2105 | c
|
|---|
| 2106 | include 'SIZE'
|
|---|
| 2107 | include 'CTIMER'
|
|---|
| 2108 | include 'PARALLEL'
|
|---|
| 2109 | include 'TOPOL'
|
|---|
| 2110 | include 'GEOM'
|
|---|
| 2111 |
|
|---|
| 2112 | integer*8 glo_num(1),ngv
|
|---|
| 2113 | integer vertex(0:1,0:1,1),nx
|
|---|
| 2114 | logical ifcenter
|
|---|
| 2115 |
|
|---|
| 2116 | integer edge(0:1,0:1,2,lelt),enum(4,lelt)
|
|---|
| 2117 | common /scrmg/ edge,enum
|
|---|
| 2118 |
|
|---|
| 2119 | parameter (nsafe=8) ! OFTEN, nsafe=2 suffices
|
|---|
| 2120 | integer etuple(4,4*lelt*nsafe),ind(4*lelt*nsafe)
|
|---|
| 2121 | common /scrns/ ind,etuple
|
|---|
| 2122 |
|
|---|
| 2123 | integer gvf(4),aa(3),key(3),e,eg
|
|---|
| 2124 | logical ifij
|
|---|
| 2125 |
|
|---|
| 2126 | integer*8 igv,ig0
|
|---|
| 2127 | integer*8 ngvv,ngve,ngvs,ngvi,ngvm
|
|---|
| 2128 | integer*8 n_on_edge,n_on_face,n_in_interior
|
|---|
| 2129 | integer*8 i8glmax
|
|---|
| 2130 | c
|
|---|
| 2131 | c
|
|---|
| 2132 | c memory check...
|
|---|
| 2133 | c
|
|---|
| 2134 | ny = nx
|
|---|
| 2135 | nz = 1
|
|---|
| 2136 | nxyz = nx*ny*nz
|
|---|
| 2137 | c
|
|---|
| 2138 | key(1)=1
|
|---|
| 2139 | key(2)=2
|
|---|
| 2140 | key(3)=3
|
|---|
| 2141 | c
|
|---|
| 2142 | c Count number of unique vertices
|
|---|
| 2143 | nlv = 2**ldim
|
|---|
| 2144 | ngvv = iglmax(vertex,nlv*nel)
|
|---|
| 2145 | ngv = ngvv
|
|---|
| 2146 | c
|
|---|
| 2147 | c Assign hypercube ordering of vertices.
|
|---|
| 2148 | do e=1,nel
|
|---|
| 2149 | do j=0,1
|
|---|
| 2150 | do i=0,1
|
|---|
| 2151 | c Local to global node number (vertex)
|
|---|
| 2152 | il = 1 + (nx-1)*i + nx*(nx-1)*j
|
|---|
| 2153 | ile = il + nx*ny*(e-1)
|
|---|
| 2154 | glo_num(ile) = vertex(i,j,e)
|
|---|
| 2155 | enddo
|
|---|
| 2156 | enddo
|
|---|
| 2157 | enddo
|
|---|
| 2158 | if (nx.eq.2) return
|
|---|
| 2159 | c
|
|---|
| 2160 | c Assign edge labels by bounding vertices.
|
|---|
| 2161 | do e=1,nel
|
|---|
| 2162 | do j=0,1
|
|---|
| 2163 | do i=0,1
|
|---|
| 2164 | edge(i,j,1,e) = vertex(i,j,e) ! r-edge
|
|---|
| 2165 | edge(j,i,2,e) = vertex(i,j,e) ! s-edge
|
|---|
| 2166 | enddo
|
|---|
| 2167 | enddo
|
|---|
| 2168 | enddo
|
|---|
| 2169 |
|
|---|
| 2170 | c Sort edges by bounding vertices.
|
|---|
| 2171 | do i=0,4*nel-1
|
|---|
| 2172 | if (edge(0,i,1,1).gt.edge(1,i,1,1)) then
|
|---|
| 2173 | kswap = edge(0,i,1,1)
|
|---|
| 2174 | edge(0,i,1,1) = edge(1,i,1,1)
|
|---|
| 2175 | edge(1,i,1,1) = kswap
|
|---|
| 2176 | endif
|
|---|
| 2177 | etuple(3,i+1) = edge(0,i,1,1)
|
|---|
| 2178 | etuple(4,i+1) = edge(1,i,1,1)
|
|---|
| 2179 | enddo
|
|---|
| 2180 |
|
|---|
| 2181 | c Assign a number (rank) to each unique edge
|
|---|
| 2182 | m = 4
|
|---|
| 2183 | n = 4*nel
|
|---|
| 2184 | nmax = 4*lelt*nsafe ! nsafe for crystal router factor of safety
|
|---|
| 2185 |
|
|---|
| 2186 | call gbtuple_rank(etuple,m,n,nmax,cr_h,nid,np,ind)
|
|---|
| 2187 | do i=1,4*nel
|
|---|
| 2188 | enum(i,1) = etuple(3,i)
|
|---|
| 2189 | enddo
|
|---|
| 2190 | n_unique_edges = iglmax(enum,4*nel)
|
|---|
| 2191 |
|
|---|
| 2192 | c= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
|
|---|
| 2193 | c Assign global vertex numbers to SEM nodes on each edge
|
|---|
| 2194 | n_on_edge = nx-2
|
|---|
| 2195 | do e=1,nel
|
|---|
| 2196 |
|
|---|
| 2197 | iedg_loc = 0
|
|---|
| 2198 |
|
|---|
| 2199 | c Edges 1-2
|
|---|
| 2200 | do j=0,1
|
|---|
| 2201 | igv = ngv + n_on_edge*(enum(iedg_loc+1,e)-1)
|
|---|
| 2202 | i0 = nx*(nx-1)*j
|
|---|
| 2203 | i0e = i0 + nxyz*(e-1)
|
|---|
| 2204 | if (glo_num(i0e+1).lt.glo_num(i0e+nx)) then
|
|---|
| 2205 | do i=2,nx-1 ! std forward case
|
|---|
| 2206 | glo_num(i0e+i) = igv + i-1
|
|---|
| 2207 | enddo
|
|---|
| 2208 | else
|
|---|
| 2209 | do i=2,nx-1 ! backward case
|
|---|
| 2210 | glo_num(i0e+i) = igv + 1 + n_on_edge-(i-1)
|
|---|
| 2211 | enddo
|
|---|
| 2212 | endif
|
|---|
| 2213 | iedg_loc = iedg_loc + 1
|
|---|
| 2214 | enddo
|
|---|
| 2215 | c
|
|---|
| 2216 | c Edges 3-4
|
|---|
| 2217 | do i=0,1
|
|---|
| 2218 | igv = ngv + n_on_edge*(enum(iedg_loc+1,e)-1)
|
|---|
| 2219 | i0 = 1+(nx-1)*i
|
|---|
| 2220 | i0e = i0 + nxyz*(e-1)
|
|---|
| 2221 | if (glo_num(i0e).lt.glo_num(i0e+nx*(nx-1))) then
|
|---|
| 2222 | do j=2,nx-1 ! std forward case
|
|---|
| 2223 | glo_num(i0e+(j-1)*nx) = igv + j-1
|
|---|
| 2224 | enddo
|
|---|
| 2225 | else
|
|---|
| 2226 | do j=2,nx-1 ! backward case
|
|---|
| 2227 | glo_num(i0e+(j-1)*nx) = igv + 1 + n_on_edge-(j-1)
|
|---|
| 2228 | enddo
|
|---|
| 2229 | endif
|
|---|
| 2230 | iedg_loc = iedg_loc + 1
|
|---|
| 2231 | enddo
|
|---|
| 2232 | enddo
|
|---|
| 2233 |
|
|---|
| 2234 | ngve = n_unique_edges*n_on_edge
|
|---|
| 2235 | ngv = ngv + ngve
|
|---|
| 2236 | c
|
|---|
| 2237 | c Finally, number interiors
|
|---|
| 2238 | c
|
|---|
| 2239 | n_in_interior = (nx-2)*(ny-2)
|
|---|
| 2240 | ngvi = n_in_interior*nelgt
|
|---|
| 2241 | if (ifcenter) then
|
|---|
| 2242 | do e=1,nel
|
|---|
| 2243 | ig0 = ngv + n_in_interior*(lglel(e)-1)
|
|---|
| 2244 | l = 0
|
|---|
| 2245 | do j=2,ny-1
|
|---|
| 2246 | do i=2,nx-1
|
|---|
| 2247 | l = l+1
|
|---|
| 2248 | glo_num(i+nx*(j-1)+nxyz*(e-1)) = ig0+l
|
|---|
| 2249 | enddo
|
|---|
| 2250 | enddo
|
|---|
| 2251 | enddo
|
|---|
| 2252 | ngv = ngv + ngvi
|
|---|
| 2253 | else
|
|---|
| 2254 | do e=1,nel
|
|---|
| 2255 | l = 0
|
|---|
| 2256 | do j=2,ny-1
|
|---|
| 2257 | do i=2,nx-1
|
|---|
| 2258 | l = l+1
|
|---|
| 2259 | glo_num(i+nx*(j-1)+nxyz*(e-1)) = 0
|
|---|
| 2260 | enddo
|
|---|
| 2261 | enddo
|
|---|
| 2262 | enddo
|
|---|
| 2263 | endif
|
|---|
| 2264 |
|
|---|
| 2265 | c
|
|---|
| 2266 | c Quick check on maximum #dofs:
|
|---|
| 2267 | m = nxyz*nelt
|
|---|
| 2268 | ngvm = i8glmax(glo_num,m)
|
|---|
| 2269 | ngvv = ngvv + ngve ! number of unique ids w/o interior
|
|---|
| 2270 | ngvi = ngvi + ngvv ! total number of unique ids
|
|---|
| 2271 | if (nio.eq.0) write(6,1) nx,ngvv,ngvi,ngv,ngvm
|
|---|
| 2272 | 1 format(' setvert2d:',i4,4i12)
|
|---|
| 2273 | c
|
|---|
| 2274 | return
|
|---|
| 2275 | end
|
|---|
| 2276 | c-----------------------------------------------------------------------
|
|---|
| 2277 | c
|
|---|
| 2278 | subroutine map_to_crs(a,na,b,nb,if3d,w,ldw)
|
|---|
| 2279 | c
|
|---|
| 2280 | c Input: b
|
|---|
| 2281 | c Output: a
|
|---|
| 2282 | c
|
|---|
| 2283 | real a(1),b(1),w(1)
|
|---|
| 2284 | logical if3d
|
|---|
| 2285 | c
|
|---|
| 2286 | parameter(lx=40)
|
|---|
| 2287 | real za(lx),zb(lx)
|
|---|
| 2288 | c
|
|---|
| 2289 | real iba(lx*lx),ibat(lx*lx)
|
|---|
| 2290 | save iba,ibat
|
|---|
| 2291 | c
|
|---|
| 2292 | integer nao,nbo
|
|---|
| 2293 | save nao,nbo
|
|---|
| 2294 | data nao,nbo / -9, -9/
|
|---|
| 2295 | c
|
|---|
| 2296 | if (na.gt.lx.or.nb.gt.lx) then
|
|---|
| 2297 | write(6,*)'ERROR: increase lx in map_to_crs to max:',na,nb
|
|---|
| 2298 | call exitt
|
|---|
| 2299 | endif
|
|---|
| 2300 | c
|
|---|
| 2301 | if (na.ne.nao .or. nb.ne.nbo) then
|
|---|
| 2302 | nao = na
|
|---|
| 2303 | nbo = nb
|
|---|
| 2304 | call zwgll(za,w,na)
|
|---|
| 2305 | call zwgll(zb,w,nb)
|
|---|
| 2306 | call igllm(iba,ibat,zb,za,nb,na,nb,na)
|
|---|
| 2307 | endif
|
|---|
| 2308 | c
|
|---|
| 2309 | call specmpn(a,na,b,nb,iba,ibat,if3d,w,ldw)
|
|---|
| 2310 |
|
|---|
| 2311 | return
|
|---|
| 2312 | end
|
|---|
| 2313 | c-----------------------------------------------------------------------
|
|---|
| 2314 | subroutine check_p_bc(glo_num,nx,ny,nz,nel)
|
|---|
| 2315 |
|
|---|
| 2316 | include 'SIZE'
|
|---|
| 2317 | include 'TOTAL'
|
|---|
| 2318 |
|
|---|
| 2319 | integer*8 glo_num(nx,ny,nz,nel)
|
|---|
| 2320 | integer*8 gmn
|
|---|
| 2321 |
|
|---|
| 2322 | integer e,f,fo,ef,efo
|
|---|
| 2323 | integer eface0(6)
|
|---|
| 2324 | save eface0
|
|---|
| 2325 | data eface0 / 4,2,1,3,5,6 /
|
|---|
| 2326 |
|
|---|
| 2327 | ifld = 2
|
|---|
| 2328 | if (ifflow) ifld = 1
|
|---|
| 2329 |
|
|---|
| 2330 | nface=2*ldim
|
|---|
| 2331 | do e=1,nelt
|
|---|
| 2332 | do f=1,nface,2
|
|---|
| 2333 | fo = f+1
|
|---|
| 2334 | ef = eface0(f)
|
|---|
| 2335 | efo = eface0(fo)
|
|---|
| 2336 | if (cbc(ef,e,ifld).eq.'p '.and.cbc(efo,e,ifld).eq.'p ') then
|
|---|
| 2337 | if (f.eq.1) then ! r=-/+1
|
|---|
| 2338 | do k=1,nz
|
|---|
| 2339 | do j=1,ny
|
|---|
| 2340 | gmn = min(glo_num(1,j,k,e),glo_num(nx,j,k,e))
|
|---|
| 2341 | glo_num(1 ,j,k,e) = gmn
|
|---|
| 2342 | glo_num(nx,j,k,e) = gmn
|
|---|
| 2343 | enddo
|
|---|
| 2344 | enddo
|
|---|
| 2345 | elseif (f.eq.3) then ! s=-/+1
|
|---|
| 2346 | do k=1,nz
|
|---|
| 2347 | do i=1,nx
|
|---|
| 2348 | gmn = min(glo_num(i,1,k,e),glo_num(i,ny,k,e))
|
|---|
| 2349 | glo_num(i,1 ,k,e) = gmn
|
|---|
| 2350 | glo_num(i,ny,k,e) = gmn
|
|---|
| 2351 | enddo
|
|---|
| 2352 | enddo
|
|---|
| 2353 | else
|
|---|
| 2354 | do j=1,ny
|
|---|
| 2355 | do i=1,nx
|
|---|
| 2356 | gmn = min(glo_num(i,j,1,e),glo_num(i,j,nz,e))
|
|---|
| 2357 | glo_num(i,j,1 ,e) = gmn
|
|---|
| 2358 | glo_num(i,j,nz,e) = gmn
|
|---|
| 2359 | enddo
|
|---|
| 2360 | enddo
|
|---|
| 2361 | endif
|
|---|
| 2362 | endif
|
|---|
| 2363 | enddo
|
|---|
| 2364 | enddo
|
|---|
| 2365 |
|
|---|
| 2366 | return
|
|---|
| 2367 | end
|
|---|
| 2368 | c-----------------------------------------------------------------------
|
|---|
