KTH framework for Nek5000 toolboxes; testing version  0.0.1
meshsmoother.f
Go to the documentation of this file.
1  subroutine meshsmoother
2  include 'SIZE'
3  include 'TOTAL'
4 
5  parameter(ndfsbc=1) !number of boundary conditions
6  character*3 dfsbc(ndfsbc)
7  save dfsbc
8  data dfsbc /'W '/ !BCs listed here
9 
10  idftyp = 0 !distance function - 0 -> exponential, 1-> tanh
11  alpha = 15. !Input for wall distance function
12  beta = 0.1 !
13 
14  nouter = 50 !total loops around laplacian and optimizer smoothing
15  nlap = 20 !number of laplacian iterations in each loop
16  nopt = 20 !number of optimization iterations in each loop
17 
18  mtyp = 1 !metric type
19 
20  call smoothmesh(mtyp,nouter,nlap,nopt,ndfsbc,dfsbc,idftyp,alpha,
21  $ beta)
22 
23  return
24  end
25 c-----------------------------------------------------------------------
26  subroutine smoothmesh(mtyp,nouter,nlap,nopt,nbc,dcbc,
27  $ idftyp,alpha,beta)
28  include 'SIZE'
29  include 'TOTAL'
30  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
31  parameter(nxyzc=lxc*lyc*lzc,nxyz=lx1*ly1*lz1)
32 
33  real xmc(lxc,lyc,lzc,lelv),
34  $ ymc(lxc,lyc,lzc,lelv),
35  $ zmc(lxc,lyc,lzc,lelv),
36  $ mltc(lxc,lyc,lzc,lelv)
37  common /coarsemesh/ xmc,ymc,zmc
38 
39  real x8(2**ldim,lelv*(2**ldim)),
40  $ y8(2**ldim,lelv*(2**ldim)),
41  $ z8(2**ldim,lelv*(2**ldim)),
42  $ nodmask(2**ldim,lelv*(2**ldim)),
43  $ dis((2**ldim)*lelv*(2**ldim)),
44  $ mlt((2**ldim)*lelv*(2**ldim))
45  common /coarsemesh8/ x8,y8,z8
46 
47  real dx(nxyzc,lelv),dy(nxyzc,lelv),dz(nxyzc,lelv),
48  $ xbp(nxyzc,lelv),ybp(nxyzc,lelv),zbp(nxyzc,lelv)
49  common / msmbackup / dx,dy,dz,xbp,ybp,zbp
50 
51  real dxm(nxyz,lelv),dym(nxyz,lelv),dzm(nxyz,lelv),
52  $ dxn(nxyz,lelv),dyn(nxyz,lelv),dzn(nxyz,lelv)
53 
54  integer opt,gshl,gshlc,lapinv,optinv
55  real alpha,beta,etstart, etend,f1sav,f1,f1pre
56  character*3 dcbc(nbc)
57 c -------------------
58 c -------------------
59 ccc INITIALIZE VARIABLES
60  lapinv = 0
61  optinv = 0
62 
63  if (nid.eq.0.and.loglevel.ge.5)
64  $ write(6,*) 'SMOOTHER-check original mesh'
65  call fix_geom
66 c Copy original mesh from xm1 to dxm
67  call copy(dxm,xm1,lx1*ly1*lz1*nelv)
68  call copy(dym,ym1,lx1*ly1*lz1*nelv)
69  if (ldim.eq.3) call copy(dzm,zm1,lx1*ly1*lz1*nelv)
70 
71 ccc Generate interpolators for going to lx1 = 3
72  call gen_int_lx1_to_3
73 
74 ccc COPY THE ORIGINAL MESH TO dx,dy,dz vectors
75  call copy(dx,xmc,lxc*lyc*lzc*nelv)
76  call copy(dy,ymc,lxc*lyc*lzc*nelv)
77  if (ldim.eq.3) call copy(dz,zmc,lxc*lyc*lzc*nelv)
78 ccc COPY THE ORIGINAL MESH FOR BACKUP IF MESH BECOMES INVERTED
79  call copy(xbp,xmc,lxc*lyc*lzc*nelt)
80  call copy(ybp,ymc,lxc*lyc*lzc*nelt)
81  if (ldim.eq.3) call copy(zbp,zmc,lxc*lyc*lzc*nelt)
82 
83  call xmtox8(xmc,x8)
84  call xmtox8(ymc,y8)
85  if (ldim.eq.3) call xmtox8(zmc,z8)
86 
87 ccc CREATE MASK
88  call genmask(nodmask,mlt,gshl,mltc,gshlc)
89 
90 ccc CONSTRUCT WEIGHT FUNCTION
91  call disfun(dis,idftyp,alpha,beta,dcbc,nbc)
92 
93  etstart=dnekclock()
94 ccc GET INITIAL ENERGY
95  call getglobsum(x8,y8,z8,mlt,gshl,2**ldim,1,f1sav)
96  if (nid.eq.0)
97  $ write(6,'(A,1p1e13.6)') 'SMOOTHER-initial energy',f1sav
98  f1 = f1sav
99 
100 ccc START SMOOTHING HERE
101  do j=1,nouter
102  f1pre = f1
103  if (nid.eq.0.and.loglevel.ge.2)
104  $ write(6,'(A,I5)') 'SMOOTHER-iteration',j
105  mtyp = 1 !if mtyp = 1, jacobian, 2 = l^2, 3 - scale jacobian
106  if (nlap.gt.0) then
107  call fastlap(nlap,nodmask,mlt,gshl,dis,lapinv,mltc,gshlc)
108  if (lapinv.eq.1) nlap = 0
109  if (lapinv.eq.1) call restbackmesh
110  endif
111  if (nopt.gt.0) then
112  call optcg(nopt,nodmask,mlt,gshl,dis,mtyp,optinv,mltc,gshlc)
113  if (optinv.eq.1) call restbackmesh !xbp->xm1 and xbp->x8
114  if (optinv.eq.1) nopt = 0
115  endif
116 
117  call getglobsum(x8,y8,z8,mlt,gshl,2**ldim,mtyp,f1)
118 
119  call xmtox8(xmc,x8)
120  call xmtox8(ymc,y8)
121  if (ldim.eq.3) call xmtox8(zmc,z8)
122 
123  call xmctoxm1
124 
125  if (nid.eq.0) write(6,'(A,I7,1p1e13.5)') 'SMOOTHER-energy',j,f1
126  if (nid.eq.0.and.loglevel.ge.2) write(6,*) 'loop complete',j
127  if (f1.ge.f1pre) goto 5001 !mesh didn't improve since last iteration
128  if (nopt.eq.0.and.nlap.eq.0) goto 5001 !no iterations left to do
129  enddo
130  5001 continue
131 
132  etend=dnekclock()
133 ccc RESTORE BOUNDARY LAYER
134  call restbndrlayer(dx,dy,dz,dis,mltc,gshlc) !dx,dy,dz is now actually smooth-coarse
135  call fix_geom
136 
137 ccc Solve the Laplace's equation to morph the interior mesh to match
138 ccc the actual surface mesh
139  call opsub3(dxn,dyn,dzn,dxm,dym,dzm,xm1,ym1,zm1)
140  call my_mv_mesh(dxn,dyn,dzn)
141  call opadd2(xm1,ym1,zm1,dxn,dyn,dzn)
142 
143  call getglobsum(x8,y8,z8,mlt,gshl,2**ldim,mtyp,f1)
144  if (nid.eq.0) then
145  write(6,'(A,1p1e13.6)') 'SMOOTHER-initial energy ',f1sav
146  write(6,'(A,1p1e13.6)') 'SMOOTHER-final energy ',f1
147  write(6,'(A,f5.2)') 'SMOOTHER-improve % ',100.*(f1sav-f1)/f1sav
148  write(6,'(A,1p1e13.6)') 'SMOOTHER-time taken ',etend-etstart
149  endif
150 
151  call geom_reset(1) ! recompute Jacobians, etc.
152 
153 ccc OUTPUT THE MESH
154 c call gen_rea_full(1) !output the rea for smooth mesh
155 
156  return
157  end
158 c-----------------------------------------------------------------------
159  subroutine gen_int_lx1_to_3 !interpolate mesh from lx1 to lx1=3
160  include 'SIZE'
161  include 'TOTAL'
162  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
163  real dxmc(3,3),dymc(3,3),dzmc(3,3)
164  real dxtmc(3,3),dytmc(3,3),dztmc(3,3)
165  common /coarseders/ dxmc,dymc,dzmc,dxtmc,dytmc,dztmc
166  real ixtmc3(lx1,3),ixmc3(3,lx1)
167  real ixtmf3(3,lx1),ixmf3(lx1,3)
168  common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
169  real zgmc(3),wxmc(3) !points and weights
170  real zgmf(lx1),wxmf(lx1) !points and weights
171  real w(lx1,lx1)
172 
173 c Generate GLL points and weight
174  call zwgll(zgmc,wxmc,3)
175  call zwgll(zgmf,wxmf,lx1)
176 c Generate derivative matrices with transpose operators
177  call dgll(dxmc,dxtmc,zgmc,3,3)
178  call dgll(dymc,dytmc,zgmc,3,3)
179  call rzero(dzmc,3*3)
180  call rzero(dztmc,3*3)
181  if (ldim.eq.3) call dgll(dzmc,dztmc,zgmc,3,3)
182 c Generate interpolation matrices
183  call igllm (ixmc3,ixtmc3,zgmf,zgmc,nx1,3,nx1,3)
184  call igllm (ixmf3,ixtmf3,zgmc,zgmf,3,nx1,3,nx1)
185 
186 c Interpolate mesh to lx1=3
187  call xm1toxmc
188 
189  return
190  end
191 c-----------------------------------------------------------------------
192  subroutine int_fine_to_coarse_2d(u1,u2,n1,n2)
193  include 'SIZE'
194 c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
195 c u1 is size n1xn1 and u2 is size n2xn2
196  real ixtmc3(lx1,3),ixmc3(3,lx1)
197  real ixtmf3(3,lx1),ixmf3(lx1,3)
198  common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
199  real u1(n1,n1),u2(n2,n2)
200  real w(20,20)
201 
202  if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
203  $ in int_fine_to_coarse and related routines '
204  if (lx1.gt.20) call exitt
205 
206  call mxm(ixmc3,n2,u1,n1,w,n1)
207  call mxm(w,n2,ixtmc3,n1,u2,n2)
208 
209  return
210  end
211 c-----------------------------------------------------------------------
212  subroutine int_fine_to_coarse_3d(u1,u2,n1,n2)
213  include 'SIZE'
214 c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
215 c u1 is size n1xn1 and u2 is size n2xn2
216  real ixtmc3(lx1,3),ixmc3(3,lx1)
217  real ixtmf3(3,lx1),ixmf3(lx1,3)
218  common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
219  real u1(n1,n1,n1),u2(n2,n2,n2)
220  real w(20*20*20)
221  real v(20*20*20)
222 
223  if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
224  $ in int_fine_to_coarse and related routines '
225  if (lx1.gt.20) call exitt
226 
227  mm = n1*n1
228  mn = n1*n2
229  nn = n2*n2
230 
231  call mxm(ixmc3,n2,u1,n1,v ,mm)
232  iv=1
233  iw=1
234  do k=1,n1
235  call mxm(v(iv),n2,ixtmc3,n1,w(iw),n2)
236  iv = iv+mn
237  iw = iw+nn
238  enddo
239  call mxm(w,nn,ixtmc3,n1,u2,n2)
240 
241  return
242  end
243 c-----------------------------------------------------------------------
244  subroutine int_coarse_to_fine_2d(u1,u2,n1,n2)
245  include 'SIZE'
246 c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
247 c u1 is size n1xn1 and u2 is size n2xn2
248 c u1 is fine, u2 is coarse
249  real ixtmc3(lx1,3),ixmc3(3,lx1)
250  real ixtmf3(3,lx1),ixmf3(lx1,3)
251  common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
252  real zgmc(3),wxmc(3) !points and weights
253  real zgmf(lx1),wxmf(lx1) !points and weights
254  common /coarswz/ zgmc,wxmc,zgmf,wxmf
255  real u1(n1,n1),u2(n2,n2)
256  real w(20,20)
257 
258  if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
259  $ in int_coarse_to_fine and related routines '
260  if (lx1.gt.20) call exitt
261 
262  call mxm(ixmf3,n1,u2,n2,w,n2)
263  call mxm(w,n1,ixtmf3,n2,u1,n1)
264 
265  return
266  end
267 c-----------------------------------------------------------------------
268  subroutine int_coarse_to_fine_3d(u1,u2,n1,n2)
269  include 'SIZE'
270 c Interpolate u1 to u2 using interpolation matrices in12 and in12^T
271 c u1 is size n1xn1 and u2 is size n2xn2
272  real ixtmc3(lx1,3),ixmc3(3,lx1)
273  real ixtmf3(3,lx1),ixmf3(lx1,3)
274  common /coarseints/ ixmc3,ixtmc3,ixmf3,ixtmf3
275  real u1(n1,n1,n1),u2(n2,n2,n2)
276  real w(20*20*20)
277  real v(20*20*20)
278 
279  if (lx1.gt.20) write(6,*) 'SMOOTHER-increase size of work array
280  $ in int_coarse_to_fine and related routines '
281  if (lx1.gt.20) call exitt
282 
283  mm = n2*n2
284  mn = n2*n1
285  nn = n1*n1
286 
287  call mxm(ixmf3,n1,u2,n2,v ,mm)
288  iv=1
289  iw=1
290  do k=1,n2
291  call mxm(v(iv),n1,ixtmf3,n2,w(iw),n1)
292  iv = iv+mn
293  iw = iw+nn
294  enddo
295  call mxm(w,nn,ixtmf3,n2,u1,n1)
296 
297  return
298  end
299 c-----------------------------------------------------------------------
300  subroutine genbackupmesh !xm1->xbp backup the mesh
301  include 'SIZE'
302  include 'TOTAL'
303 ! backsup the mesh
304  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
305  real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
306  $ zmc(lxc,lyc,lzc,lelv)
307  common /coarsemesh/ xmc,ymc,zmc
308 
309  parameter(lt = lxc*lyc*lzc*lelv)
310  real dx(lt),dy(lt),dz(lt),xbp(lt),ybp(lt),zbp(lt)
311  common / msmbackup / dx,dy,dz,xbp,ybp,zbp
312 
313  call copy(xbp,xmc,lxc*lyc*lzc*nelt)
314  call copy(ybp,ymc,lxc*lyc*lzc*nelt)
315  if (ldim.eq.3) call copy(zbp,zmc,lxc*lyc*lzc*nelt)
316 
317  return
318  end
319 c-----------------------------------------------------------------------
320  subroutine restbackmesh !xbp->xm1 and xbp->x8
321  include 'SIZE'
322  include 'TOTAL'
323  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
324  real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
325  $ zmc(lxc,lyc,lzc,lelv)
326  common /coarsemesh/ xmc,ymc,zmc
327 
328  real x8(2,2,ldim-1,lelv*(2**ldim)),y8(2,2,ldim-1,lelv*(2**ldim))
329  real z8(2,2,ldim-1,lelv*(2**ldim))
330  common /coarsemesh8/ x8,y8,z8
331 
332  parameter(lt = lxc*lyc*lzc*lelv)
333  real dx(lt),dy(lt),dz(lt),xbp(lt),ybp(lt),zbp(lt)
334  common / msmbackup / dx,dy,dz,xbp,ybp,zbp
335 
336  call copy(xmc,xbp,lxc*lyc*lzc*nelt)
337  call copy(ymc,ybp,lxc*lyc*lzc*nelt)
338  if (ldim.eq.3) call copy(zmc,zbp,lxc*lyc*lzc*nelt)
339  call xmtox8(xmc,x8)
340  call xmtox8(ymc,y8)
341  if (ldim.eq.3) call xmtox8(zmc,z8)
342 
343  return
344  end
345 c-----------------------------------------------------------------------
346  subroutine optcg(itmax,nodmask,mlt,gshl,dis,opt,optinv,mltc,gshlc)
347  include 'SIZE'
348  include 'TOTAL'
349 
350  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
351  parameter(nxyzc=lxc*lyc*lzc,nxyz=lx1*ly1*lz1)
352 
353  real xmc(lxc,lyc,lzc,lelv),
354  $ ymc(lxc,lyc,lzc,lelv),
355  $ zmc(lxc,lyc,lzc,lelv),
356  $ mltc(lxc,lyc,lzc,lelv)
357  common /coarsemesh/ xmc,ymc,zmc
358 
359  real x8(2**ldim,lelv*(2**ldim)),
360  $ y8(2**ldim,lelv*(2**ldim)),
361  $ z8(2**ldim,lelv*(2**ldim)),
362  $ nodmask(2**ldim,lelv*(2**ldim)),
363  $ dis(2**ldim,lelv*(2**ldim)),
364  $ mlt(2**ldim,lelv*(2**ldim)),
365  $ dx(2**ldim,lelv*(2**ldim)),
366  $ dy(2**ldim,lelv*(2**ldim)),
367  $ dz(2**ldim,lelv*(2**ldim)),
368  $ dfdx(2**ldim,lelv*(2**ldim),ldim),
369  $ dfdxn(2**ldim,lelv*(2**ldim),ldim),
370  $ sk(2**ldim,lelv*(2**ldim),ldim)
371  common /coarsemesh8/ x8,y8,z8
372 
373  integer iter,gshl,siz,eg,opt,optinv,kerr,gshlc
374  real f2,lambda,num,den,scale2
375  real alpha,bk,num1,den1,wrk1
376  real hval(2**ldim,lelv*(2**ldim))
377 
378  optinv = 0 !initialize to 0
379 
380  if (nid.eq.0.and.loglevel.ge.2)
381  $ write(6,*) 'Optimization loop started'
382 
383  if (opt.eq.1) siz = 2**ldim !for jacobian
384  if (opt.eq.2) siz = 4+8*(ldim-2) !for len
385  n1 = 2**ldim
386  n2 = nelv*(2**ldim)*n1
387 
388  call get_nodscale(hval,x8,y8,z8,gshl)
389 c if (nid.eq.0.and.loglevel.ge.4)
390 c $ write(6,'(A,1p1e13.5)') 'perturbation amount is ',hval
391 
392  call gradf(f2,dfdx,x8,y8,z8,mlt,gshl,siz,opt,hval)
393  do i=1,ldim
394  call copy(sk(1,1,i),dfdx(1,1,i),n2)
395  call cmult(sk(1,1,i),-1.,n2)
396  enddo
397 
398  do iter=1,itmax
399 c do line search at this point
400  call
401  $ dolsalpha(x8,y8,z8,sk,alpha,siz,opt,mlt,gshl,nodmask,iter)
402 
403  if (alpha.eq.0) goto 5002
404 
405  call col4(dx,sk(1,1,1),nodmask,dis,n2)
406  call col4(dy,sk(1,1,2),nodmask,dis,n2)
407  if (ldim.eq.3) call col4(dz,sk(1,1,ldim),nodmask,dis,n2)
408 
409  call add2s2(x8,dx,alpha,n2)
410  call add2s2(y8,dy,alpha,n2)
411  if (ldim.eq.3) call add2s2(z8,dz,alpha,n2)
412 
413 c Gradf at new positions
414  call gradf(f2,dfdxn,x8,y8,z8,mlt,gshl,siz,opt,hval)
415 
416 c get Bk = dfdxn'*dfdxn / (dfdx'*dfdx)
417  num1 = 0.
418  den1 = 0.
419  do k=1,ldim
420  do j=1,nelv*(2**ldim)
421  do i=1,2**ldim
422  num1 = dfdxn(i,j,k)*mlt(i,j)*dfdxn(i,j,k) + num1
423  den1 = dfdx(i,j,k)*mlt(i,j)*dfdx(i,j,k) + den1
424  enddo
425  enddo
426  enddo
427  call gop(num1,wrk1,'+ ',1)
428  call gop(den1,wrk1,'+ ',1)
429  bk = num1/den1
430 
431  do k=1,ldim
432  do j=1,nelv*(2**ldim)
433  do i=1,2**ldim
434  sk(i,j,k) = -dfdxn(i,j,k) + bk*sk(i,j,k)
435  dfdx(i,j,k) = dfdxn(i,j,k)
436  enddo
437  enddo
438  enddo
439 
440  if (mod(iter,5).eq.0.or.iter.eq.itmax) then
441  call x8toxm(xmc,x8)
442  call x8toxm(ymc,y8)
443  if (ldim.eq.3) call x8toxm(zmc,z8)
444  call fixcurs(mltc,gshlc)
445  call xmtox8(xmc,x8)
446  call xmtox8(ymc,y8)
447  if (ldim.eq.3) call xmtox8(zmc,z8)
448  endif
449  enddo
450 
451  5002 continue
452  if (alpha.eq.0) then
453  optinv = 2
454  call x8toxm(xmc,x8)
455  call x8toxm(ymc,y8)
456  if (ldim.eq.3) call x8toxm(zmc,z8)
457  call fixcurs(mltc,gshlc)
458  call xmtox8(xmc,x8)
459  call xmtox8(ymc,y8)
460  if (ldim.eq.3) call xmtox8(zmc,z8)
461  endif
462 
463  call glmapm1chkinv(kerr)
464  if (kerr.gt.0) then
465  optinv = 1 !Terminate smoother
466  if (nid.eq.0) write(6,*) 'Terminating optimizer'
467  else
468  call genbackupmesh
469  endif
470 
471  return
472  end
473 c-----------------------------------------------------------------------
474  subroutine
475  $ dolsalpha(xe,ye,ze,sk,alpha,siz,opt,mlt,gshl,nodmask,iter)
476  include 'SIZE'
477  include 'TOTAL'
478  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
479  real xmc(lxc,lyc,lzc,lelv),
480  $ ymc(lxc,lyc,lzc,lelv),
481  $ zmc(lxc,lyc,lzc,lelv)
482  common /coarsemesh/ xmc,ymc,zmc
483  real xe(2**ldim,lelv*(2**ldim)),
484  $ ye(2**ldim,lelv*(2**ldim)),
485  $ ze(2**ldim,lelv*(2**ldim)),
486  $ dx(2**ldim,lelv*(2**ldim)),
487  $ dy(2**ldim,lelv*(2**ldim)),
488  $ dz(2**ldim,lelv*(2**ldim)),
489  $ xs(2**ldim,lelv*(2**ldim)),
490  $ ys(2**ldim,lelv*(2**ldim)),
491  $ zs(2**ldim,lelv*(2**ldim)),
492  $ nodmask(2**ldim,lelv*(2**ldim)),
493  $ mlt(2**ldim,lelv*(2**ldim)),
494  $ sk((2**ldim),lelv*(2**ldim),ldim),
495  $ jac(2**ldim,lelv*(2**ldim))
496  real alpha,wrk1,f1,f2
497  integer siz,opt,z,gshl,invflag
498  integer n1,n2,chk1,chk2,tstop,maxcount,iter
499 
500  common /lsinfo / lssteps
501  real lssteps(5)
502  integer icalld
503  save icalld
504  data icalld /0/
505  integer idx
506  save idx
507  data idx /0/
508 
509  if (icalld.eq.0) then
510  call rzero(lssteps,5)
511  icalld = 1
512  idx = 0
513  endif
514 
515  if (idx.lt.6) then
516  alpha = 1.0
517  else
518  dumsum = vlsum(lssteps,5)
519  alpha = 2.*dumsum
520  endif
521  alphasav = alpha
522  maxcount = 20
523  chk1 = 0
524  chk2 = 0
525  tstop = 0
526  n1 = 2**ldim
527  n2 = nelv*(2**ldim)*n1
528 
529  call copy(xs(1,1),xe(1,1),n2)
530  call copy(ys(1,1),ye(1,1),n2)
531  if (ldim.eq.3) call copy(zs(1,1),ze(1,1),n2)
532 
533  call getglobsum(xs,ys,zs,mlt,gshl,siz,opt,f1)
534  z = 0
535  do while (tstop.eq.0.and.z.le.maxcount)
536  z = z+1
537  call col3c(dx,sk(1,1,1),nodmask,alpha,n2)
538  call col3c(dy,sk(1,1,2),nodmask,alpha,n2)
539  if (ldim.eq.3) call col3c(dz,sk(1,1,ldim),nodmask,alpha,n2)
540 
541  call add3(xs,xe,dx,n2)
542  call add3(ys,ye,dy,n2)
543  if (ldim.eq.3) call add3(zs,ze,dz,n2)
544 
545  call getglobsum(xs,ys,zs,mlt,gshl,siz,opt,f2)
546 
547  if (f2.lt.f1) then
548  chk1 = 1
549  endif
550  call gop(chk1,wrk1,'m ',1)
551 
552  chk2 = 0
553  if (chk1.eq.1) then
554  call x8toxm(xmc,xs)
555  call x8toxm(ymc,ys)
556  if (ldim.eq.3) call x8toxm(zmc,zs)
557  call glmapm1chkinv(invflag)
558  if (invflag.eq.0) chk2 = 1
559  endif
560  call gop(chk2,wrk1,'m ',1)
561 
562  if (chk1.eq.1.and.chk2.eq.1.) then
563  alphasav = alpha
564  tstop = 1.
565  else
566  chk1 = 0.
567  chk2 = 0.
568  alpha = alpha/(10.)
569  endif
570  enddo
571 
572  if (idx.lt.6) then
573  idx = idx+1
574  lssteps(idx) = alphasav
575  else
576  lssteps(1) = lssteps(2)
577  lssteps(2) = lssteps(3)
578  lssteps(3) = lssteps(4)
579  lssteps(4) = lssteps(5)
580  lssteps(5) = alphasav
581  endif
582 
583  if (tstop.eq.1) then
584  alpha = alphasav
585  if (nid.eq.0.and.loglevel.ge.3) write(6,101) iter,f2
586  101 format(i5,' glob_phi ',1p1e13.6)
587 c write(6,*) z,'number of ls steps'
588  else
589  alpha = 0.
590  if (nid.eq.0.and.loglevel.ge.4)
591  $ write(6,*) 'SMOOTHER-line-search alpha set to 0'
592  endif
593 
594  return
595  end
596 c-----------------------------------------------------------------------
597  subroutine getglobsum(xe,ye,ze,mlt,gshl,siz,opt,fl)
598  include 'SIZE'
599  include 'TOTAL'
600  real xe(2**ldim,lelv*(2**ldim))
601  real ye(2**ldim,lelv*(2**ldim))
602  real ze(2**ldim,lelv*(2**ldim))
603  real mlt(2**ldim,lelv*(2**ldim))
604  integer gshl,siz,e,opt
605  real f1,fl,wrk1
606 
607  fl = 0.
608  do e=1,nelv*(2**ldim)
609  if (opt.eq.1)
610  $ call get_jac(f1,xe(1,e),ye(1,e),ze(1,e),siz,e,1)
611  if (opt.eq.2) call get_len(f1,xe(1,e),ye(1,e),ze(1,e),siz)
612  fl = fl+f1
613  enddo
614  call gop(fl,wrk1,'+ ',1)
615 
616  return
617  end
618 c-----------------------------------------------------------------------
619  subroutine disfun(dis,funtype,alpha,beta,dcbc,nbc)
620  include 'SIZE'
621  include 'TOTAL'
622  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
623  real dd1(lx1*ly1*lz1*lelv),
624  $ dd2(lx1*ly1*lz1*lelv),
625  $ ddd(lx1*ly1*lz1*lelv),
626  $ ddc(lxc*lyc*lzc,lelv),
627  $ dis((2**ldim)*lelv*(2**ldim))
628  integer i,nbc,j,funtype
629  character*3 dcbc(nbc)
630  real dscale,dmax,alpha,beta,dum2
631 
632  if (nid.eq.0.and.loglevel.ge.5)
633  $ write(6,*) 'calculate distance function'
634 
635  call rone (dd1,lx1*ly1*lz1*nelv)
636  call sm_cheap_dist(dd1,1,dcbc(1)) ! Distance function scaling
637 
638  if (nbc.ge.2) then
639  do i=2,nbc
640  call rone (dd2,lx1*ly1*lz1*nelv)
641  call sm_cheap_dist(dd2,1,dcbc(i)) ! Distance function scaling
642  do j=1,lx1*ly1*lz1*nelv
643  dd1(j) = min(dd1(j),dd2(j))
644  enddo
645  enddo
646  endif
647 
648  dmax = glamax(dd1,lx1*ly1*lz1*nelv)
649  dscale = 1./dmax
650  call cmult(dd1,dscale,lx1*ly1*lz1*nelv) !normalized 0 to 1
651 c call outpost(dd1,vy,vz,pr,t,' ')
652 
653  nxyz = lx1*ly1*lz1
654  if (ldim.eq.2) then
655  do i=1,nelv
656  call int_fine_to_coarse_2d(dd1((i-1)*nxyz+1),ddc(1,i),lx1,lxc)
657  enddo
658  else
659  do i=1,nelv
660  call int_fine_to_coarse_3d(dd1((i-1)*nxyz+1),ddc(1,i),lx1,lxc)
661  enddo
662  endif
663  call xmtox8(ddc,dis)
664 c
665  if (funtype.eq.0) then
666  do i=1,(2**ldim)*nelv*(2**ldim)
667  dis(i) = (1-exp(-dis(i)/beta))
668  enddo
669  elseif (funtype.eq.1) then
670  do i=1,(2**ldim)*nelv*(2**ldim)
671  dis(i) = 0.5*(tanh(alpha*(dis(i)-beta))+1)
672  enddo
673  else
674  call exitti('Please set the funtype to 0 or 1$',funtype)
675  endif
676 
677  if (nid.eq.0.and.loglevel.ge.5)
678  $ write(6,'(1p1e13.4,A)') dmax,'max disfun'
679 
680  return
681  end
682 c-----------------------------------------------------------------------
683  subroutine restbndrlayer(dx,dy,dz,dis,mltc,gshlc)
684  include 'SIZE'
685  include 'TOTAL'
686 c dis*smoothmesh + (1-dis)*original mesh
687  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
688  real xmc(lxc,lyc,lzc,lelv),
689  $ ymc(lxc,lyc,lzc,lelv),
690  $ zmc(lxc,lyc,lzc,lelv)
691  common /coarsemesh/ xmc,ymc,zmc
692 
693  real dx(lxc*lyc*lzc,lelv),
694  $ dy(lxc*lyc*lzc,lelv),
695  $ dz(lxc*lyc*lzc,lelv),
696  $ dis((2**ldim)*lelv*(2**ldim)),
697  $ dis2(lxc,lyc,lzc,lelv),
698  $ mltc(lxc,lyc,lzc,lelv)
699 
700  integer gshlc
701 
702  call x8toxm(dis2,dis)
703  dum2 = glamax(dis2,lxc*lyc*lzc*nelv)
704  dum2 = 1./dum2
705  call cmult(dis2,dum2,lxc*lyc*lzc*nelv)
706 
707  call col2(xmc,dis2,lxc*lyc*lzc*nelv) !w*xs
708  call col2(ymc,dis2,lxc*lyc*lzc*nelv)
709  if (ndim.eq.3) call col2(zmc,dis2,lxc*lyc*lzc*nelv)
710 
711  call add2(xmc,dx,lxc*lyc*lzc*nelv) !w*xs+xo
712  call add2(ymc,dy,lxc*lyc*lzc*nelv)
713  if (ndim.eq.3) call add2(zmc,dz,lxc*lyc*lzc*nelv)
714 
715  call col2(dx,dis2,lxc*lyc*lzc*nelv) !w*xo
716  call col2(dy,dis2,lxc*lyc*lzc*nelv)
717  if (ndim.eq.3) call col2(dz,dis2,lxc*lyc*lzc*nelv)
718 
719  call sub2(xmc,dx,lxc*lyc*lzc*nelv) !w*xs+xo-w*xo=w*xs+xo(1-w)
720  call sub2(ymc,dy,lxc*lyc*lzc*nelv) !this is now store in xmc..zmc
721  if (ndim.eq.3) call sub2(zmc,dz,lxc*lyc*lzc*nelv)
722 
723  call sub2(dx,xmc,lxc*lyc*lzc*nelv) !dx=xo-xmc
724  call sub2(dy,ymc,lxc*lyc*lzc*nelv)
725  if (ndim.eq.3) call sub2(dz,zmc,lxc*lyc*lzc*nelv)
726 
727  call fixcurs(mltc,gshlc)
728  call xmctoxm1
729 
730  return
731  end
732 c-----------------------------------------------------------------------
733  subroutine masklayers(nodmask,nlayers)
734  include 'SIZE'
735  include 'TOTAL'
736  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
737  real nodmask(2**ldim,lelv*(2**ldim))
738  real d(lx1*ly1*lz1,lelv)
739  real vcmask(lxc,lyc,lzc,lelv)
740  integer e,f,i,j,nlayers,valm,c
741 
742  character*3 tcbc(5)
743  save tcbc
744  data tcbc /'W ','v ','V ','mv ','MV ' /
745 
746  call x8toxm(v1mask,nodmask)
747 
748  call rone(d,nx1*ny1*nz1*nelv)
749  do e=1,nelv
750  do f=1,2*ldim
751  do c=1,5
752  if(cbc(f,e,1).eq.tcbc(c)) call facev(d,e,f,zero,nx1,ny1,nz1)
753  enddo
754  enddo
755  enddo
756  call dsop(d,'mul',nx1,ny1,nz1)
757 
758  do i=1,nlayers
759  call dsop(d,'mul',nx1,ny1,nz1)
760  do e=1,nelv
761  val = glamin(d(1,e),lx1**ldim)
762  call cmult(d(1,e),val,lx1**ldim)
763  enddo
764  enddo
765 
766  call dsop(d,'mul',nx1,ny1,nz1)
767  call col2(v1mask,d,lx1*ly1*lz1*nelv)
768  do i=1,nelv
769  call int_fine_to_coarse_2d(v1mask(1,1,1,e),vcmask(1,1,1,e),lx1,3)
770  enddo
771  call xmtox8(vcmask,nodmask)
772 
773  return
774  end
775 c-----------------------------------------------------------------------
776  subroutine fixcurs(mltc,gshlc)
777  include 'SIZE'
778  include 'TOTAL'
779 c This routine fixes the curved edges post smoothing
780  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
781  real xmc(lxc,lyc,lzc,lelv),
782  $ ymc(lxc,lyc,lzc,lelv),
783  $ zmc(lxc,lyc,lzc,lelv)
784  common /coarsemesh/ xmc,ymc,zmc
785 
786  common /ctmp0/ zg(3)
787  real w1(lxc*lyc*lzc*lelv),w2(lxc*lyc*lzc*lelv)
788  real vcmask(lxc,lyc,lzc,lelv)
789  integer f,e,nedge,n,nfaces,gshlc
790  integer ke(3,12)
791  real xyz(3,3)
792  real xd(lxc*lyc*lzc),
793  $ yd(lxc*lyc*lzc),
794  $ zd(lxc*lyc*lzc),
795  $ mltc(lxc*lyc*lzc*lelv)
796  real nxv,nyv,nzv
797 
798  save ke
799  data ke / 1, 2, 3, 3, 6, 9, 9, 8, 7, 7, 4, 1
800  $ , 19,20,21, 21,24,27, 27,26,25, 25,22,19
801  $ , 1,10,19, 3,12,21, 9,18,27, 7,16,25 /
802  integer kin(7)
803 
804 c START BY LOOPING OVER EACH ELEMENT AND THEN OVER EACH EDGE
805  nedge = 4 + 8*(ldim-2)
806  nfaces = 2*ldim
807 
808  n = lxc*lyc*lzc*nelv
809  call rone (vcmask,n)
810  do e=1,nelv
811  do f=1,nfaces
812  if (cbc(f,e,1).ne.'E '.and.cbc(f,e,1).ne.' ')
813  $ call facev (vcmask,e,f,0.0,lxc,lyc,lzc)
814  enddo
815  enddo
816  call fgslib_gs_op(gshlc,vcmask,1,2,0)
817 
818  do e=1,nelv
819  do j=1,nedge
820  call rzero(xyz,9)
821  do i=1,3
822  xyz(1,i) = xmc(ke(i,j),1,1,e)
823  xyz(2,i) = ymc(ke(i,j),1,1,e)
824  if (ldim.eq.3) xyz(3,i) = zmc(ke(i,j),1,1,e)
825  enddo
826 
827  if (vcmask(ke(2,j),1,1,e).gt.0.5) then
828  call fixedgs(xyz,nxv,nyv,nzv)
829  xmc(ke(2,j),1,1,e) = nxv
830  ymc(ke(2,j),1,1,e) = nyv
831  if (ldim.eq.3) zmc(ke(2,j),1,1,e) = nzv
832  endif
833  enddo
834 
835 c fix face center and body center
836  zg(1) = -1.
837  zg(2) = 0.
838  zg(3) = 1.
839  call copy(xd,xmc(1,1,1,e),lxc*lyc*lzc)
840  call copy(yd,ymc(1,1,1,e),lxc*lyc*lzc)
841  if (ldim.eq.3) call copy(zd,zmc(1,1,1,e),lxc*lyc*lzc)
842  call gh_face_extend(xd,zg,3,2,w1,w2)
843  call gh_face_extend(yd,zg,3,2,w1,w2)
844  if (ldim.eq.3) call gh_face_extend(zd,zg,3,2,w1,w2)
845  kin(1) = 5
846  kin(2) = 11
847  kin(3) = 13
848  kin(4) = 15
849  kin(5) = 17
850  kin(6) = 23
851  kin(7) = 14
852  do i=1,(ldim-2)*6+1
853  xmc(kin(i),1,1,e) = xd(kin(i))
854  ymc(kin(i),1,1,e) = yd(kin(i))
855  if (ldim.eq.3) zmc(kin(i),1,1,e) = zd(kin(i))
856  enddo
857  enddo
858 
859  return
860  end
861 c-----------------------------------------------------------------------
862  subroutine fixedgs(xyz,nx,ny,nz)
863  real xyz(3,0:2) ! Coordinates
864  real x0(3),x1(3),x2(3),v2(3),v1(3),l02,l01,n1i,u1(3),u2(3)
865  real nx,ny,nz
866  real h,tol,h2
867 
868  call copy(x0,xyz(1,0),3)
869  call copy(x1,xyz(1,1),3)
870  call copy(x2,xyz(1,2),3)
871  ! ^ x2
872  call sub3(v1,x1,x0,3) ! v1=x1-x0 / \ v2 p1 = x0+v2*dot
873  call sub3(v2,x2,x0,3) ! v2=x2-x0 x1 .<-- x0
874  ! v1
875  l02 = vlsc2(v2,v2,3) ! || x2-x0 ||
876  if (l02.gt.0) then
877  l02 = sqrt(l02)
878  scl = 1./l02
879  call cmult2(u2,v2,scl,3) ! Unit tangent
880  endif
881 
882  l01 = vlsc2(v1,v1,3) ! || x1-x0 ||
883  if (l01.gt.0) then
884  l01 = sqrt(l01)
885  scl = 1./l01
886  call cmult2(u1,v1,scl,3) ! Unit tangent
887  endif
888 
889  dot = vlsc2(v1,u2,3)
890 
891  if (dot.le.0) then
892  write(6,*) 'SMOOTHER-ERROR 1 IN SHIFT - YOU SHOULD ABORT'
893  return
894  elseif (dot.gt.l02) then
895  write(6,*) 'SMOOTHER-ERROR 2 IN SHIFT - YOU SHOULD ABORT'
896  return
897  endif
898  h = 0.
899  h2 = 0.
900  tol = 1.e-8
901  do i=1,3
902  p1i = x0(i) + u2(i)*dot ! Projection of x1 onto [x0,x2]
903  n1i = x1(i) - p1i ! Normal vector
904  h = h + n1i**2
905  h2 = h2+ (x2(i)-x0(i))**2
906  xmi = 0.5*(x0(i)+x2(i))
907  x1(i) = xmi + n1i ! X1 point shifted to be centered at midpoint
908  enddo
909  if (h.le.h2*tol) then
910  x1(1) = 0.5*(x0(1)+x2(1))
911  x1(2) = 0.5*(x0(2)+x2(2))
912  x1(3) = 0.5*(x0(3)+x2(3))
913  endif
914 
915 
916  nx = x1(1)
917  ny = x1(2)
918  nz = x1(3)
919 
920  return
921  end
922 c-----------------------------------------------------------------------
923  subroutine gradf(f2,dfdx,x8,y8,z8,mlt,gshl,siz,opt,h)
924  include 'SIZE'
925  include 'TOTAL'
926  real x8(2**ldim,lelv*(2**ldim)),
927  $ y8(2**ldim,lelv*(2**ldim)),
928  $ z8(2**ldim,lelv*(2**ldim)),
929  $ mlt(2**ldim,lelv*(2**ldim)),
930  $ dfdx(2**ldim,lelv*(2**ldim),ldim)
931  real xt(2**ldim),
932  $ yt(2**ldim),
933  $ zt(2**ldim)
934  integer siz,opt,vertex,gshl,e,eg,e0,f
935  real par(siz)
936  real f1,fl,gl,f2,h(2**ldim,lelv*(2**ldim))
937 
938  f1 = 0
939  do e=1,nelv*(2**ldim)
940  if (opt.eq.1)
941  $ call get_jac(fl,x8(1,e),y8(1,e),z8(1,e),siz,e,0)
942  if (opt.eq.2) call get_len(fl,x8(1,e),y8(1,e),z8(1,e),siz)
943  f1 = f1+fl
944  call copy(xt,x8(1,e),2**ldim)
945  call copy(yt,y8(1,e),2**ldim)
946  call copy(zt,z8(1,e),2**ldim)
947  do j=1,2**ldim
948  xt(j) = x8(j,e)+h(j,e)
949  if (opt.eq.1) call get_jac(gl,xt,yt,zt,siz,e,1)
950  if (opt.eq.2) call get_len(gl,xt,yt,zt,siz)
951  xt(j) = x8(j,e)-h(j,e)
952  if (opt.eq.1) call get_jac(fl,xt,yt,zt,siz,e,1)
953  if (opt.eq.2) call get_len(fl,xt,yt,zt,siz)
954  xt(j) = x8(j,e)
955  dfdx(j,e,1) = (gl-fl)/(2.*h(j,e))
956 
957  yt(j) = y8(j,e)+h(j,e)
958  if (opt.eq.1) call get_jac(gl,xt,yt,zt,siz,e,2)
959  if (opt.eq.2) call get_len(gl,xt,yt,zt,siz)
960  yt(j) = y8(j,e)-h(j,e)
961  if (opt.eq.1) call get_jac(fl,xt,yt,zt,siz,e,2)
962  if (opt.eq.2) call get_len(fl,xt,yt,zt,siz)
963  yt(j) = y8(j,e)
964  dfdx(j,e,2) = (gl-fl)/(2.*h(j,e))
965 
966  if (ldim.eq.3) then
967  zt(j) = z8(j,e)+h(j,e)
968  if (opt.eq.1) call get_jac(gl,xt,yt,zt,siz,e,3)
969  if (opt.eq.2) call get_len(gl,xt,yt,zt,siz)
970  zt(j) = z8(j,e)-h(j,e)
971  if (opt.eq.1) call get_jac(fl,xt,yt,zt,siz,e,3)
972  if (opt.eq.2) call get_len(fl,xt,yt,zt,siz)
973  zt(j) = z8(j,e)
974  dfdx(j,e,ldim) = (gl-fl)/(2.*h(j,e))
975  endif
976  enddo
977  enddo
978 
979  call fgslib_gs_op(gshl,dfdx(1,1,1),1,1,0)
980  call fgslib_gs_op(gshl,dfdx(1,1,2),1,1,0)
981  if (ldim.eq.3) call fgslib_gs_op(gshl,dfdx(1,1,ldim),1,1,0)
982 
983  f2 = glsum(f1,1)
984 
985  return
986  end
987 c-----------------------------------------------------------------------
988  subroutine get_jac(val,x,y,z,siz,el,dire)
989  include 'SIZE'
990  include 'TOTAL'
991 c This routine compiles the jacobian matrix and then does the
992 c frobenius norm of the matrix times that of the inverse
993  integer i,j,k,siz,el,dire
994 
995  real par(siz),det(siz)
996  real jm(ldim,ldim),jin(ldim,ldim),frn(ldim**2)
997  real x(2**ldim),y(2**ldim),z(2**ldim),jac(2**ldim)
998  real fr1,fr2,sum1,dumc,val
999 
1000  integer bzindx(24),czindx(24)
1001  integer penalty
1002  SAVE bzindx
1003  DATA bzindx / 2,3,5, 1,4,6, 4,1,7, 3,2,8,
1004  $ 6,7,1, 5,8,2, 8,5,3, 7,6,4 /
1005 c bzindx tells which node is node connected to in r,s,t direction
1006 c example: node 1 is connected to 2,3,5; 2 to 1,4,6 and so on
1007  SAVE czindx
1008  DATA czindx / 1,1,1, -1,1,1, 1,-1,1, -1,-1,1,
1009  $ 1,1,-1, -1,1,-1, 1,-1,-1, -1,-1,-1 /
1010 c this tells which node is to be subtracted.
1011 
1012  if (siz.ne.2**ldim) then
1013  write(6,*) 'siz value wrong input into get_jac'
1014  write(6,*) siz,'sizval',2**ldim,'should be this'
1015  call exitt
1016  endif
1017 
1018  do i=1,2**ldim
1019  ind1 = (i-1)*3 !tells where to look in b/czindx
1020  do j=1,ldim
1021  ind2 = ind1+j
1022  jm(1,j) = 0.5*czindx(ind2)*(x(bzindx(ind2))-x(i))
1023  jm(2,j) = 0.5*czindx(ind2)*(y(bzindx(ind2))-y(i))
1024  if (ldim.eq.3)
1025  $ jm(ldim,j) = 0.5*czindx(ind2)*(z(bzindx(ind2))-z(i))
1026  enddo
1027 c jacobian matrix has been calculated at this point.
1028 c now calculate determinant
1029  if (ldim.eq.3) then
1030  jac(i)=jm(1,1)*(jm(2,2)*jm(ldim,ldim)-jm(2,ldim)*jm(ldim,2))+
1031  $ jm(1,ldim)*(jm(2,1)*jm(ldim,2)-jm(2,2)*jm(ldim,1)) +
1032  $ jm(1,2)*(jm(2,ldim)*jm(ldim,1)-jm(2,1)*jm(ldim,ldim))
1033  else
1034  jac(i)=jm(1,1)*jm(2,2)-jm(2,1)*jm(1,2)
1035  endif
1036 c calculate inverse matrix
1037  if (ldim.eq.3) then
1038  jin(1,1) = jm(2,2)*jm(ldim,ldim)-jm(ldim,2)*jm(2,ldim)
1039  jin(2,1) = jm(2,ldim)*jm(ldim,1)-jm(ldim,ldim)*jm(2,1)
1040  jin(ldim,1) = jm(2,1)*jm(ldim,2)-jm(ldim,1)*jm(2,2)
1041 
1042  jin(1,2) = jm(1,ldim)*jm(ldim,2)-jm(ldim,ldim)*jm(1,2)
1043  jin(2,2) = jm(1,1)*jm(ldim,ldim)-jm(ldim,1)*jm(1,ldim)
1044  jin(ldim,2) = jm(1,2)*jm(ldim,1)-jm(ldim,2)*jm(1,1)
1045 
1046  jin(1,ldim) = jm(1,2)*jm(2,ldim)-jm(2,2)*jm(1,ldim)
1047  jin(2,ldim) = jm(1,ldim)*jm(2,1)-jm(2,ldim)*jm(1,1)
1048  jin(ldim,ldim) = jm(1,1)*jm(2,2)-jm(2,1)*jm(1,2)
1049  else
1050  jin(1,1) = jm(2,2)
1051  jin(2,1) = -jm(2,1)
1052 
1053  jin(1,2) = -jm(1,2)
1054  jin(2,2) = jm(1,1)
1055  endif
1056 
1057  dumc = 1/jac(i)
1058  call cmult(jin,dumc,ldim**2) !scale inverse by inv det
1059 
1060  call rzero(frn,ldim**2)
1061  call col3(frn,jm,jm,ldim**2) !square the entries
1062  sum1 = vlsum(frn,ldim**2)
1063  fr1 = sqrt(sum1) !squareroot
1064 
1065  call rzero(frn,ldim**2)
1066  call col3(frn,jin,jin,ldim**2)
1067  sum1 = vlsum(frn,ldim**2)
1068  fr2 = sqrt(sum1)
1069 
1070  par(i) = fr1*fr2
1071  par(i) = (par(i)/ldim)**2
1072 
1073  enddo
1074 
1075  val = vlsum(par,siz)
1076  val = val/siz !normalize to 1 for each element
1077  val=abs(val)
1078 
1079  return
1080  end
1081 c-----------------------------------------------------------------------
1082  subroutine get_len(val,x,y,z,siz)
1083  include 'SIZE'
1084  include 'TOTAL'
1085  integer siz
1086  real x(2**ldim),y(2**ldim),z(2**ldim)
1087  real par(siz),xm,ym,zm,val
1088 
1089  integer azindx(24)
1090  SAVE azindx
1091  DATA azindx / 1,2 ,2,4, 3,4, 1,3, 2,6, 6,8,
1092  $ 4,8, 5,6, 7,8, 5,7, 1,5, 3,7 /
1093 c azindx tells what nodes make up each edge
1094 
1095  do i=1,siz
1096  ind = (i-1)*2
1097  xm = x(azindx(ind+1))-x(azindx(ind+2))
1098  ym = y(azindx(ind+1))-y(azindx(ind+2))
1099  zm = z(azindx(ind+1))-z(azindx(ind+2))
1100  par(i) = sqrt(xm*xm+ym*ym+zm*zm)
1101  par(i) = par(i)**2
1102  enddo
1103 
1104  val = vlsum(par,siz)
1105 
1106  return
1107  end
1108 c-----------------------------------------------------------------------
1109  subroutine get_nodscale(scalek,x8,y8,z8,gshl)
1110  include 'SIZE'
1111  include 'TOTAL'
1112  real x8(2**ldim,lelv*(2**ldim)),
1113  $ y8(2**ldim,lelv*(2**ldim)),
1114  $ z8(2**ldim,lelv*(2**ldim))
1115  real curval,xm,ym,zm,work1(1),dum,wrk1
1116 c real scalek
1117  real scalek(2**ldim,lelv*(2**ldim))
1118  integer bzindx(24),e,gshl
1119  DATA bzindx / 2,3,5, 1,4,6, 4,1,7, 3,2,8,
1120  $ 6,7,1, 5,8,2, 8,5,3, 7,6,4 /
1121 c bzindx tells what node is connected to what node
1122 
1123  n1 = nelv*(2**ldim)
1124  n2 = n1*(2**ldim)
1125  curval = 1e+10
1126  call rone(scalek,n2)
1127  call cmult(scalek,curval,n2)
1128 
1129  do e=1,nelv*(2**ldim)
1130  do i=1,2**ldim
1131  curval = 1e+10
1132  ind = (i-1)*3
1133  do j=1,ldim
1134  xm = x8(i,e)-x8(bzindx(ind+j),e)
1135  ym = y8(i,e)-y8(bzindx(ind+j),e)
1136  zm = 0.
1137  if (ldim.eq.3) zm = z8(i,e)-z8(bzindx(ind+j),e)
1138  dum = sqrt(xm*xm+ym*ym+zm*zm)
1139  curval = min(dum,curval)
1140  enddo
1141  scalek(i,e) = curval
1142  enddo
1143  enddo
1144 c
1145  fac = 1.e-2
1146  call cmult(scalek,fac,n2)
1147  call fgslib_gs_op(gshl,scalek,1,3,0)
1148 
1149  return
1150  end
1151 c-----------------------------------------------------------------------
1152  subroutine genmask(nodmask,mlt,gshl,mltc,gshlc)
1153  include 'SIZE'
1154  include 'TOTAL'
1155  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1156  integer gs_handle
1157  integer vertex(1)
1158  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1159  common /ivrtx/ vertex
1160  integer*8 glo_num(lxc*lyc*lzc,lelv),ngv
1161  integer*8 glo_numk2(2**ldim,lelv*(2**ldim))
1162 
1163  real mlt(2**ldim,lelv*(2**ldim)),
1164  $ nodmask(2**ldim,lelv*(2**ldim))
1165  real mltc(lxc*lyc*lzc*lelv),
1166  $ wrkmask(lxc*lyc*lzc*lelv)
1167  integer gshlc,gshl,e,eg,e0,f
1168 
1169  integer zindx(64)
1170  SAVE zindx
1171  DATA zindx / 1, 2, 4, 5, 10, 11, 13, 14,
1172  $ 2, 3, 5, 6, 11, 12, 14, 15,
1173  $ 4, 5, 7, 8, 13, 14, 16, 17,
1174  $ 5, 6, 8, 9, 14, 15, 17, 18,
1175  $ 10, 11, 13, 14, 19, 20, 22, 23,
1176  $ 11, 12, 14, 15, 20, 21, 23, 24,
1177  $ 13, 14, 16, 17, 22, 23, 25, 26,
1178  $ 14, 15, 17, 18, 23, 24, 26, 27 /
1179 
1180 
1181  call setupds_center(gs_handle,lxc,lyc,lzc,nelv,
1182  $ nelgv,vertex,glo_num)
1183 
1184 c setup the mask first so that it can be distribute as well
1185  zero = 0.
1186  call rone(wrkmask,lxc*lyc*lzc*nelv)
1187  do e=1,nelv
1188  do f=1,2*ldim
1189  if(cbc(f,e,1).ne.'E '.and.cbc(f,e,1).ne.' ')then
1190  call facev(wrkmask,e,f,zero,lxc,lyc,lzc)
1191  endif
1192  enddo
1193  enddo
1194 
1195  n = (lxc*lyc*lzc)*nelv
1196  call rone(mltc,n)
1197  call fgslib_gs_setup(gshlc,glo_num,n,nekcomm,np)
1198  call fgslib_gs_op(gshlc,mltc,1,1,0)
1199 
1200  call fgslib_gs_op(gshlc,wrkmask,1,2,0)
1201  call xmtox8(wrkmask,nodmask)
1202 
1203  n = 0
1204  do e=1,nelv !each element
1205  do j=1,2**ldim !broken down into 4 quads/8hexes
1206  n = n+1
1207  ind1 = (j-1)*8
1208  do k=1,2**ldim
1209  glo_numk2(k,n) = glo_num(zindx(ind1+k),e)
1210  enddo
1211  enddo
1212  enddo
1213 
1214  n = (2**ldim)*(nelv*(2**ldim))
1215  call fgslib_gs_setup(gshl,glo_numk2,n,nekcomm,np)
1216  call rone(mlt,n)
1217  call fgslib_gs_op(gshl,mlt,1,1,0) ! '+'
1218  xmlt = glmax(mlt,n)
1219  call invcol1(mlt,n)
1220 
1221  call fgslib_gs_op(gshl,nodmask,1,2,0)
1222 
1223  return
1224  end
1225 c-----------------------------------------------------------------------
1226  subroutine setupds_center(gs_handle,nx,ny,nz,nel,melg,
1227  $ vertex,glo_num)
1228  include 'SIZE'
1229  include 'INPUT'
1230  include 'PARALLEL'
1231  include 'NONCON'
1232  integer gs_handle
1233  integer vertex(1)
1234  integer*8 glo_num(1),ngv
1235  common /nekmpi/ mid,mp,nekcomm,nekgroup,nekreal
1236 
1237  n = nx*ny*nz*nel
1238  call set_vert(glo_num,ngv,nx,nel,vertex,.true.)
1239  call fgslib_gs_setup(gs_handle,glo_num,n,nekcomm,mp)
1240  return
1241  end
1242 c-----------------------------------------------------------------------
1243  subroutine fastlap(iter,nodmask,mlt,gshl,dis,lapinv,mltc,gshlc)
1244  include 'SIZE'
1245  include 'TOTAL'
1246  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1247  real xmc(lxc,lyc,lzc,lelv),
1248  $ ymc(lxc,lyc,lzc,lelv),
1249  $ zmc(lxc,lyc,lzc,lelv)
1250  common /coarsemesh/ xmc,ymc,zmc
1251 
1252  real x8(2,2,ldim-1,lelv*(2**ldim)),
1253  $ y8(2,2,ldim-1,lelv*(2**ldim)),
1254  $ z8(2,2,ldim-1,lelv*(2**ldim))
1255  common /coarsemesh8/ x8,y8,z8
1256 
1257  real dxx(2**ldim,lelv*(2**ldim)),
1258  $ dyy(2**ldim,lelv*(2**ldim)),
1259  $ dzz(2**ldim,lelv*(2**ldim)),
1260  $ dis(2**ldim,lelv*(2**ldim)),
1261  $ mlt(2**ldim,lelv*(2**ldim)),
1262  $ nodmask(2**ldim,lelv*(2**ldim)),
1263  $ mltc(lxc*lyc*lzc,lelv)
1264  integer z,e,f,gshl,lapinv,kerr,gshlc
1265  real xbar,ybar,zbar,sfac
1266 
1267  if (nid.eq.0.and.loglevel.ge.2)
1268  $ write(6,*) 'laplacian smoothing started'
1269 
1270  n = 2**ldim
1271  sfac = 0.01
1272  do k=1,iter
1273  do e=1,nelv*(2**ldim)
1274  xbar = vlsum(x8(1,1,1,e),n)/(n)
1275  ybar = vlsum(y8(1,1,1,e),n)/(n)
1276  if (ldim.eq.3) zbar = vlsum(z8(1,1,1,e),n)/(n)
1277  do i=1,n
1278  dxx(i,e) = sfac*(xbar-x8(i,1,1,e))
1279  dyy(i,e) = sfac*(ybar-y8(i,1,1,e))
1280  if (ldim.eq.3) dzz(i,e) = sfac*(zbar-z8(i,1,1,e))
1281  enddo
1282  enddo
1283 
1284  n2 = nelv*(2**ldim)*(n)
1285 
1286  call col2(dxx,nodmask,n2)
1287  call col2(dxx,dis,n2)
1288  call dsavg_general(dxx,mlt,gshl)
1289  call add2(x8,dxx,n2)
1290 
1291  call col2(dyy,nodmask,n2)
1292  call col2(dyy,dis,n2)
1293  call dsavg_general(dyy,mlt,gshl)
1294  call add2(y8,dyy,n2)
1295 
1296  if (ldim.eq.3) then
1297  call col2(dzz,nodmask,n2)
1298  call col2(dzz,dis,n2)
1299  call dsavg_general(dzz,mlt,gshl)
1300  call add2(z8,dzz,n2)
1301  endif
1302 
1303  dxm =glamax(dxx,n2)
1304  dym =glamax(dyy,n2)
1305  if (ldim.eq.3) then
1306  dzm =glamax(dzz,n2)
1307  else
1308  dzm = 0.
1309  endif
1310  if (nio.eq.0.and.loglevel.ge.3) write(6,1) k,dxm,dym,dzm
1311  1 format(i5,1p3e12.3,' dxm')
1312 
1313  enddo
1314 
1315  call x8toxm(xmc,x8)
1316  call x8toxm(ymc,y8)
1317  if (ldim.eq.3) call x8toxm(zmc,z8)
1318  call fixcurs(mltc,gshlc)
1319  call xmtox8(xmc,x8)
1320  call xmtox8(ymc,y8)
1321  if (ldim.eq.3) call xmtox8(zmc,z8)
1322 
1323  call glmapm1chkinv(kerr)
1324  if (kerr.gt.0.) then
1325  lapinv = 1 !Terminate Laplacian smoothing
1326  if (nid.eq.0.and.loglevel.ge.2) write(6,*) 'terminating Laplacian'
1327  else
1328  call genbackupmesh
1329  endif
1330 
1331  return
1332  end
1333 c-----------------------------------------------------------------------
1334  subroutine dsavg_general(u,mlt,gshl)
1335  include 'SIZE'
1336  include 'TOTAL'
1337  integer gshl
1338  real u(2**ldim,lelv*(2**ldim))
1339  real mlt(2**ldim,lelv*(2**ldim))
1340 
1341  call fgslib_gs_op(gshl,u,1,1,0) ! '+'
1342  call col2(u,mlt,(2**ldim)*nelv*(2**ldim))
1343 
1344  return
1345  end
1346 c-----------------------------------------------------------------------
1347  subroutine xmtox8(xd,x8)
1348  include 'SIZE'
1349  include 'TOTAL'
1350  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1351  real x8(2**ldim,lelv*(2**ldim))
1352  real xd(lxc*lyc*lzc,lelv)
1353  integer e,k,n,j,ind1
1354  integer zindx(64)
1355  SAVE zindx
1356  DATA zindx / 1, 2, 4, 5, 10, 11, 13, 14,
1357  $ 2, 3, 5, 6, 11, 12, 14, 15,
1358  $ 4, 5, 7, 8, 13, 14, 16, 17,
1359  $ 5, 6, 8, 9, 14, 15, 17, 18,
1360  $ 10, 11, 13, 14, 19, 20, 22, 23,
1361  $ 11, 12, 14, 15, 20, 21, 23, 24,
1362  $ 13, 14, 16, 17, 22, 23, 25, 26,
1363  $ 14, 15, 17, 18, 23, 24, 26, 27 /
1364 
1365  n = 0
1366  do e=1,nelv !each element
1367  do j=1,2**ldim !broken down into 4 quads/8hexes
1368  n = n+1
1369  ind1 = (j-1)*8
1370  do k=1,2**ldim
1371  x8(k,n) = xd(zindx(ind1+k),e)
1372  enddo
1373  enddo
1374  enddo
1375 
1376  return
1377  end
1378 c-----------------------------------------------------------------------
1379  subroutine x8toxm(xd,x8)
1380  include 'SIZE'
1381  include 'TOTAL'
1382  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1383  real x8(2**ldim,lelv*(2**ldim))
1384  real xd(lxc*lyc*lzc,lelv)
1385  integer zindx((2**3)**2)
1386  integer e,k,n,j,ind1
1387  DATA zindx / 1, 2, 4, 5, 10, 11, 13, 14,
1388  $ 2, 3, 5, 6, 11, 12, 14, 15,
1389  $ 4, 5, 7, 8, 13, 14, 16, 17,
1390  $ 5, 6, 8, 9, 14, 15, 17, 18,
1391  $ 10, 11, 13, 14, 19, 20, 22, 23,
1392  $ 11, 12, 14, 15, 20, 21, 23, 24,
1393  $ 13, 14, 16, 17, 22, 23, 25, 26,
1394  $ 14, 15, 17, 18, 23, 24, 26, 27 /
1395 
1396 
1397  n = 0
1398  do e=1,nelv !each element
1399  do j=1,2**ldim !broken down into 4 quads/8hexes
1400  n = n+1
1401  ind1 = (j-1)*8
1402  do k=1,2**ldim
1403  xd(zindx(ind1+k),e) = x8(k,n)
1404  enddo
1405  enddo
1406  enddo
1407 
1408  return
1409  end
1410 c-----------------------------------------------------------------------
1411  subroutine col3c(a,b,c,d,n)
1412  real a(1),b(1),c(1),d
1413 
1414  do i=1,n
1415  a(i)=b(i)*c(i)*d
1416  enddo
1417 
1418  return
1419  end
1420 c-----------------------------------------------------------------------
1421  subroutine glmapm1chkinv(kerr)
1422  include 'SIZE'
1423  include 'GEOM'
1424  include 'INPUT'
1425  include 'SOLN'
1426 C
1427 C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
1428 C share the same array structure in Scratch Common /SCRNS/.
1429  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1430  parameter(nxc=3,nyc=3,nzc=1+(ldim-2)*(nxc-1))
1431  real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
1432  $ zmc(lxc,lyc,lzc,lelv)
1433  common /coarsemesh/ xmc,ymc,zmc
1434 C
1435  real XRMc(LXc,LYc,LZc,LELv)
1436  $ , YRMc(LXc,LYc,LZc,LELv)
1437  $ , xsmc(lxc,lyc,lzc,lelv)
1438  $ , ysmc(lxc,lyc,lzc,lelv)
1439  $ , xtmc(lxc,lyc,lzc,lelv)
1440  $ , ytmc(lxc,lyc,lzc,lelv)
1441  $ , zrmc(lxc,lyc,lzc,lelv)
1442  real ZSMc(LXc,LYc,LZc,LELv)
1443  $ , ZTMc(LXc,LYc,LZc,LELv)
1444  $ , jacmc(lxc,lyc,lzc,lelv)
1445  real rxmc(lxc,lyc,lzc,lelv)
1446  $ , rymc(lxc,lyc,lzc,lelv)
1447  $ , rzmc(lxc,lyc,lzc,lelv)
1448  $ , sxmc(lxc,lyc,lzc,lelv)
1449  $ , symc(lxc,lyc,lzc,lelv)
1450  $ , szmc(lxc,lyc,lzc,lelv)
1451  $ , txmc(lxc,lyc,lzc,lelv)
1452  $ , tymc(lxc,lyc,lzc,lelv)
1453  $ , tzmc(lxc,lyc,lzc,lelv)
1454  integer kerr,ierr
1455 C
1456  nxyc = nxc*nyc
1457  nyzc = nyc*nzc
1458  nxyzc = nxc*nyc*nzc
1459  ntotc = nxyzc*nelv
1460 C
1461  CALL xyzrstc (xrmc,yrmc,zrmc,xsmc,ysmc,zsmc,xtmc,ytmc,ztmc,
1462  $ ifaxis)
1463 
1464 
1465 C
1466  IF (ndim.EQ.2) THEN
1467  CALL rzero (jacmc,ntotc)
1468  CALL addcol3 (jacmc,xrmc,ysmc,ntotc)
1469  CALL subcol3 (jacmc,xsmc,yrmc,ntotc)
1470  CALL copy (rxmc,ysmc,ntotc)
1471  CALL copy (rymc,xsmc,ntotc)
1472  CALL chsign (rymc,ntotc)
1473  CALL copy (sxmc,yrmc,ntotc)
1474  CALL chsign (sxmc,ntotc)
1475  CALL copy (symc,xrmc,ntotc)
1476  CALL rzero (rzmc,ntotc)
1477  CALL rzero (szmc,ntotc)
1478  CALL rone (tzmc,ntotc)
1479  ELSE
1480  CALL rzero (jacmc,ntotc)
1481  CALL addcol4 (jacmc,xrmc,ysmc,ztmc,ntotc)
1482  CALL addcol4 (jacmc,xtmc,yrmc,zsmc,ntotc)
1483  CALL addcol4 (jacmc,xsmc,ytmc,zrmc,ntotc)
1484  CALL subcol4 (jacmc,xrmc,ytmc,zsmc,ntotc)
1485  CALL subcol4 (jacmc,xsmc,yrmc,ztmc,ntotc)
1486  CALL subcol4 (jacmc,xtmc,ysmc,zrmc,ntotc)
1487  CALL ascol5 (rxmc,ysmc,ztmc,ytmc,zsmc,ntotc)
1488  CALL ascol5 (rymc,xtmc,zsmc,xsmc,ztmc,ntotc)
1489  CALL ascol5 (rzmc,xsmc,ytmc,xtmc,ysmc,ntotc)
1490  CALL ascol5 (sxmc,ytmc,zrmc,yrmc,ztmc,ntotc)
1491  CALL ascol5 (symc,xrmc,ztmc,xtmc,zrmc,ntotc)
1492  CALL ascol5 (szmc,xtmc,yrmc,xrmc,ytmc,ntotc)
1493  CALL ascol5 (txmc,yrmc,zsmc,ysmc,zrmc,ntotc)
1494  CALL ascol5 (tymc,xsmc,zrmc,xrmc,zsmc,ntotc)
1495  CALL ascol5 (tzmc,xrmc,ysmc,xsmc,yrmc,ntotc)
1496  ENDIF
1497 C
1498  kerr = 0
1499  DO 500 ie=1,nelv
1500  CALL chkjacinv(jacmc(1,1,1,ie),nxyzc,ie,xmc(1,1,1,ie),
1501  $ ymc(1,1,1,ie),zmc(1,1,1,ie),ndim,ierr)
1502  if (ierr.ne.0) kerr = kerr+1
1503  500 CONTINUE
1504  kerr = iglsum(kerr,1)
1505 
1506  RETURN
1507  END
1508 C-----------------------------------------------------------------------
1509  subroutine chkjacinv(jac,n,iel,X,Y,Z,ND,IERR)
1510 c
1511  include 'SIZE'
1512  include 'PARALLEL'
1513 C
1514 C Check the array JAC for a change in sign.
1515 C
1516  REAL JAC(N),x(1),y(1),z(1)
1517  integer ierr
1518 c
1519  ierr = 0
1520  sign = jac(1)
1521  DO 100 i=2,n
1522  IF (sign*jac(i).LE.0.0) THEN
1523  ierr = 1
1524  ENDIF
1525  100 CONTINUE
1526 
1527  RETURN
1528  END
1529 C-----------------------------------------------------------------------
1530  subroutine xyzrstc(xrmc,yrmc,zrmc,xsmc,ysmc,zsmc,
1531  $ XTMc,YTMc,ZTMc,IFAXIS)
1532 C-----------------------------------------------------------------------
1533 C
1534 C Compute global-to-local derivatives on mesh 1.
1535 C
1536 C-----------------------------------------------------------------------
1537  include 'SIZE'
1538 C
1539  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1540  dimension xrmc(lxc,lyc,lzc,1),yrmc(lxc,lyc,lzc,1)
1541  $ , zrmc(lxc,lyc,lzc,1),xsmc(lxc,lyc,lzc,1)
1542  $ , ysmc(lxc,lyc,lzc,1),zsmc(lxc,lyc,lzc,1)
1543  $ , xtmc(lxc,lyc,lzc,1),ytmc(lxc,lyc,lzc,1)
1544  $ , ztmc(lxc,lyc,lzc,1)
1545  LOGICAL IFAXIS
1546  real dxmc(3,3),dymc(3,3),dzmc(3,3)
1547  real dxtmc(3,3),dytmc(3,3),dztmc(3,3)
1548  common /coarseders/ dxmc,dymc,dzmc,dxtmc,dytmc,dztmc
1549  real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
1550  $ zmc(lxc,lyc,lzc,lelv)
1551  common /coarsemesh/ xmc,ymc,zmc
1552 C
1553  nxyc=lxc*lyc
1554  nyzc=lyc*lzc
1555 C
1556  DO 100 iel=1,nelv
1557 C
1558  CALL mxm (dxmc,lxc,xmc(1,1,1,iel),lxc,xrmc(1,1,1,iel),nyzc)
1559  CALL mxm (dxmc,lxc,ymc(1,1,1,iel),lxc,yrmc(1,1,1,iel),nyzc)
1560  CALL mxm (dxmc,lxc,zmc(1,1,1,iel),lxc,zrmc(1,1,1,iel),nyzc)
1561 C
1562  DO 10 iz=1,lzc
1563  CALL mxm (xmc(1,1,iz,iel),lxc,dytmc,lyc,xsmc(1,1,iz,iel),lyc)
1564  CALL mxm (ymc(1,1,iz,iel),lxc,dytmc,lyc,ysmc(1,1,iz,iel),lyc)
1565  CALL mxm (zmc(1,1,iz,iel),lxc,dytmc,lyc,zsmc(1,1,iz,iel),lyc)
1566  10 CONTINUE
1567 C
1568  IF (ldim.EQ.3) THEN
1569  CALL mxm (xmc(1,1,1,iel),nxyc,dztmc,lzc,xtmc(1,1,1,iel),lzc)
1570  CALL mxm (ymc(1,1,1,iel),nxyc,dztmc,lzc,ytmc(1,1,1,iel),lzc)
1571  CALL mxm (zmc(1,1,1,iel),nxyc,dztmc,lzc,ztmc(1,1,1,iel),lzc)
1572  ELSE
1573  CALL rzero (xtmc(1,1,1,iel),nxyc)
1574  CALL rzero (ytmc(1,1,1,iel),nxyc)
1575  CALL rone (ztmc(1,1,1,iel),nxyc)
1576  ENDIF
1577 C
1578  100 CONTINUE
1579 C
1580  RETURN
1581  END
1582 C-----------------------------------------------------------------------
1583  subroutine xm1toxmc
1584  include 'SIZE'
1585  include 'TOTAL'
1586  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1587  real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
1588  $ zmc(lxc,lyc,lzc,lelv)
1589  common /coarsemesh/ xmc,ymc,zmc
1590  integer e
1591 
1592  if (ldim.eq.2) then
1593  do e=1,nelv
1594  call int_fine_to_coarse_2d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
1595  call int_fine_to_coarse_2d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
1596  enddo
1597  else
1598  do e=1,nelv
1599  call int_fine_to_coarse_3d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
1600  call int_fine_to_coarse_3d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
1601  call int_fine_to_coarse_3d(zm1(1,1,1,e),zmc(1,1,1,e),lx1,3)
1602  enddo
1603  endif
1604 
1605  RETURN
1606  END
1607 C-----------------------------------------------------------------------
1608  subroutine xmctoxm1
1609  include 'SIZE'
1610  include 'TOTAL'
1611  parameter(lxc=3,lyc=3,lzc=1+(ldim-2)*(lxc-1))
1612  real xmc(lxc,lyc,lzc,lelv),ymc(lxc,lyc,lzc,lelv),
1613  $ zmc(lxc,lyc,lzc,lelv)
1614  common /coarsemesh/ xmc,ymc,zmc
1615  integer e
1616 
1617  if (ldim.eq.2) then
1618  do e=1,nelv
1619  call int_coarse_to_fine_2d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
1620  call int_coarse_to_fine_2d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
1621  enddo
1622  else
1623  do e=1,nelv
1624  call int_coarse_to_fine_3d(xm1(1,1,1,e),xmc(1,1,1,e),lx1,3)
1625  call int_coarse_to_fine_3d(ym1(1,1,1,e),ymc(1,1,1,e),lx1,3)
1626  call int_coarse_to_fine_3d(zm1(1,1,1,e),zmc(1,1,1,e),lx1,3)
1627  enddo
1628  endif
1629 
1630  RETURN
1631  END
1632 C-----------------------------------------------------------------------
1633  subroutine my_mv_mesh(dx,dy,dz)
1634  include 'SIZE'
1635  include 'TOTAL'
1636 
1637  parameter(lt = lx1*ly1*lz1*lelv)
1638  common /mrthoi/ napprx(2),nappry(2),napprz(2)
1639  common /mrthov/ apprx(lt,0:mxprev)
1640  $ , appry(lt,0:mxprev)
1641  $ , apprz(lt,0:mxprev)
1642  common /mstuff/ d(lt),h1(lt),h2(lt),mask(lt)
1643  real mask
1644  real mask2(lx1,ly1,lz1,lelv)
1645  real dx(1),dy(1),dz(1)
1646 
1647  integer e,f
1648 
1649  integer icalld
1650  save icalld
1651  data icalld /0/
1652 
1653  n = nx1*ny1*nz1*nelv
1654  nface = 2*ndim
1655 
1656  napprx(1)=0
1657  nappry(1)=0
1658  napprz(1)=0
1659 
1660  nxz = nx1*nz1
1661  nxyz = nx1*ny1*nz1
1662 
1663  volbl = 0. ! Volume of elements in boundary layer
1664  do e=1,nelv
1665  do f=1,nface
1666  if (cbc(f,e,1).eq.'W ') then
1667  srfbl = srfbl + vlsum(area(1,1,f,e),nxz )
1668  volbl = volbl + vlsum(bm1(1,1,1,e),nxyz)
1669  endif
1670  enddo
1671  enddo
1672  srfbl = glsum(srfbl,1) ! Sum over all processors
1673  volbl = glsum(volbl,1)
1674 
1675  delta = volbl / srfbl ! Avg thickness of b.l. elements
1676 
1677  call rone (h1,n)
1678  call rzero(h2,n)
1679  call sm_cheap_dist(d,1,'W ')
1680 
1681  deltap = 5*delta ! Protected b.l. thickness
1682  do i=1,n
1683  arg = -(d(i)/deltap)**2
1684  h1(i) = h1(i) + 9*exp(arg)
1685  enddo
1686 
1687  call rone(mask,n)
1688  do e=1,nelv
1689  do f=1,nface
1690  if (cbc(f,e,1).ne.'E '.and.cbc(f,e,1).ne.' ')
1691  $ call facev(mask,e,f,0.,nx1,ny1,nz1)
1692  enddo
1693  enddo
1694  call dsop(mask,'mul',nx1,ny1,nz1)
1695  call copy(mask2,mask,n)
1696  call chsign(mask2,n)
1697  call cadd(mask2,1.,n)
1698  call col2(dx,mask2,n)
1699  call col2(dy,mask2,n)
1700  call col2(dz,mask2,n)
1701  tol = -1.e-6
1702  call laplaceh('mshx',dx,h1,h2,mask,vmult,1,tol,500,apprx,napprx)
1703  call laplaceh('mshy',dy,h1,h2,mask,vmult,1,tol,500,appry,nappry)
1704  if (if3d)
1705  $ call laplaceh('mshz',dz,h1,h2,mask,vmult,1,tol,500,apprz,napprz)
1706 
1707  return
1708  end
1709 c-----------------------------------------------------------------------
1710  subroutine laplaceh
1711  $ (name,u,h1,h2,mask,mult,ifld,tli,maxi,approx,napprox)
1712  include 'SIZE'
1713  include 'TOTAL'
1714  include 'CTIMER'
1715 c
1716  character*4 name
1717  real u(1),h1(1),h2(1),mask(1),mult(1),approx (1)
1718  integer napprox(1)
1719 
1720  parameter(lt=lx1*ly1*lz1*lelv)
1721  common /scruz/ r(lt),ub(lt)
1722 
1723  logical ifstdh
1724  character*4 cname
1725  character*6 name6
1726 
1727  logical ifwt,ifvec
1728 
1729  call chcopy(cname,name,4)
1730  call capit (cname,4)
1731 
1732  call blank (name6,6)
1733  call chcopy(name6,name,4)
1734  ifwt = .true.
1735  ifvec = .false.
1736  isd = 1
1737  imsh = 1
1738  nel = nelfld(ifld)
1739 
1740  n = nx1*ny1*nz1*nel
1741 
1742  call copy (ub,u,n) ! ub = u on boundary
1743  call dsavg(ub) ! Make certain ub is in H1
1744  ! _
1745  call axhelm (r,ub,h1,h2,1,1) ! r = A*ub
1746 
1747  do i=1,n ! _
1748  r(i)=-r(i)*mask(i) ! r = -M*A*ub
1749  enddo
1750 
1751  call dssum (r,nx1,ny1,nz1) ! dssum rhs
1752 
1753  call project1
1754  $ (r,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
1755 
1756  tol = -abs(tli)
1757  if (nel.eq.nelv) then
1758  call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,binvm1)
1759  else
1760  call hmhzpf (name,u,r,h1,h2,mask,mult,imsh,tol,maxi,isd,bintm1)
1761  endif
1762 
1763  call project2
1764  $ (u,n,approx,napprox,h1,h2,mask,mult,ifwt,ifvec,name6)
1765 
1766  call add2(u,ub,n)
1767 
1768  return
1769  end
1770 c-----------------------------------------------------------------------
1771  subroutine sm_cheap_dist(d,ifld,b)
1772 c this is a copy of the routine in navier5.f, modified to be less
1773 c verbose
1774 
1775 
1776  include 'SIZE'
1777  include 'GEOM' ! Coordinates
1778  include 'INPUT' ! cbc()
1779  include 'TSTEP' ! nelfld
1780  include 'PARALLEL' ! gather-scatter handle for field "ifld"
1781 
1782  real d(lx1,ly1,lz1,lelt)
1783 
1784  character*3 b ! Boundary condition of interest
1785 
1786  integer e,eg,f
1787 
1788  nel = nelfld(ifld)
1789  n = lx1*ly1*lz1*nel
1790 
1791  call domain_size(xmin,xmax,ymin,ymax,zmin,zmax)
1792 
1793  xmn = min(xmin,ymin)
1794  xmx = max(xmax,ymax)
1795  if (if3d) xmn = min(xmn ,zmin)
1796  if (if3d) xmx = max(xmx ,zmax)
1797 
1798  big = 10*(xmx-xmn)
1799  call cfill(d,big,n)
1800 
1801  nface = 2*ldim
1802  do e=1,nel ! Set d=0 on walls
1803  do f=1,nface
1804  if (cbc(f,e,ifld).eq.b) call facev(d,e,f,0.,lx1,ly1,lz1)
1805  enddo
1806  enddo
1807 
1808  do ipass=1,10000
1809  dmax = 0
1810  nchange = 0
1811  do e=1,nel
1812  do k=1,lz1
1813  do j=1,ly1
1814  do i=1,lx1
1815  i0=max( 1,i-1)
1816  j0=max( 1,j-1)
1817  k0=max( 1,k-1)
1818  i1=min(lx1,i+1)
1819  j1=min(ly1,j+1)
1820  k1=min(lz1,k+1)
1821  do kk=k0,k1
1822  do jj=j0,j1
1823  do ii=i0,i1
1824 
1825  if (if3d) then
1826  dtmp = d(ii,jj,kk,e) + dist3d(
1827  $ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e),zm1(ii,jj,kk,e)
1828  $ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e),zm1(i ,j ,k ,e))
1829  else
1830  dtmp = d(ii,jj,kk,e) + dist2d(
1831  $ xm1(ii,jj,kk,e),ym1(ii,jj,kk,e)
1832  $ ,xm1(i ,j ,k ,e),ym1(i ,j ,k ,e))
1833  endif
1834 
1835  if (dtmp.lt.d(i,j,k,e)) then
1836  d(i,j,k,e) = dtmp
1837  nchange = nchange+1
1838  dmax = max(dmax,d(i,j,k,e))
1839  endif
1840  enddo
1841  enddo
1842  enddo
1843 
1844  enddo
1845  enddo
1846  enddo
1847  enddo
1848  call fgslib_gs_op(gsh_fld(ifld),d,1,3,0) ! min over all elements
1849  nchange = iglsum(nchange,1)
1850  dmax = glmax(dmax,1)
1851  if (nio.eq.0.and.loglevel.ge.3) write(6,1) ipass,nchange,dmax,b
1852  1 format(i9,i12,1pe12.4,' max distance b: ',a3)
1853  if (nchange.eq.0) goto 1000
1854  enddo
1855  1000 return
1856  end
1857 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
gop
subroutine gop(x, w, op, n)
Definition: comm_mpi.f:201
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
dnekclock
real *8 function dnekclock()
Definition: comm_mpi.f:393
facev
subroutine facev(a, ie, iface, val, nx, ny, nz)
Definition: connect1.f:1077
dsop
subroutine dsop(u, op, nx, ny, nz)
Definition: dssum.f:101
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
dot
real function dot(V1, V2, N)
Definition: genxyz.f:885
axhelm
subroutine axhelm(au, u, helm1, helm2, imesh, isd)
Definition: hmholtz.f:73
dsavg
subroutine dsavg(u)
Definition: ic.f:1878
capit
subroutine capit(lettrs, n)
Definition: ic.f:1648
geom_reset
subroutine geom_reset(icall)
Definition: ic.f:1802
cadd
subroutine cadd(a, const, n)
Definition: math.f:326
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
ascol5
subroutine ascol5(a, b, c, d, e, n)
Definition: math.f:153
sub2
subroutine sub2(a, b, n)
Definition: math.f:164
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
iglsum
function iglsum(a, n)
Definition: math.f:926
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
glamin
real function glamin(a, n)
Definition: math.f:887
blank
subroutine blank(A, N)
Definition: math.f:19
add3
subroutine add3(a, b, c, n)
Definition: math.f:644
col4
subroutine col4(a, b, c, d, n)
Definition: math.f:120
vlsum
real function vlsum(vec, n)
Definition: math.f:417
vlsc2
real function vlsc2(x, y, n)
Definition: math.f:739
subcol3
subroutine subcol3(a, b, c, n)
Definition: math.f:186
glsum
function glsum(x, n)
Definition: math.f:861
subcol4
subroutine subcol4(a, b, c, d, n)
Definition: math.f:197
sub3
subroutine sub3(a, b, c, n)
Definition: math.f:175
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
glamax
real function glamax(a, n)
Definition: math.f:874
chcopy
subroutine chcopy(a, b, n)
Definition: math.f:281
cfill
subroutine cfill(a, b, n)
Definition: math.f:244
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
glmax
function glmax(a, n)
Definition: math.f:960
chsign
subroutine chsign(a, n)
Definition: math.f:305
setupds_center
subroutine setupds_center(gs_handle, nx, ny, nz, nel, melg, vertex, glo_num)
Definition: meshsmoother.f:1228
genbackupmesh
subroutine genbackupmesh
Definition: meshsmoother.f:301
meshsmoother
subroutine meshsmoother
Definition: meshsmoother.f:2
chkjacinv
subroutine chkjacinv(jac, n, iel, X, Y, Z, ND, IERR)
Definition: meshsmoother.f:1510
int_fine_to_coarse_3d
subroutine int_fine_to_coarse_3d(u1, u2, n1, n2)
Definition: meshsmoother.f:213
smoothmesh
subroutine smoothmesh(mtyp, nouter, nlap, nopt, nbc, dcbc, idftyp, alpha, beta)
Definition: meshsmoother.f:28
int_coarse_to_fine_2d
subroutine int_coarse_to_fine_2d(u1, u2, n1, n2)
Definition: meshsmoother.f:245
my_mv_mesh
subroutine my_mv_mesh(dx, dy, dz)
Definition: meshsmoother.f:1634
restbackmesh
subroutine restbackmesh
Definition: meshsmoother.f:321
gen_int_lx1_to_3
subroutine gen_int_lx1_to_3
Definition: meshsmoother.f:160
disfun
subroutine disfun(dis, funtype, alpha, beta, dcbc, nbc)
Definition: meshsmoother.f:620
sm_cheap_dist
subroutine sm_cheap_dist(d, ifld, b)
Definition: meshsmoother.f:1772
laplaceh
subroutine laplaceh(name, u, h1, h2, mask, mult, ifld, tli, maxi, approx, napprox)
Definition: meshsmoother.f:1712
fixedgs
subroutine fixedgs(xyz, nx, ny, nz)
Definition: meshsmoother.f:863
get_len
subroutine get_len(val, x, y, z, siz)
Definition: meshsmoother.f:1083
dsavg_general
subroutine dsavg_general(u, mlt, gshl)
Definition: meshsmoother.f:1335
getglobsum
subroutine getglobsum(xe, ye, ze, mlt, gshl, siz, opt, fl)
Definition: meshsmoother.f:598
xmctoxm1
subroutine xmctoxm1
Definition: meshsmoother.f:1609
masklayers
subroutine masklayers(nodmask, nlayers)
Definition: meshsmoother.f:734
col3c
subroutine col3c(a, b, c, d, n)
Definition: meshsmoother.f:1412
restbndrlayer
subroutine restbndrlayer(dx, dy, dz, dis, mltc, gshlc)
Definition: meshsmoother.f:684
x8toxm
subroutine x8toxm(xd, x8)
Definition: meshsmoother.f:1380
int_fine_to_coarse_2d
subroutine int_fine_to_coarse_2d(u1, u2, n1, n2)
Definition: meshsmoother.f:193
xyzrstc
subroutine xyzrstc(xrmc, yrmc, zrmc, xsmc, ysmc, zsmc, XTMc, YTMc, ZTMc, IFAXIS)
Definition: meshsmoother.f:1532
optcg
subroutine optcg(itmax, nodmask, mlt, gshl, dis, opt, optinv, mltc, gshlc)
Definition: meshsmoother.f:347
get_nodscale
subroutine get_nodscale(scalek, x8, y8, z8, gshl)
Definition: meshsmoother.f:1110
xmtox8
subroutine xmtox8(xd, x8)
Definition: meshsmoother.f:1348
fastlap
subroutine fastlap(iter, nodmask, mlt, gshl, dis, lapinv, mltc, gshlc)
Definition: meshsmoother.f:1244
dolsalpha
subroutine dolsalpha(xe, ye, ze, sk, alpha, siz, opt, mlt, gshl, nodmask, iter)
Definition: meshsmoother.f:476
genmask
subroutine genmask(nodmask, mlt, gshl, mltc, gshlc)
Definition: meshsmoother.f:1153
int_coarse_to_fine_3d
subroutine int_coarse_to_fine_3d(u1, u2, n1, n2)
Definition: meshsmoother.f:269
glmapm1chkinv
subroutine glmapm1chkinv(kerr)
Definition: meshsmoother.f:1422
xm1toxmc
subroutine xm1toxmc
Definition: meshsmoother.f:1584
get_jac
subroutine get_jac(val, x, y, z, siz, el, dire)
Definition: meshsmoother.f:989
fixcurs
subroutine fixcurs(mltc, gshlc)
Definition: meshsmoother.f:777
gradf
subroutine gradf(f2, dfdx, x8, y8, z8, mlt, gshl, siz, opt, h)
Definition: meshsmoother.f:924
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
opadd2
subroutine opadd2(a1, a2, a3, b1, b2, b3)
Definition: navier1.f:2350
opsub3
subroutine opsub3(a1, a2, a3, b1, b2, b3, c1, c2, c3)
Definition: navier1.f:2370
project1
subroutine project1(b, n, rvar, ivar, h1, h2, msk, w, ifwt, ifvec, name6)
Definition: navier4.f:637
hmhzpf
subroutine hmhzpf(name, u, r, h1, h2, mask, mult, imesh, tli, maxit, isd, bi)
Definition: navier4.f:514
project2
subroutine project2(x, n, rvar, ivar, h1, h2, msk, w, ifwt, ifvec, name6)
Definition: navier4.f:1130
domain_size
subroutine domain_size(xmin, xmax, ymin, ymax, zmin, zmax)
Definition: navier5.f:2948
dist2d
function dist2d(a, b, x, y)
Definition: navier5.f:2938
fix_geom
subroutine fix_geom
Definition: navier5.f:2322
dist3d
function dist3d(a, b, c, x, y, z)
Definition: navier5.f:2928
gh_face_extend
subroutine gh_face_extend(x, zg, n, gh_type, e, v)
Definition: navier5.f:2422
set_vert
subroutine set_vert(glo_num, ngv, nx, nel, vertex, ifcenter)
Definition: navier8.f:4
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
dgll
subroutine dgll(D, DT, Z, NZ, lzd)
Definition: speclib.f:801
cmult2
subroutine cmult2(A, B, CONST, N)
Definition: subs2.f:341