KTH framework for Nek5000 toolboxes; testing version  0.0.1
mvmesh.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine cbcmesh
3 C
4 C Generate boundary conditions (CBC arrays) for mesh solver
5 C
6  include 'SIZE'
7  include 'GEOM'
8  include 'INPUT'
9  include 'TSTEP'
10 C
11  CHARACTER CBM*1,CBF*3,CBT*3,CB*3
12 C
13  ifld = 0
14  nface = 2*ldim
15  ifmelt = .false.
16  IF (iftmsh(ifld)) ifmelt=.true.
17 C
18  DO 100 iel=1,nelt
19  DO 100 ifc=1,nface
20  cbm = cbc(ifc,iel,0)
21  cbf = cbc(ifc,iel,1)
22  cbt = cbc(ifc,iel,2)
23 C
24  IF (cbt(1:1).EQ.'M') THEN
25  ifld = 2
26  cb = cbt
27  GOTO 200
28  ENDIF
29  IF (cbf(1:1).EQ.'M' .OR. cbf(1:1).EQ.'m') THEN
30  ifld = 1
31  cb = cbf
32 C IF (CBF.EQ.'mv ' .AND. CBT.EQ.'E ') THEN
33 C IFTMSH(0)=.TRUE.
34 C ENDIF
35  IF (cbf.EQ.'mv ' .AND. cbm.EQ.'+' ) THEN
36  cb = 'mvn'
37  cbc(ifc,iel,1) = cb
38  ENDIF
39  GOTO 200
40  ENDIF
41  IF (cbf(1:1).EQ.'V' .OR. cbf(1:1).EQ.'v' .OR.
42  $ cbf(1:1).EQ.'W' ) THEN
43  ifld = 1
44  cb = 'FIX'
45  IF (ifmelt .OR. cbm.EQ.'+') cb='SYM'
46  GOTO 200
47  ENDIF
48  IF (cbt.EQ.'T ' .OR. cbt.EQ.'t ') THEN
49  ifld = 2
50  cb = 'FIX'
51  IF (cbm.EQ.'+') cb='SYM'
52  GOTO 200
53  ENDIF
54  IF (cbf.EQ.'P ' .OR. cbf.EQ.'E ') THEN
55  ifld = 1
56  cb = cbf
57  IF (cbm.EQ.'-') cb='FIX'
58  IF (cbm.EQ.'+') cb='SYM'
59  GOTO 200
60  ENDIF
61  IF (cbt.EQ.'P ' .OR. cbt.EQ.'E ') THEN
62  ifld = 2
63  cb = cbt
64  IF (cbm.EQ.'-') cb='FIX'
65  IF (cbm.EQ.'+') cb='SYM'
66  GOTO 200
67  ENDIF
68  ifld = 1
69  IF (cbf.EQ.' ') ifld = 2
70  cb = 'SYM'
71  IF (cbm.EQ.'-') cb = 'FIX'
72 C
73  200 cbc(ifc,iel,0) = cb
74  DO 250 i=1,5
75  250 bc(i,ifc,iel,0)=bc(i,ifc,iel,ifld)
76  100 CONTINUE
77 C
78  return
79  end
80 c-----------------------------------------------------------------------
81  subroutine admeshv
82 C
83  include 'SIZE'
84  include 'SOLN'
85  include 'TSTEP'
86 C
87  COMMON /scruz/ fm1(lx1,ly1,lz1,lelt)
88  $ , fm2(lx1,ly1,lz1,lelt)
89  $ , fm3(lx1,ly1,lz1,lelt)
90  $ , phi(lx1,ly1,lz1,lelt)
91 C
92  ntot1=lx1*ly1*lz1*nelv
93 C
94  CALL rzero (fm1,ntot1)
95  CALL rzero (fm2,ntot1)
96  CALL rzero (fm3,ntot1)
97 C
98  CALL divws (fm1,vx,phi,nelv,1)
99  CALL divws (fm2,vy,phi,nelv,2)
100  CALL add2 (bfx,fm1,ntot1)
101  CALL add2 (bfy,fm2,ntot1)
102  IF (ldim.EQ.3) THEN
103  CALL divws (fm3,vz,phi,nelv,3)
104  CALL add2 (bfz,fm3,ntot1)
105  ENDIF
106 C
107  return
108  end
109 c-----------------------------------------------------------------------
110  subroutine admesht
111 C
112  include 'SIZE'
113  include 'SOLN'
114  include 'TSTEP'
115 C
116  COMMON /scruz/ fmt(lx1,ly1,lz1,lelt)
117  $ , phi(lx1,ly1,lz1,lelt)
118 C
119  ifld = 0
120  nel = nelfld(ifld)
121  ntot1= lx1*ly1*lz1*nel
122 C
123  CALL rzero (fmt,ntot1)
124  CALL divws (fmt,t(1,1,1,1,ifield-1),phi,nel,1)
125  CALL addcol3 (bq(1,1,1,1,ifield-1),fmt,vtrans(1,1,1,1,ifield),
126  & ntot1)
127 C
128  return
129  end
130 c-----------------------------------------------------------------------
131  subroutine divws (fms,sfv,phi,nel,idir)
132 C
133  include 'SIZE'
134  include 'GEOM'
135  include 'MASS'
136  include 'MVGEOM'
137  include 'WZ'
138  include 'INPUT'
139 C
140  COMMON /scrsf/ phr(lx1,ly1,lz1,lelt)
141  $ , phs(lx1,ly1,lz1,lelt)
142  $ , pht(lx1,ly1,lz1,lelt)
143 
144  dimension fms(lx1,ly1,lz1,1)
145  $ , sfv(lx1,ly1,lz1,1)
146  $ , phi(lx1,ly1,lz1,1)
147 C
148  nxyz1 = lx1*ly1*lz1
149  ntot1 = nxyz1*nel
150 C
151  CALL col3 (phi,sfv,wx,ntot1)
152  CALL urst (phi,phr,phs,pht,nel)
153  CALL addcol3 (fms,rxm1,phr,ntot1)
154  CALL addcol3 (fms,sxm1,phs,ntot1)
155  IF (ldim.EQ.3) CALL addcol3 (fms,txm1,pht,ntot1)
156 C
157  CALL col3 (phi,sfv,wy,ntot1)
158  CALL urst (phi,phr,phs,pht,nel)
159  CALL addcol3 (fms,rym1,phr,ntot1)
160  CALL addcol3 (fms,sym1,phs,ntot1)
161  IF (ldim.EQ.3) CALL addcol3 (fms,tym1,pht,ntot1)
162 C
163  IF (ldim.EQ.3) THEN
164  CALL col3 (phi,sfv,wz,ntot1)
165  CALL urst (phi,phr,phs,pht,nel)
166  CALL addcol3 (fms,rzm1,phr,ntot1)
167  CALL addcol3 (fms,szm1,phs,ntot1)
168  CALL addcol3 (fms,tzm1,pht,ntot1)
169  ENDIF
170 C
171  CALL col2 (fms,bm1,ntot1)
172  CALL invcol2 (fms,jacm1,ntot1)
173 C
174  IF (ifaxis) CALL axifms (fms,sfv,phi,nel,idir)
175 C
176  return
177  end
178 c-----------------------------------------------------------------------
179  subroutine axifms (fms,sfv,phi,nel,idir)
180 C
181  include 'SIZE'
182  include 'DXYZ'
183  include 'GEOM'
184  include 'INPUT'
185  include 'MASS'
186  include 'MVGEOM'
187  include 'WZ'
188  COMMON /scrsf/ phr(lx1,ly1,lz1,lelt)
189  $ , phs(lx1,ly1,lz1,lelt)
190  $ , pht(lx1,ly1,lz1,lelt)
191 C
192  dimension fms(lx1,ly1,lz1,1)
193  $ , phi(lx1,ly1,lz1,1)
194  $ , sfv(lx1,ly1,lz1,1)
195  $ , wys(lx1)
196  equivalence(wys(1),pht(1,1,1,1))
197 C
198  nxyz1 = lx1*ly1*lz1
199  ntot1 = nxyz1*nel
200  CALL col3 (phi,sfv,wy,ntot1)
201 C
202  DO 100 iel=1,nel
203  IF ( ifrzer(iel) ) THEN
204  IF (idir.EQ.1) THEN
205  CALL mxm (wy(1,1,1,iel),lx1,datm1,ly1,wys,1)
206  DO 220 ix=1,lx1
207  fms(ix,1,1,iel)= fms(ix,1,1,iel) + wxm1(ix)*wam1(1)*
208  $ wys(ix)*sfv(ix,1,1,iel)*jacm1(ix,1,1,iel)
209  220 CONTINUE
210  ENDIF
211  DO 320 ix=1,lx1
212  DO 320 iy=2,ly1
213  fms(ix,iy,1,iel)=fms(ix,iy,1,iel) + phi(ix,iy,1,iel) *
214  $ bm1(ix,iy,1,iel) / ym1(ix,iy,1,iel)
215  320 CONTINUE
216  ELSE
217  CALL addcol4 (fms(1,1,1,iel),phi(1,1,1,iel),jacm1(1,1,1,iel),
218  $ w2cm1,nxyz1)
219  ENDIF
220  100 CONTINUE
221 C
222  return
223  end
224 c-----------------------------------------------------------------------
225  subroutine updcoor
226 C-----------------------------------------------------------------------
227 C
228 C Subroutine to update geometry for moving boundary problems
229 C
230 C-----------------------------------------------------------------------
231  include 'SIZE'
232  include 'TSTEP'
233  include 'INPUT'
234 C
235  ifield = 0
236  nel = nelfld(ifield)
237 C
238 C Update collocation points coordinates
239 C
240  CALL updxyz (nel)
241 C
242 C Shift lagged mesh velocity
243 C
244  if (.not.ifrich) call lagmshv (nel)
245 C
246  return
247  end
248 c-----------------------------------------------------------------------
249  subroutine mvbdry (nel)
250 
251 c Evaluate mesh velocities at all moving boundaries
252 
253  include 'SIZE'
254  include 'GEOM'
255  include 'INPUT'
256  include 'MVGEOM'
257  include 'SOLN'
258  include 'TSTEP'
259 
260  common /scrsf/ wvx(lx1*ly1*lz1,lelt)
261  $ , wvy(lx1*ly1*lz1,lelt)
262  $ , wvz(lx1*ly1*lz1,lelt)
263  common /scrch/ wtx(lx1*ly1*lz1,lelt)
264  $ , wty(lx1*ly1*lz1,lelt)
265  common /scrmg/ wtz(lx1*ly1*lz1,lelt)
266  $ , rnx(lx1*ly1*lz1,lelt)
267  $ , rny(lx1*ly1*lz1,lelt)
268  $ , rnz(lx1*ly1*lz1,lelt)
269  common /scruz/ dsa(lx1*ly1*lz1,lelt)
270  $ , qni(lx1*ly1*lz1,lelt)
271  $ , smt(lx1*ly1*lz1,lelt)
272  $ , ta(lx1*ly1*lz1,lelt)
273 
274  logical ifalgn,ifnorx,ifnory,ifnorz,ifdsmv,ifregw
275  character cb*3
276  integer e,f
277 
278  ifield = 0
279  nxyz1 = lx1*ly1*lz1
280  n = lx1*ly1*lz1*nel
281  nface = 2*ldim
282  call rzero3 (rnx,rny,rnz,n)
283 
284  do 100 e=1,nel
285  do 100 f=1,nface
286  cb = cbc(f,e,ifield)
287  if (cb.eq.'ms ' .or. cb.eq.'MS ' .or.
288  $ cb.eq.'msi' .or. cb.eq.'MSI' .or.
289  $ cb.eq.'mm ' .or. cb.eq.'MM ' .or.
290  $ cb.eq.'mv ' .OR. cb.eq.'mvn' .or.
291  $ cb.eq.'MLI') then
292  call facexv (unx(1,1,f,e),uny(1,1,f,e),
293  $ unz(1,1,f,e),rnx(1,e),
294  $ rny(1,e),rnz(1,e),f,1)
295  endif
296  100 continue
297 
298  call opdssum (rnx,rny,rnz)
299  call unitvec (rnx,rny,rnz,n)
300 
301  call rzero3 (wvx,wvy,wvz,n)
302  call rzero3 (wtx,wty,wtz,n)
303  do 1000 isweep=1,2
304 
305  ifregw = .false.
306  ifdsmv = .false.
307  call rzero (dsa,n)
308  call rzero (ta,n)
309 
310  if (ifflow) then
311  ifield = 1
312  call rzero (smt,n)
313  do 210 e=1,nelv
314  do 210 f=1,nface
315  cb = cbc(f,e,ifield)
316  if (cb.eq.'mv ' .or. cb.eq.'mvn' .or.
317  $ cb.eq.'mm ' .or. cb.eq.'MM ' .or.
318  $ cb.eq.'msi' .or. cb.eq.'MSI' .or.
319  $ cb.eq.'ms ' .or. cb.eq.'MS ') then
320  ifregw = .true.
321  call facec3 (wvx(1,e),wvy(1,e),wvz(1,e),
322  $ vx(1,1,1,e),vy(1,1,1,e),vz(1,1,1,e),f)
323  if (cb.ne.'mv ')
324  $ call norcmp2(wvx(1,e),wvy(1,e),wvz(1,e),e,f)
325  endif
326 c if (cb.eq.'msi' .or. cb.eq.'MSI') then
327 c ifdsmv = .true.
328 c call facsmt (smt(1,e),f)
329 c endif
330  210 continue
331 
332  call dsavg(wvx)
333  call dsavg(wvy)
334  if (if3d) call dsavg(wvz)
335 
336  if (istep.eq.0) call opcopy(wx,wy,wz,wvx,wvy,wvz)
337 c if (istep.eq.4) then
338 c call opcopy(wx,wy,wz,wvx,wvy,wvz)
339 c call outpost(vx,vy,vz,wtx,wtx,' ')
340 c call outpost(wx,wy,wz,wtx,wtx,' ')
341 c write(6,*) 'quit1',istep
342 c stop
343 c endif
344 
345  iregw = 0 ! Global handshake on logicals ifregw, ifdsmv
346  if (ifregw) iregw = 1 ! pff 9/11/07
347  iregw = iglmax(iregw,1)
348  if (iregw.eq.1) ifregw = .true.
349 
350  idsmv = 0
351  if (ifdsmv) idsmv = 1
352  idsmv = iglmax(idsmv,1)
353  if (idsmv.eq.1) ifdsmv = .true.
354 
355  ifdsmv = .false.
356  ifregw = .true.
357 
358 c if (ifdsmv) then
359 c call dssum (smt,lx1,ly1,lz1)
360 c do 215 e=1,nelv
361 c do 215 f=1,nface
362 c cb = cbc(f,e,ifield)
363 c if (cb.eq.'msi' .or. cb.eq.'MSI') then
364 c call facec3 (wtx(1,e),wty(1,e),wtz(1,e),
365 c $ vx(1,1,1,e),vy(1,1,1,e),vz(1,1,1,e),f)
366 c call facemv (wtx(1,e),wty(1,e),wtz(1,e),
367 c $ rnx(1,e),rny(1,e),rnz(1,e),
368 c $ smt(1,e),f)
369 c endif
370 c 215 continue
371 c call opdssum(wtx,wty,wtz)
372 c endif
373 
374  endif
375 C
376  if (ifmelt .and. istep.gt.0) then
377  ifield = 2
378  call rzero (smt,n)
379  call cqnet (qni,ta,nel)
380  do 220 e=1,nelt
381  do 220 f=1,nface
382  cb = cbc(f,e,ifield)
383  if (cb.eq.'MLI') call facsmt (smt(1,e),f)
384  if (cb.eq.'MLI' .or. cb.eq.'MCI') then
385  call facexs (area(1,1,f,e),ta,f,1)
386  call add2 (dsa(1,e),ta,nxyz1)
387  endif
388  220 continue
389  call dssum (smt,lx1,ly1,lz1)
390  call dssum (dsa,lx1,ly1,lz1)
391  do 280 e=1,nelt
392  do 280 f=1,nface
393  cb = cbc(f,e,ifield)
394  if (cb.eq.'MLI') then
395  rhola = -0.5 * bc(5,f,e,ifield)
396  call facemt (wtx(1,e),wty(1,e),wtz(1,e),
397  $ rnx(1,e),rny(1,e),rnz(1,e),
398  $ qni(1,e),dsa(1,e),smt(1,e),
399  $ rhola,f)
400  endif
401  280 continue
402  call opdssum (wtx,wty,wtz)
403  endif
404 
405  ifield = 0
406  do 330 e=1,nel
407  do 330 f=1,nface
408  cb = cbc(f,e,ifield)
409  if (cb.eq.'SYM') then
410  call chknord (ifalgn,ifnorx,ifnory,ifnorz,f,e)
411  if (ifregw) then
412  if (ifnorx) call facev (wvx,e,f,0.0,lx1,ly1,lz1)
413  if (ifnory) call facev (wvy,e,f,0.0,lx1,ly1,lz1)
414  if (ifnorz) call facev (wvz,e,f,0.0,lx1,ly1,lz1)
415  if (.not.ifalgn) call faczqn (wvx(1,e),wvy(1,e),
416  $ wvz(1,e),f,e)
417  endif
418  if (ifdsmv .or. ifmelt) then
419  if (ifnorx) call facev (wtx,e,f,0.0,lx1,ly1,lz1)
420  if (ifnory) call facev (wty,e,f,0.0,lx1,ly1,lz1)
421  if (ifnorz) call facev (wtz,e,f,0.0,lx1,ly1,lz1)
422  if (.not.ifalgn) call faczqn (wtx(1,e),wty(1,e),
423  $ wtz(1,e),f,e)
424  endif
425  endif
426  330 continue
427 
428  do 350 e=1,nel
429  do 350 f=1,nface
430  cb = cbc(f,e,ifield)
431  if (cb.eq.'FIX') then
432  if (ifregw) then
433  call facev (wvx,e,f,0.0,lx1,ly1,lz1)
434  call facev (wvy,e,f,0.0,lx1,ly1,lz1)
435  if (ldim.eq.3) call facev (wvz,e,f,0.0,lx1,ly1,lz1)
436  endif
437  if (ifdsmv .or. ifmelt) then
438  call facev (wtx,e,f,0.0,lx1,ly1,lz1)
439  call facev (wty,e,f,0.0,lx1,ly1,lz1)
440  if (ldim.eq.3) call facev (wtz,e,f,0.0,lx1,ly1,lz1)
441  endif
442  endif
443  350 continue
444 
445  if (isweep.eq.1) then
446  if (ifregw) then
447  call dsop (wvx,'MXA',lx1,ly1,lz1)
448  call dsop (wvy,'MXA',lx1,ly1,lz1)
449  if (ldim.eq.3) call dsop (wvz,'MXA',lx1,ly1,lz1)
450  endif
451  if (ifdsmv .or. ifmelt) then
452  call dsop (wtx,'MXA',lx1,ly1,lz1)
453  call dsop (wty,'MXA',lx1,ly1,lz1)
454  if (ldim.eq.3) call dsop (wtz,'MXA',lx1,ly1,lz1)
455  endif
456  else
457  if (ifregw) then
458  call dsop (wvx,'MNA',lx1,ly1,lz1)
459  call dsop (wvy,'MNA',lx1,ly1,lz1)
460  if (ldim.eq.3) call dsop (wvz,'MNA',lx1,ly1,lz1)
461  endif
462  if (ifdsmv .or. ifmelt) then
463  call dsop (wtx,'MNA',lx1,ly1,lz1)
464  call dsop (wty,'MNA',lx1,ly1,lz1)
465  if (ldim.eq.3) call dsop (wtz,'MNA',lx1,ly1,lz1)
466  endif
467  endif
468 
469  1000 continue
470 
471  call rmask (wx,wy,wz,nel)
472 c if (istep.eq.2) then
473 c call outpost(wx,wy,wz,wtx,wtx,' ')
474 c write(6,*) 'quit2'
475 c stop
476 c endif
477 
478  if (ifregw) then
479  call add2 (wx,wvx,n)
480  call add2 (wy,wvy,n)
481  if (ldim.eq.3) call add2 (wz,wvz,n)
482  endif
483  if (ifdsmv .or. ifmelt) then
484  call add2 (wx,wtx,n)
485  call add2 (wy,wty,n)
486  if (ldim.eq.3) call add2 (wz,wtz,n)
487  endif
488 
489 c if (istep.gt.4) then
490 c call outpost(wx,wy,wz,wtx,wtx,' ')
491 c write(6,*) 'quit3'
492 c stop
493 c endif
494 
495  return
496  end
497 c-----------------------------------------------------------------------
498  subroutine norcmp2(wvx,wvy,wvz,e,f)
499  include 'SIZE'
500  include 'GEOM'
501  include 'INPUT'
502 
503 
504  real wvx(lx1,ly1,lz1),wvy(lx1,ly1,lz1),wvz(lx1,ly1,lz1)
505 
506  integer e,f
507 
508  common /scruz/ r1(lx1,ly1,lz1),r2(lx1,ly1,lz1),r3(lx1,ly1,lz1)
509 
510  call facind(i0,i1,j0,j1,k0,k1,lx1,ly1,lz1,f)
511 
512  l=0
513  do k=k0,k1
514  do j=j0,j1
515  do i=i0,i1
516  l=l+1
517  scale=wvx(i,j,k)*unx(l,1,f,e)
518  $ +wvy(i,j,k)*uny(l,1,f,e)
519  $ +wvz(i,j,k)*unz(l,1,f,e)
520  wvx(i,j,k) = scale*unx(l,1,f,e)
521  wvy(i,j,k) = scale*uny(l,1,f,e)
522  wvz(i,j,k) = scale*unz(l,1,f,e)
523  enddo
524  enddo
525  enddo
526 
527 c wvxm = vlamax(wvx,lx1*ly1)
528 c write(6,*) f,e,wvxm,' w-max'
529 
530  return
531  end
532 c-----------------------------------------------------------------------
533  subroutine norcmp (wt1,wt2,wt3,rnx,rny,rnz,ifc)
534 C
535  include 'SIZE'
536  COMMON /scruz/ r1(lx1,ly1,lz1),r2(lx1,ly1,lz1),r3(lx1,ly1,lz1)
537 C
538  dimension wt1(lx1,ly1,lz1),wt2(lx1,ly1,lz1),wt3(lx1,ly1,lz1)
539  $ , rnx(lx1,ly1,lz1),rny(lx1,ly1,lz1),rnz(lx1,ly1,lz1)
540 C
541  nxyz1 = lx1*ly1*lz1
542 C
543  CALL copy (r1,wt1,nxyz1)
544  CALL copy (r2,wt2,nxyz1)
545  IF (ldim.EQ.3) CALL copy (r3,wt3,nxyz1)
546  CALL facind2 (js1,jf1,jskip1,js2,jf2,jskip2,ifc)
547 C
548  IF (ldim.EQ.2) THEN
549  DO 200 j2=js2,jf2,jskip2
550  DO 200 j1=js1,jf1,jskip1
551  wn = r1(j1,j2,1)*rnx(j1,j2,1) +
552  $ r2(j1,j2,1)*rny(j1,j2,1)
553  wt1(j1,j2,1) = wn *rnx(j1,j2,1)
554  wt2(j1,j2,1) = wn *rny(j1,j2,1)
555  200 CONTINUE
556  ELSE
557  DO 300 j2=js2,jf2,jskip2
558  DO 300 j1=js1,jf1,jskip1
559  wn = r1(j1,j2,1)*rnx(j1,j2,1) +
560  $ r2(j1,j2,1)*rny(j1,j2,1) +
561  $ r3(j1,j2,1)*rnz(j1,j2,1)
562  wt1(j1,j2,1) = wn *rnx(j1,j2,1)
563  wt2(j1,j2,1) = wn *rny(j1,j2,1)
564  wt3(j1,j2,1) = wn *rnz(j1,j2,1)
565  300 CONTINUE
566  ENDIF
567 C
568  return
569  end
570 c-----------------------------------------------------------------------
571  subroutine facemv (wt1,wt2,wt3,rnx,rny,rnz,smt,ifc)
572 C
573  include 'SIZE'
574  COMMON /scruz/ r1(lx1,ly1,lz1),r2(lx1,ly1,lz1),r3(lx1,ly1,lz1)
575 C
576  dimension wt1(lx1,ly1,lz1),wt2(lx1,ly1,lz1),wt3(lx1,ly1,lz1)
577  $ , rnx(lx1,ly1,lz1),rny(lx1,ly1,lz1),rnz(lx1,ly1,lz1)
578  $ , smt(lx1,ly1,lz1)
579 C
580  nxyz1 = lx1*ly1*lz1
581 C
582  CALL copy (r1,wt1,nxyz1)
583  CALL copy (r2,wt2,nxyz1)
584  IF (ldim.EQ.3) CALL copy (r3,wt3,nxyz1)
585  CALL facind2 (js1,jf1,jskip1,js2,jf2,jskip2,ifc)
586 C
587  IF (ldim.EQ.2) THEN
588  DO 200 j2=js2,jf2,jskip2
589  DO 200 j1=js1,jf1,jskip1
590  wn = ( r1(j1,j2,1)*rnx(j1,j2,1) +
591  $ r2(j1,j2,1)*rny(j1,j2,1) ) / smt(j1,j2,1)
592  wt1(j1,j2,1) = wn *rnx(j1,j2,1)
593  wt2(j1,j2,1) = wn *rny(j1,j2,1)
594  200 CONTINUE
595  ELSE
596  DO 300 j2=js2,jf2,jskip2
597  DO 300 j1=js1,jf1,jskip1
598  wn = ( r1(j1,j2,1)*rnx(j1,j2,1) +
599  $ r2(j1,j2,1)*rny(j1,j2,1) +
600  $ r3(j1,j2,1)*rnz(j1,j2,1) ) / smt(j1,j2,1)
601  wt1(j1,j2,1) = wn *rnx(j1,j2,1)
602  wt2(j1,j2,1) = wn *rny(j1,j2,1)
603  wt3(j1,j2,1) = wn *rnz(j1,j2,1)
604  300 CONTINUE
605  ENDIF
606 C
607  return
608  end
609 c-----------------------------------------------------------------------
610  subroutine faczqn (wt1,wt2,wt3,ifc,iel)
611 C
612  include 'SIZE'
613  include 'GEOM'
614  include 'TOPOL'
615  COMMON /scruz/ r1(lx1,ly1,lz1),r2(lx1,ly1,lz1),r3(lx1,ly1,lz1)
616 C
617  dimension wt1(lx1,ly1,lz1),wt2(lx1,ly1,lz1),wt3(lx1,ly1,lz1)
618 C
619  nxyz1 = lx1*ly1*lz1
620  CALL copy (r1,wt1,nxyz1)
621  CALL copy (r2,wt2,nxyz1)
622  IF (ldim.EQ.3) CALL copy (r3,wt3,nxyz1)
623 C
624  CALL facind2 (js1,jf1,jskip1,js2,jf2,jskip2,ifc)
625  i = 0
626 C
627  IF (ldim.EQ.2) THEN
628  DO 200 j2=js2,jf2,jskip2
629  DO 200 j1=js1,jf1,jskip1
630  i = i+1
631  w1 = r1(j1,j2,1)*t1x(i,1,ifc,iel) +
632  $ r2(j1,j2,1)*t1y(i,1,ifc,iel)
633  wt1(j1,j2,1) = w1 *t1x(i,1,ifc,iel)
634  wt2(j1,j2,1) = w1 *t1y(i,1,ifc,iel)
635  200 CONTINUE
636  ELSE
637  DO 300 j2=js2,jf2,jskip2
638  DO 300 j1=js1,jf1,jskip1
639  i = i+1
640  w1 = r1(j1,j2,1)*t1x(i,1,ifc,iel) +
641  $ r2(j1,j2,1)*t1y(i,1,ifc,iel) +
642  $ r3(j1,j2,1)*t1z(i,1,ifc,iel)
643  wt1(j1,j2,1) = w1 *t1x(i,1,ifc,iel)
644  wt2(j1,j2,1) = w1 *t1y(i,1,ifc,iel)
645  wt3(j1,j2,1) = w1 *t1z(i,1,ifc,iel)
646  300 CONTINUE
647  ENDIF
648 C
649  return
650  end
651 c-----------------------------------------------------------------------
652  subroutine facsmt (smt,ifc)
653 C
654  include 'SIZE'
655  dimension smt(lx1,ly1,lz1)
656 C
657  CALL facind2 (js1,jf1,jskip1,js2,jf2,jskip2,ifc)
658 C
659  DO 100 j2=js2,jf2,jskip2
660  DO 100 j1=js1,jf1,jskip1
661  smt(j1,j2,1)=smt(j1,j2,1) + 1.0
662  100 CONTINUE
663 C
664  return
665  end
666 c-----------------------------------------------------------------------
667  subroutine cqnet (qni,ta,nel)
668 C
669  include 'SIZE'
670  include 'INPUT'
671  include 'SOLN'
672  include 'TSTEP'
673  COMMON /scrvh/ h1(lx1,ly1,lz1,lelt)
674  $ , h2(lx1,ly1,lz1,lelt)
675 C
676  dimension qni(lx1,ly1,lz1,1)
677  $ , ta(lx1,ly1,lz1,1)
678 C
679  intloc = -1
680  imshl = 2
681  ntot1 = lx1*ly1*lz1*nel
682 C
683  CALL sethlm (h1,h2,intloc)
684  CALL axhelm (ta,t(1,1,1,1,ifield-1),h1,h2,imshl,1)
685  CALL sub3 (qni,ta,bq(1,1,1,1,ifield-1),ntot1)
686  CALL dssum (qni,lx1,ly1,lz1)
687 C
688  return
689  end
690 c-----------------------------------------------------------------------
691  subroutine facemt (w1,w2,w3,rnx,rny,rnz,qni,dsa,smt,rhola,ifc)
692 C
693  include 'SIZE'
694  include 'GEOM'
695 C
696  dimension w1(lx1,ly1,lz1)
697  $ , w2(lx1,ly1,lz1)
698  $ , w3(lx1,ly1,lz1)
699  $ , rnx(lx1,ly1,lz1)
700  $ , rny(lx1,ly1,lz1)
701  $ , rnz(lx1,ly1,lz1)
702  $ , qni(lx1,ly1,lz1)
703  $ , dsa(lx1,ly1,lz1)
704  $ , smt(lx1,ly1,lz1)
705 C
706  CALL facind2 (js1,jf1,jskip1,js2,jf2,jskip2,ifc)
707 C
708  IF (ldim.EQ.2) THEN
709  DO 200 j2=js2,jf2,jskip2
710  DO 200 j1=js1,jf1,jskip1
711  aa = qni(j1,j2,1) / ( dsa(j1,j2,1)*smt(j1,j2,1)*rhola )
712  w1(j1,j2,1) = rnx(j1,j2,1) * aa
713  w2(j1,j2,1) = rny(j1,j2,1) * aa
714  200 CONTINUE
715  ELSE
716  DO 300 j2=js2,jf2,jskip2
717  DO 300 j1=js1,jf1,jskip1
718  aa = qni(j1,j2,1) / ( dsa(j1,j2,1)*smt(j1,j2,1)*rhola )
719  w1(j1,j2,1) = rnx(j1,j2,1) * aa
720  w2(j1,j2,1) = rny(j1,j2,1) * aa
721  w3(j1,j2,1) = rnz(j1,j2,1) * aa
722  300 CONTINUE
723  ENDIF
724 C
725  return
726  end
727 c-----------------------------------------------------------------------
728  subroutine elasolv (nel)
729 C
730 C Elastostatic solver for mesh deformation
731 C
732  include 'SIZE'
733  include 'GEOM'
734  include 'INPUT'
735  include 'MVGEOM'
736  include 'SOLN'
737  include 'TSTEP'
738 C
739  COMMON /scrns/ dw1(lx1,ly1,lz1,lelt)
740  $ , dw2(lx1,ly1,lz1,lelt)
741  $ , dw3(lx1,ly1,lz1,lelt)
742  $ , aw1(lx1,ly1,lz1,lelt)
743  $ , aw2(lx1,ly1,lz1,lelt)
744  $ , aw3(lx1,ly1,lz1,lelt)
745  COMMON /scrvh/ h1(lx1,ly1,lz1,lelt)
746  $ , h2(lx1,ly1,lz1,lelt)
747  common /scruz/ prt(lx1,ly1,lz1,lelt)
748  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
749  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
750 
751  if (ifusermv) return ! Compute wx,wy,wz in userchk.
752 c
753  ntot1 = lx1*ly1*lz1*nel
754  maxit = 1000
755  matmod = -1
756  ifh2 = .false.
757  ifsolv = .true.
758  imsolv = 0
759  vnu = 0.0
760  vnu = param(47)
761  if (vnu.eq.0) vnu = 0.4
762  vnu = max(0.00,vnu)
763  vnu = min(0.499,vnu)
764 
765 C Set up elastic material constants
766 
767  ce = 1./(1. + vnu)
768  c2 = vnu * ce / (1. - 2.*vnu)
769  c3 = 0.5 * ce
770  CALL cfill (h1,c2,ntot1)
771  CALL cfill (h2,c3,ntot1)
772 
773 C Solve for interior mesh velocities
774 
775  CALL meshtol (aw1,tolmsh,nel,imsolv)
776  IF (imsolv.EQ.1) return
777 
778  CALL axhmsf (aw1,aw2,aw3,wx,wy,wz,h1,h2,matmod)
779 
780 c if (istep.eq.2) then
781 c call outpost(wx,wy,wz,h1,h2,' ')
782 c call outpost(aw1,aw2,aw3,h1,h2,' ')
783 c write(6,*) 'quit elas 1'
784 c stop
785 c endif
786 
787  CALL chsign (aw1,ntot1)
788  CALL chsign (aw2,ntot1)
789  IF (ldim.EQ.3) CALL chsign (aw3,ntot1)
790  CALL hmhzsf ('NOMG',dw1,dw2,dw3,aw1,aw2,aw3,h1,h2,
791  $ w1mask,w2mask,w3mask,wmult,tolmsh,
792  $ maxit,matmod)
793 C
794 C Update mesh velocities
795 C
796  CALL add2 (wx,dw1,ntot1)
797  CALL add2 (wy,dw2,ntot1)
798  IF (ldim.EQ.3) CALL add2 (wz,dw3,ntot1)
799 
800  ifldt = ifield
801  ifield=1
802  if (ifheat) ifield=2
803  call dsavg(wx)
804  call dsavg(wy)
805  call dsavg(wz)
806  ifield = ifldt
807 
808 c if (istep.gt.1) then
809 c ifldx = ifield
810 c ifield = 1
811 c call incomprn (wx,wy,wz,prt) ! project U onto div-free space
812 c ifield = ifldx
813 c return
814 c endif
815 
816  return
817  end
818 c-----------------------------------------------------------------------
819  subroutine meshtol (ta,tolmsh,nel,imsolv)
820 C
821  include 'SIZE'
822  include 'EIGEN'
823  include 'MVGEOM'
824  include 'TSTEP'
825  dimension ta(lx1,ly1,lz1,1)
826 C
827  ntot1 = lx1*ly1*lz1*nel
828  tolab = tolrel
829 C
830  delta = 1.0e-9
831  x = 1.0 + delta
832  y = 1.0
833  diff = abs(x - y)
834  IF (diff .EQ. 0.0) eps = 1.0e-05
835  IF (diff .GT. 0.0) eps = 1.0e-12
836 C
837  CALL opdot (ta,wx,wy,wz,wx,wy,wz,ntot1)
838 C
839  wdot = glmax(ta,ntot1)
840  wmax = sqrt(wdot)
841  IF (wmax .LT. eps) THEN
842  imsolv = 1
843  return
844  ENDIF
845 C
846  tolmsh = tolab * wmax * sqrt(eigaa)
847 C
848  return
849  end
850 c-----------------------------------------------------------------------
851  subroutine updxyz (nel)
852 C
853  include 'SIZE'
854  include 'TSTEP'
855  include 'MVGEOM'
856  include 'GEOM'
857  include 'INPUT'
858  COMMON /scrsf/ ux(lx1,ly1,lz1,lelt)
859  $ , uy(lx1,ly1,lz1,lelt)
860  $ , uz(lx1,ly1,lz1,lelt)
861  dimension abm(3)
862 C
863  ntot1 = lx1*ly1*lz1*nel
864 C
865  DO 10 i=1,nbd
866  10 abm(i) = dt*abmsh(i)
867 C
868  IF (istep.EQ.0) THEN
869  CALL copy (ux,wx,ntot1)
870  CALL copy (uy,wy,ntot1)
871  IF (ldim.EQ.3) CALL copy (uz,wz,ntot1)
872  else
873  if (ifrich) then
874  call cmult2(ux,wx,dt,ntot1)
875  call cmult2(uy,wy,dt,ntot1)
876  if (ldim.eq.3) call cmult2(uz,wz,dt,ntot1)
877  else
878  CALL cmult2 (ux,wx,abm(1),ntot1)
879  CALL cmult2 (uy,wy,abm(1),ntot1)
880  IF (ldim.EQ.3) CALL cmult2 (uz,wz,abm(1),ntot1)
881  DO 100 ilag=2,nbd
882  CALL add2s2 (ux,wxlag(1,1,1,1,ilag-1),abm(ilag),ntot1)
883  CALL add2s2 (uy,wylag(1,1,1,1,ilag-1),abm(ilag),ntot1)
884  IF (ldim.EQ.3)
885  $ CALL add2s2 (uz,wzlag(1,1,1,1,ilag-1),abm(ilag),ntot1)
886  100 CONTINUE
887  endif
888  ENDIF
889 C
890  CALL add2 (xm1,ux,ntot1)
891  CALL add2 (ym1,uy,ntot1)
892  IF (ldim.EQ.3) CALL add2 (zm1,uz,ntot1)
893 C
894  return
895  end
896 c-----------------------------------------------------------------------
897  subroutine lagmshv (nel)
898 C-----------------------------------------------------------------------
899 C
900 C Keep old mesh velocity
901 C
902 C-----------------------------------------------------------------------
903  include 'SIZE'
904  include 'INPUT'
905  include 'MVGEOM'
906  include 'TSTEP'
907 C
908  ntot1 = lx1*ly1*lz1*nel
909 C
910  DO 100 ilag=nbdinp-1,2,-1
911  CALL copy (wxlag(1,1,1,1,ilag),wxlag(1,1,1,1,ilag-1),ntot1)
912  CALL copy (wylag(1,1,1,1,ilag),wylag(1,1,1,1,ilag-1),ntot1)
913  IF (ldim.EQ.3)
914  $ CALL copy (wzlag(1,1,1,1,ilag),wzlag(1,1,1,1,ilag-1),ntot1)
915  100 CONTINUE
916 C
917  CALL copy (wxlag(1,1,1,1,1),wx,ntot1)
918  CALL copy (wylag(1,1,1,1,1),wy,ntot1)
919  IF (ldim.EQ.3)
920  $ CALL copy (wzlag(1,1,1,1,1),wz,ntot1)
921 C
922  return
923  end
924 c-----------------------------------------------------------------------
925  subroutine facec3 (a1,a2,a3,b1,b2,b3,ifc)
926 C
927 C Copy the face (IFC) of B1,B2,B3 to A1,A2,A3.
928 C IFACE is the input in the pre-processor ordering scheme.
929 C
930  include 'SIZE'
931  dimension a1(lx1,ly1,lz1)
932  $ , a2(lx1,ly1,lz1)
933  $ , a3(lx1,ly1,lz1)
934  $ , b1(lx1,ly1,lz1)
935  $ , b2(lx1,ly1,lz1)
936  $ , b3(lx1,ly1,lz1)
937 C
938  CALL facind2 (js1,jf1,jskip1,js2,jf2,jskip2,ifc)
939 C
940  DO 100 j2=js2,jf2,jskip2
941  DO 100 j1=js1,jf1,jskip1
942  a1(j1,j2,1)=b1(j1,j2,1)
943  a2(j1,j2,1)=b2(j1,j2,1)
944  a3(j1,j2,1)=b3(j1,j2,1)
945  100 CONTINUE
946  return
947  end
948 c-----------------------------------------------------------------------
949  subroutine ptbgeom
950 C-----------------------------------------------------------------------
951 C
952 C Subroutine to impose perturbation to geometry before solution
953 C for moving boundary problems
954 C
955 C-----------------------------------------------------------------------
956  include 'SIZE'
957  include 'GEOM'
958  include 'MVGEOM'
959  include 'SOLN'
960  include 'TSTEP'
961  include 'INPUT'
962  COMMON /scruz/ xm3(lx3,ly3,lz3,lelt)
963  $ , ym3(lx3,ly3,lz3,lelt)
964  $ , zm3(lx3,ly3,lz3,lelt)
965 C
966  IF (istep .EQ. 0) return
967  ifield = 0
968  nel = nelfld(ifield)
969  ntot1 = lx1*ly1*lz1*nel
970  imesh = 1
971  IF ( iftmsh(ifield) ) imesh = 2
972 C
973  CALL ibdgeom (nel)
974  CALL elasolv (nel)
975  CALL updxyz (nel)
976  CALL geom1 (xm3,ym3,zm3)
977  CALL geom2
978  CALL updmsys (0)
979  CALL volume
980  CALL setinvm
981  CALL lagmass
982 C
983  CALL rzero (wx,ntot1)
984  CALL rzero (wy,ntot1)
985  IF (ldim.EQ.3) CALL rzero (wz,ntot1)
986 C
987  return
988  end
989 c-----------------------------------------------------------------------
990  subroutine ibdgeom (nel)
991 C
992 C Routine to evaluate mesh velocities at all moving boundaries
993 C
994  include 'SIZE'
995  include 'GEOM'
996  include 'INPUT'
997  include 'PARALLEL'
998  include 'MVGEOM'
999  include 'TSTEP'
1000 C
1001  CHARACTER CB*1
1002 C
1003  nface = 2*ldim
1004  ntot1 = lx1*ly1*lz1*nel
1005 C
1006  CALL rzero (wx,ntot1)
1007  CALL rzero (wy,ntot1)
1008  CALL rzero (wz,ntot1)
1009 C
1010  DO 1000 isweep=1,2
1011 C
1012  ifld = 0
1013  DO 110 iel=1,nel
1014  DO 110 ifc=1,nface
1015  ieg = lglel(iel)
1016  cb = cbc(ifc,iel,ifld)
1017  IF (cb.EQ.'M' .OR. cb.EQ.'m') THEN
1018  CALL facind (kx1,kx2,ky1,ky2,kz1,kz2,lx1,ly1,lz1,ifc)
1019  DO 140 iz=kz1,kz2
1020  DO 140 iy=ky1,ky2
1021  DO 140 ix=kx1,kx2
1022  CALL inigeom (wx(ix,iy,iz,iel),wy(ix,iy,iz,iel),
1023  $ wz(ix,iy,iz,iel),xm1(ix,iy,iz,iel),
1024  $ ym1(ix,iy,iz,iel),zm1(ix,iy,iz,iel),
1025  $ ifc,ieg)
1026  140 CONTINUE
1027  ENDIF
1028  110 CONTINUE
1029 C
1030  IF (isweep.EQ.1) THEN
1031  CALL dsop (wx,'MXA',lx1,ly1,lz1)
1032  CALL dsop (wy,'MXA',lx1,ly1,lz1)
1033  IF (ldim.EQ.3) CALL dsop (wz,'MXA',lx1,ly1,lz1)
1034  ELSE
1035  CALL dsop (wx,'MNA',lx1,ly1,lz1)
1036  CALL dsop (wy,'MNA',lx1,ly1,lz1)
1037  IF (ldim.EQ.3) CALL dsop (wz,'MNA',lx1,ly1,lz1)
1038  ENDIF
1039 C
1040  1000 CONTINUE
1041 C
1042  return
1043  end
1044 c-----------------------------------------------------------------------
1045  subroutine inigeom (ux,uy,uz,x,y,z,iside,iel)
1046 C
1047  include 'SIZE'
1048  include 'TSTEP'
1049 C
1050  ux = 0.0
1051  uy = 0.0
1052  uz = 0.0
1053 C
1054  return
1055  end
1056 c-----------------------------------------------------------------------
1057  subroutine quickmv
1058  include 'SIZE'
1059  include 'TOTAL'
1060  include 'ZPER'
1061 c
1062  if (if3d) then
1063  call quickmv3d
1064  else
1065  call quickmv2d
1066  endif
1067  return
1068  end
1069 c-----------------------------------------------------------------------
1070  subroutine quickmv2d
1071  include 'SIZE'
1072  include 'TOTAL'
1073  include 'ZPER'
1074 c
1075  integer e,ex,ey,ez,eg
1076  common /surfa/ zsurf(lx1,lz1,lelx,lely)
1077  $ , wsurf(lx1,lz1,lelx,lely)
1078  real nxs,nys,nzs
1079 c
1080  icount = 0
1081  do ex=1,nelx
1082  do ix=1,lx1
1083  zsurf(ix,1,ex,1) = -1.e20
1084  wsurf(ix,1,ex,1) = -1.e20
1085  ey=nely
1086  eg = ex + nelx*(ey-1)
1087  mid = gllnid(eg)
1088  e = gllel(eg)
1089  if (mid.eq.nid) then
1090  zsurf(ix,1,ex,1) = ym1(ix,ly1,1,e)
1091  vxs = vx(ix,ly1,1,e)
1092  vys = vy(ix,ly1,1,e)
1093  nxs = unx(ix,1,3,e) ! Face 3 is on top in 2D
1094  nys = uny(ix,1,3,e)
1095  gamma_s = (vxs*nxs + vys*nys)/(nys)
1096  wsurf(ix,1,ex,1) = gamma_s ! vertical component of wsurf
1097  endif
1098  zsurf(ix,1,ex,1) = glmax(zsurf(ix,1,ex,1),1)
1099  wsurf(ix,1,ex,1) = glmax(wsurf(ix,1,ex,1),1)
1100  icount = icount+1
1101 
1102 c write(6,6) ex,e,ix,xm1(ix,ly1,1,e),ym1(ix,ly1,1,e)
1103 c $ ,vxs,vys,nxs,nys,gamma_s,wsurf(ix,1,ex,1),zsurf(ix,1,ex,1)
1104 c 6 format(3i3,1p9e12.4,' srf')
1105 
1106  enddo
1107  enddo
1108  zmin = glmin(ym1,lx1*ly1*lz1*nelv)
1109 c
1110  do ex=1,nelx
1111  do ix=1,lx1
1112  do ey=1,nely
1113  eg = ex + nelx*(ey-1)
1114  mid = gllnid(eg)
1115  e = gllel(eg)
1116  if (mid.eq.nid) then
1117  do iy=1,ly1
1118  wy(ix,iy,1,e) = wsurf(ix,1,ex,1)
1119  $ * (ym1(ix,iy,1,e)-zmin)/(zsurf(ix,1,ex,1)-zmin)
1120  enddo
1121  endif
1122  enddo
1123  enddo
1124  enddo
1125 c
1126  n = lx1*ly1*lz1*nelv
1127  call rzero(wx,n)
1128  call dsavg(wy)
1129 
1130 c call opcopy(vx,vy,vz,wx,wy,wz)
1131 c call outpost(vx,vy,vz,pr,t,' ')
1132 c call exitt
1133 
1134  return
1135  end
1136 c-----------------------------------------------------------------------
1137  subroutine quickmv3d
1138  include 'SIZE'
1139  include 'TOTAL'
1140  include 'ZPER'
1141 c
1142  integer e,ex,ey,ez,eg
1143  common /surfa/ zsurf(lx1,lz1,lelx,lely)
1144  $ , wsurf(lx1,lz1,lelx,lely)
1145  real nxs,nys,nzs
1146 c
1147  icount = 0
1148  do ey=1,nely
1149  do ex=1,nelx
1150  do iy=1,ly1
1151  do ix=1,lx1
1152  zsurf(ix,iy,ex,ey) = -1.e20
1153  wsurf(ix,iy,ex,ey) = -1.e20
1154  ez = nelz
1155  eg = ex + nelx*(ey-1) + nelx*nely*(ez-1)
1156  mid = gllnid(eg)
1157  e = gllel(eg)
1158  if (mid.eq.nid) then
1159  zsurf(ix,iy,ex,ey) = zm1(ix,iy,lz1,e)
1160  vxs = vx(ix,iy,lz1,e)
1161  vys = vy(ix,iy,lz1,e)
1162  vzs = vz(ix,iy,lz1,e)
1163  nxs = unx(ix,iy,6,e) ! Face 6 is on top in 3D
1164  nys = uny(ix,iy,6,e)
1165  nzs = unz(ix,iy,6,e)
1166  gamma_s = (vxs*nxs+vys*nys+vzs*nzs)/(nzs)
1167  wsurf(ix,iy,ex,ey) = gamma_s ! vertical component of wsurf
1168  endif
1169  zsurf(ix,iy,ex,ey) = glmax(zsurf(ix,iy,ex,ey),1)
1170  wsurf(ix,iy,ex,ey) = glmax(wsurf(ix,iy,ex,ey),1)
1171  icount = icount+1
1172  enddo
1173  enddo
1174  enddo
1175  enddo
1176 
1177  n = lx1*ly1*lz1*nelv
1178  zmin = glmin(zm1,n)
1179 c
1180  do ey=1,nely
1181  do ex=1,nelx
1182  do iy=1,ly1
1183  do ix=1,lx1
1184  do ez=1,nelz
1185  eg = ex + nelx*(ey-1) + nelx*nely*(ez-1)
1186  mid = gllnid(eg)
1187  e = gllel(eg)
1188  if (mid.eq.nid) then
1189  do iz=1,lz1
1190  wz(ix,iy,iz,e) = wsurf(ix,iy,ex,ey)
1191  $ * (zm1(ix,iy,iz,e)-zmin)/(zsurf(ix,iy,ex,ey)-zmin)
1192  enddo
1193  endif
1194  enddo
1195  enddo
1196  enddo
1197  enddo
1198  enddo
1199 c
1200  n = lx1*ly1*lz1*nelv
1201  call rzero(wx,n)
1202  call rzero(wy,n)
1203  call dsavg(wz)
1204 c
1205  return
1206  end
1207 c-----------------------------------------------------------------------
unitvec
subroutine unitvec(X, Y, Z, N)
Definition: bdry.f:1831
chknord
subroutine chknord(IFALGN, IFNORX, IFNORY, IFNORZ, IFC, IEL)
Definition: bdry.f:185
rzero3
subroutine rzero3(A, B, C, N)
Definition: bdry.f:1821
facind2
subroutine facind2(JS1, JF1, JSKIP1, JS2, JF2, JSKIP2, IFC)
Definition: bdry.f:2098
setinvm
subroutine setinvm
Definition: coef.f:1334
geom1
subroutine geom1(xm3, ym3, zm3)
Definition: coef.f:361
volume
subroutine volume
Definition: coef.f:1003
geom2
subroutine geom2
Definition: coef.f:823
lagmass
subroutine lagmass
Definition: coef.f:1315
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
scale
subroutine scale(xyzl, nl)
Definition: connect2.f:602
gllel
integer function gllel(ieg)
Definition: dprocmap.f:183
gllnid
integer function gllnid(ieg)
Definition: dprocmap.f:161
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
dsavg
subroutine dsavg(u)
Definition: ic.f:1878
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
glmin
function glmin(a, n)
Definition: math.f:973
invcol2
subroutine invcol2(a, b, n)
Definition: math.f:73
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
add2
subroutine add2(a, b, n)
Definition: math.f:622
col2
subroutine col2(a, b, n)
Definition: math.f:564
add2s2
subroutine add2s2(a, b, c1, n)
Definition: math.f:690
copy
subroutine copy(a, b, n)
Definition: math.f:260
addcol4
subroutine addcol4(a, b, c, d, n)
Definition: math.f:142
iglmax
function iglmax(a, n)
Definition: math.f:913
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
chsign
subroutine chsign(a, n)
Definition: math.f:305
quickmv2d
subroutine quickmv2d
Definition: mvmesh.f:1071
cbcmesh
subroutine cbcmesh
Definition: mvmesh.f:3
facsmt
subroutine facsmt(smt, ifc)
Definition: mvmesh.f:653
admeshv
subroutine admeshv
Definition: mvmesh.f:82
ptbgeom
subroutine ptbgeom
Definition: mvmesh.f:950
quickmv3d
subroutine quickmv3d
Definition: mvmesh.f:1138
facec3
subroutine facec3(a1, a2, a3, b1, b2, b3, ifc)
Definition: mvmesh.f:926
facemv
subroutine facemv(wt1, wt2, wt3, rnx, rny, rnz, smt, ifc)
Definition: mvmesh.f:572
cqnet
subroutine cqnet(qni, ta, nel)
Definition: mvmesh.f:668
facemt
subroutine facemt(w1, w2, w3, rnx, rny, rnz, qni, dsa, smt, rhola, ifc)
Definition: mvmesh.f:692
faczqn
subroutine faczqn(wt1, wt2, wt3, ifc, iel)
Definition: mvmesh.f:611
updxyz
subroutine updxyz(nel)
Definition: mvmesh.f:852
meshtol
subroutine meshtol(ta, tolmsh, nel, imsolv)
Definition: mvmesh.f:820
divws
subroutine divws(fms, sfv, phi, nel, idir)
Definition: mvmesh.f:132
mvbdry
subroutine mvbdry(nel)
Definition: mvmesh.f:250
norcmp2
subroutine norcmp2(wvx, wvy, wvz, e, f)
Definition: mvmesh.f:499
axifms
subroutine axifms(fms, sfv, phi, nel, idir)
Definition: mvmesh.f:180
ibdgeom
subroutine ibdgeom(nel)
Definition: mvmesh.f:991
norcmp
subroutine norcmp(wt1, wt2, wt3, rnx, rny, rnz, ifc)
Definition: mvmesh.f:534
updcoor
subroutine updcoor
Definition: mvmesh.f:226
inigeom
subroutine inigeom(ux, uy, uz, x, y, z, iside, iel)
Definition: mvmesh.f:1046
lagmshv
subroutine lagmshv(nel)
Definition: mvmesh.f:898
quickmv
subroutine quickmv
Definition: mvmesh.f:1058
admesht
subroutine admesht
Definition: mvmesh.f:111
elasolv
subroutine elasolv(nel)
Definition: mvmesh.f:729
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
opdssum
subroutine opdssum(a, b, c)
Definition: navier1.f:2582
opcopy
subroutine opcopy(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2474
rnz
Definition: crs_amg.c:1093
axhmsf
subroutine axhmsf(au1, au2, au3, u1, u2, u3, h1, h2, matmod)
Definition: subs1.f:1241
hmhzsf
subroutine hmhzsf(name, u1, u2, u3, r1, r2, r3, h1, h2, rmask1, rmask2, rmask3, rmult, tol, maxit, matmod)
Definition: subs1.f:1096
sethlm
subroutine sethlm(h1, h2, intloc)
Definition: subs1.f:1014
urst
subroutine urst(u, ur, us, ut, nel)
Definition: subs1.f:1490
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341
facexv
subroutine facexv(A1, A2, A3, B1, B2, B3, IFACE1, IOP)
Definition: subs2.f:170
opdot
subroutine opdot(DP, A1, A2, A3, B1, B2, B3, N)
Definition: subs2.f:80
updmsys
subroutine updmsys(IFLD)
Definition: subs2.f:712
facexs
subroutine facexs(A, B, IFACE1, IOP)
Definition: subs2.f:125
rmask
subroutine rmask(R1, R2, R3, NEL)
Definition: subs2.f:1801