KTH framework for Nek5000 toolboxes; testing version  0.0.1
navier8.f
Go to the documentation of this file.
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-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
igl_running_sum
function igl_running_sum(in)
Definition: comm_mpi.f:700
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
facind
subroutine facind(kx1, kx2, ky1, ky2, kz1, kz2, nx, ny, nz, iface)
Definition: connect1.f:1028
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
get_vert
subroutine get_vert
Definition: map2.f:126
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
icopy
subroutine icopy(a, b, n)
Definition: math.f:289
isort
subroutine isort(a, ind, n)
Definition: math.f:1273
copy
subroutine copy(a, b, n)
Definition: math.f:260
blank
subroutine blank(A, N)
Definition: math.f:19
iglmax
function iglmax(a, n)
Definition: math.f:913
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
ltrunc
function ltrunc(string, l)
Definition: math.f:494
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
iglmax_ms
function iglmax_ms(a, n)
Definition: math.f:1729
rone
subroutine rone(a, n)
Definition: math.f:230
invcol1
subroutine invcol1(a, n)
Definition: math.f:62
rzero
subroutine rzero(a, n)
Definition: math.f:208
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
iunswap
subroutine iunswap(b, ind, n, temp)
Definition: navier6.f:135
check_p_bc
subroutine check_p_bc(glo_num, nx, ny, nz, nel)
Definition: navier8.f:2315
irank
subroutine irank(A, IND, N)
Definition: navier8.f:984
map_c_to_f_l2_bilin
subroutine map_c_to_f_l2_bilin(uf, uc, w)
Definition: navier8.f:1080
crs_solve_h1
subroutine crs_solve_h1(uf, vf)
Definition: navier8.f:1341
maph1_to_l2
subroutine maph1_to_l2(a, na, b, nb, if3d, w, ldw)
Definition: navier8.f:1127
map_c_to_f_h1_bilin
subroutine map_c_to_f_h1_bilin(uf, uc)
Definition: navier8.f:1407
get_vertex
subroutine get_vertex
Definition: navier8.f:1629
a_crs_3d
subroutine a_crs_3d(a, h1, h2, xc, yc, zc, ie)
Definition: navier8.f:651
map_f_to_c_l2_bilin
subroutine map_f_to_c_l2_bilin(uc, uf, w)
Definition: navier8.f:1103
tuple_sort
subroutine tuple_sort(a, lda, n, key, nkey, ind, aa)
Definition: navier8.f:389
bindec
subroutine bindec(bin_in)
Definition: navier8.f:712
set_vert
subroutine set_vert(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:4
get_local_crs
subroutine get_local_crs(a, lda, nxc, h1, h2, w, ldw)
Definition: navier8.f:510
a_crs_2d
subroutine a_crs_2d(a, h1, h2, x, y, ie)
Definition: navier8.f:828
iftuple_altb
logical function iftuple_altb(a, b, key, nkey)
Definition: navier8.f:473
specmpn
subroutine specmpn(b, nb, a, na, ba, ab, if3d, w, ldw)
Definition: navier8.f:941
map_to_crs
subroutine map_to_crs(a, na, b, nb, if3d, w, ldw)
Definition: navier8.f:2279
gen_crs_basis2
subroutine gen_crs_basis2(b, j)
Definition: navier8.f:1584
get_local_a_tet
subroutine get_local_a_tet(a, x, y, z, kt, ie)
Definition: navier8.f:731
a_crs_enriched
subroutine a_crs_enriched(a, h1, h2, x1, y1, z1, nxc, if3d, ie)
Definition: navier8.f:560
iftuple_ialtb
logical function iftuple_ialtb(a, b, key, nkey)
Definition: navier8.f:453
out_se0
subroutine out_se0(se2crs, nx, nel, name)
Definition: navier8.f:1311
crs_solve_l2
subroutine crs_solve_l2(uf, vf)
Definition: navier8.f:38
set_up_h1_crs
subroutine set_up_h1_crs
Definition: navier8.f:84
maph1_to_l2t
subroutine maph1_to_l2t(b, nb, a, na, if3d, w, ldw)
Definition: navier8.f:1166
map_m_to_n
subroutine map_m_to_n(a, na, b, nb, if3d, w, ldw)
Definition: navier8.f:903
out_se1
subroutine out_se1(se2crs, nx, name)
Definition: navier8.f:1281
setvert3d
subroutine setvert3d(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:1744
ifacev_redef
subroutine ifacev_redef(a, ie, iface, val, nx, ny, nz)
Definition: navier8.f:1063
map_f_to_c_h1_bilin
subroutine map_f_to_c_h1_bilin(uc, uf)
Definition: navier8.f:1456
irank_vec_tally
subroutine irank_vec_tally(ind, nn, a, m, n, key, nkey, key2, aa)
Definition: navier8.f:1204
setvert2d
subroutine setvert2d(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:2103
gen_crs_basis
subroutine gen_crs_basis(b, j)
Definition: navier8.f:1540
get_local_crs_galerkin
subroutine get_local_crs_galerkin(a, ncl, nxc, h1, h2, w1, w2)
Definition: navier8.f:1499
irank_vecn
subroutine irank_vecn(ind, nn, a, m, n, key, nkey, aa)
Definition: navier8.f:1642
iranku
subroutine iranku(r, input, n, w, ind)
Definition: navier8.f:1033
iftuple_ianeb
logical function iftuple_ianeb(a, b, key, nkey)
Definition: navier8.f:493
set_h1_basis_bilin
subroutine set_h1_basis_bilin
Definition: navier8.f:1390
set_mat_ij
subroutine set_mat_ij(ia, ja, n, ne)
Definition: navier8.f:247
gbtuple_rank
subroutine gbtuple_rank(tuple, m, n, nmax, cr_h, nid, np, ind)
Definition: navier8.f:1684
ituple_sort
subroutine ituple_sort(a, lda, n, key, nkey, ind, aa)
Definition: navier8.f:327
irank_vec
subroutine irank_vec(ind, nn, a, m, n, key, nkey, aa)
Definition: navier8.f:263
set_jl_crs_mask
subroutine set_jl_crs_mask(n, mask, se_to_gcrs)
Definition: navier8.f:238
zwgl
subroutine zwgl(Z, W, NP)
Definition: speclib.f:91
zwgll
subroutine zwgll(Z, W, NP)
Definition: speclib.f:108
igllm
subroutine igllm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2)
Definition: speclib.f:1102