KTH framework for Nek5000 toolboxes; testing version  0.0.1
genxyz.f
Go to the documentation of this file.
1 c-----------------------------------------------------------------------
2  subroutine arcsrf(xml,yml,zml,nxl,nyl,nzl,ie,isid)
3  include 'SIZE'
4  include 'GEOM'
5  include 'INPUT'
6  include 'TOPOL'
7  include 'WZ'
8 C
9 C ....note..... CTMP1 is used in this format in several subsequent routines
10 C
11  COMMON /ctmp1/ h(lx1,3,2),xcrved(lx1),ycrved(ly1),zcrved(lz1)
12  $ , zgml(lx1,3),work(3,lx1,lz1)
13  dimension xml(nxl,nyl,nzl,1),yml(nxl,nyl,nzl,1),zml(nxl,nyl,nzl,1)
14  LOGICAL IFGLJ
15 C
16  ifglj = .false.
17  IF (ifaxis .AND. ifrzer(ie) .AND. (isid.EQ.2 .OR. isid.EQ.4))
18  $ifglj = .true.
19 C
20  pt1x = xc(isid,ie)
21  pt1y = yc(isid,ie)
22  IF(isid.EQ.4) THEN
23  pt2x = xc(1,ie)
24  pt2y = yc(1,ie)
25  ELSE IF(isid.EQ.8) THEN
26  pt2x = xc(5,ie)
27  pt2y = yc(5,ie)
28  ELSE
29  pt2x = xc(isid+1,ie)
30  pt2y = yc(isid+1,ie)
31  ENDIF
32 C
33 C Find slope of perpendicular
34  radius=curve(1,isid,ie)
35  gap=sqrt( (pt1x-pt2x)**2 + (pt1y-pt2y)**2 )
36  IF (abs(2.0*radius).LE.gap*1.00001) THEN
37  write(6,10) radius,isid,ie,gap
38  10 FORMAT(//,2x,'ERROR: Too small a radius (',g11.3
39  $ ,') specified for side',i2,' of element',i4,': '
40  $ ,g11.3,/,2x,'ABORTING during mesh generation.')
41  call exitt
42  ENDIF
43  xs = pt2y-pt1y
44  ys = pt1x-pt2x
45 C Make length Radius
46  xys=sqrt(xs**2+ys**2)
47 C Find Center
48  dtheta = abs(asin(0.5*gap/radius))
49  pt12x = (pt1x + pt2x)/2.0
50  pt12y = (pt1y + pt2y)/2.0
51  xcenn = pt12x - xs/xys * radius*cos(dtheta)
52  ycenn = pt12y - ys/xys * radius*cos(dtheta)
53  theta0 = atan2((pt12y-ycenn),(pt12x-xcenn))
54  IF (ifglj) THEN
55  fac = sign(1.0,radius)
56  theta1 = theta0 - fac*dtheta
57  theta2 = theta0 + fac*dtheta
58  ENDIF
59 C Compute perturbation of geometry
60  isid1 = mod1(isid,4)
61  IF (ifglj) THEN
62  i1 = isid/2
63  i2 = 2 - isid/4
64  DO 15 iy=1,nyl
65  ang = h(iy,2,i1)*theta1 + h(iy,2,i2)*theta2
66  xcrved(iy)=xcenn + abs(radius)*cos(ang)
67  $ - (h(iy,2,i1)*pt1x + h(iy,2,i2)*pt2x)
68  ycrved(iy)=ycenn + abs(radius) * sin(ang)
69  $ - (h(iy,2,i1)*pt1y + h(iy,2,i2)*pt2y)
70  15 CONTINUE
71  ELSE
72  DO 20 ix=1,nxl
73  ixt=ix
74  IF (isid1.GT.2) ixt=nxl+1-ix
75  r=zgml(ix,1)
76  IF (radius.LT.0.0) r=-r
77  xcrved(ixt) = xcenn + abs(radius) * cos(theta0 + r*dtheta)
78  $ - ( h(ix,1,1)*pt1x + h(ix,1,2)*pt2x )
79  ycrved(ixt) = ycenn + abs(radius) * sin(theta0 + r*dtheta)
80  $ - ( h(ix,1,1)*pt1y + h(ix,1,2)*pt2y )
81  20 CONTINUE
82  ENDIF
83 C Points all set, add perturbation to current mesh.
84  isid1 = mod1(isid,4)
85  isid1 = eface1(isid1)
86  izt = (isid-1)/4+1
87  iyt = isid1-2
88  ixt = isid1
89  IF (isid1.LE.2) THEN
90  CALL addtnsr(xml(1,1,1,ie),h(1,1,ixt),xcrved,h(1,3,izt)
91  $ ,nxl,nyl,nzl)
92  CALL addtnsr(yml(1,1,1,ie),h(1,1,ixt),ycrved,h(1,3,izt)
93  $ ,nxl,nyl,nzl)
94  ELSE
95  CALL addtnsr(xml(1,1,1,ie),xcrved,h(1,2,iyt),h(1,3,izt)
96  $ ,nxl,nyl,nzl)
97  CALL addtnsr(yml(1,1,1,ie),ycrved,h(1,2,iyt),h(1,3,izt)
98  $ ,nxl,nyl,nzl)
99  ENDIF
100  return
101  end
102 c-----------------------------------------------------------------------
103  subroutine defsrf(xml,yml,zml,nxl,nyl,nzl,ie,iface1,ccv)
104  include 'SIZE'
105  include 'TOPOL'
106  include 'GEOM'
107  include 'WZ'
108  COMMON /ctmp1/ h(lx1,3,2),xcrved(lx1),ycrved(ly1),zcrved(lz1)
109  $ , zgml(lx1,3),work(3,lx1,lz1)
110 C
111  dimension xml(nxl,nyl,nzl,1),yml(nxl,nyl,nzl,1),zml(nxl,nyl,nzl,1)
112  dimension x1(3),x2(3),x3(3),dx(3)
113  dimension iopp(3),nxx(3)
114  CHARACTER*1 CCV
115 C
116  CALL dsset(nxl,nyl,nzl)
117  iface = eface1(iface1)
118  js1 = skpdat(1,iface)
119  jf1 = skpdat(2,iface)
120  jskip1 = skpdat(3,iface)
121  js2 = skpdat(4,iface)
122  jf2 = skpdat(5,iface)
123  jskip2 = skpdat(6,iface)
124 C
125  iopp(1) = nxl-1
126  iopp(2) = nxl*(nyl-1)
127  iopp(3) = nxl*nyl*(nzl-1)
128  nxx(1) = nxl
129  nxx(2) = nyl
130  nxx(3) = nzl
131  idir = 2*mod(iface,2) - 1
132  ifc2 = (iface+1)/2
133  delt = 0.0
134 C DELT = SIDE(4,IFACE,IE)
135 C
136 C Compute surface deflection and perturbation due to face IFACE
137 C
138  DO 200 j2=js2,jf2,jskip2
139  DO 200 j1=js1,jf1,jskip1
140  jopp = j1 + iopp(ifc2)*idir
141  x2(1) = xml(j1,j2,1,ie)
142  x2(2) = yml(j1,j2,1,ie)
143  x2(3) = zml(j1,j2,1,ie)
144  x1(1) = xml(jopp,j2,1,ie)
145  x1(2) = yml(jopp,j2,1,ie)
146  x1(3) = zml(jopp,j2,1,ie)
147  CALL intrsc(x3,x2,x1,delt,ie,iface1)
148 C
149  dx(1) = x3(1)-x2(1)
150  dx(2) = x3(2)-x2(2)
151  dx(3) = x3(3)-x2(3)
152 C
153  nxs = nxx(ifc2)
154  joff = (j1-jopp)/(nxs-1)
155  DO 100 ix = 2,nxs
156  j = jopp + joff*(ix-1)
157  zeta = 0.5*(zgml(ix,ifc2) + 1.0)
158  xml(j,j2,1,ie) = xml(j,j2,1,ie)+dx(1)*zeta
159  yml(j,j2,1,ie) = yml(j,j2,1,ie)+dx(2)*zeta
160  zml(j,j2,1,ie) = zml(j,j2,1,ie)+dx(3)*zeta
161  100 CONTINUE
162  200 CONTINUE
163 C
164  return
165  end
166 c-----------------------------------------------------------------------
167  subroutine intrsc(x3,x2,x1,delt,ie,iface)
168 C
169  dimension x1(3),x2(3),x3(3)
170  COMMON /srfcei/ iel,ifce
171  COMMON /srfcer/ x0(3),dx(3)
172  COMMON /srfcel/ succes
173  LOGICAL SUCCES
174  COMMON /tolrnc/ tol
175 C
176 C Load parameters for surface function FNC
177 C
178  iel = ie
179  ifce = iface
180  x0(1) = x1(1)
181  x0(2) = x1(2)
182  x0(3) = x1(3)
183  dx(1) = x2(1)-x1(1)
184  dx(2) = x2(2)-x1(2)
185  dx(3) = x2(3)-x1(3)
186  dist = sqrt( dx(1)**2 + dx(2)**2 + dx(3)**2 )
187 C
188 C Initial guess for bracket is given by size of element face (DELT).
189 C
190  eta2 = 1.0
191  eta1 = eta2 - delt/dist
192  eta1 = max(eta1,0.0)
193  CALL zbrac(eta1,eta2,succes)
194 C
195  tol = 1.0e-5
196  tolsrf = tol*(eta2-eta1)
197  IF (succes) eta3 = zbrent(eta1,eta2,tolsrf)
198 C
199  x3(1) = x0(1) + dx(1)*eta3
200  x3(2) = x0(2) + dx(2)*eta3
201  x3(3) = x0(3) + dx(3)*eta3
202 C
203  return
204  end
205 c-----------------------------------------------------------------------
206  subroutine zbrac(x1,x2,succes)
207 C
208 C Given a function FNC and an initial guess (X1,X2), the routine
209 C expands the range geometrically until a root is bracketed by the
210 C returned range (X1,X2) [SUCCES=.TRUE.], or until the range becomes
211 C unacceptably large [SUCCES=.FALSE.].
212 C ( Numerical Recipes, p. 245; pff 9 Aug 1989 09:00:20 )
213 C
214  parameter(factor=1.08,ntry=50)
215  LOGICAL SUCCES
216 C
217  succes = .true.
218 C
219  IF (x1.EQ.x2) x1 = .99*x1
220  IF (x1.EQ.0.0) x1 = 1.0e-04
221 C
222  f1 = fnc(x1)
223  f2 = fnc(x2)
224  DO 100 j=1,ntry
225  IF (f1*f2.LT.0.0) return
226  IF (abs(f1).LT.abs(f2)) THEN
227  x1 = x1 + factor*(x1-x2)
228  f1 = fnc(x1)
229  ELSE
230  x2 = x2 + factor*(x2-x1)
231  f2 = fnc(x2)
232  ENDIF
233  100 CONTINUE
234  succes = .false.
235  return
236  end
237  FUNCTION zbrent(X1,X2,TOL)
238 C
239 C Using the Van Wijngaarden-Dekker-Brent Method, find the root
240 C of a function FNC known to lie between X1 and X2. The root,
241 C returned as ZBRENT, will be refined until its accuracy is TOL.
242 C
243  parameter(itmax=100,eps=3.0e-8)
244 C
245  a = x1
246  b = x2
247  fa = fnc(a)
248  fb = fnc(b)
249  IF (fb*fa.GT.0.0) GOTO 9000
250  fc=fb
251 C
252  DO 1000 iter=1,itmax
253  IF (fb*fc.GT.0.0) THEN
254  c = a
255  fc = fa
256  d = b-a
257  e = d
258  ENDIF
259  IF (abs(fc).LT.abs(fb)) THEN
260  a = b
261  b = c
262  c = a
263  fa = fb
264  fb = fc
265  fc = fa
266  ENDIF
267  tol1 = 2.0*eps*abs(b)+0.5*tol
268  xm = 0.5*(c-b)
269  IF (abs(xm).LE.tol1.OR.fb.EQ.0.0) THEN
270  zbrent = b
271  return
272  ENDIF
273  IF (abs(e).GT.tol1.AND. abs(fa).GT.abs(fb)) THEN
274 C
275 C Attempt inverse quadratic interpolation
276 C
277  s=fb/fa
278  IF (a.EQ.c) THEN
279  p=2.0*xm*s
280  q=1.0-s
281  ELSE
282  q=fa/fc
283  r=fb/fc
284  p=s*( 2.0*xm*q*(q-r) - (b-a)*(r-1.0) )
285  q=(q-1.0)*(r-1.0)*(s-1.0)
286  ENDIF
287 C
288 C Check whether in bounds...
289 C
290  IF (p.GT.0.0) q = -q
291  p = abs(p)
292  IF (2.0*p.LT.min(3.0*xm*q-abs(tol1*q),abs(e*q))) THEN
293 C Accept interpolation.
294  e=d
295  d=p/q
296  ELSE
297 C Interpolation failed, use bisection.
298  d=xm
299  e=d
300  ENDIF
301  ELSE
302 C Bounds decreasing too slowly, use bisection.
303  d=xm
304  e=d
305  ENDIF
306  a=b
307  fa=fb
308  IF (abs(d).GT.tol1) THEN
309 C Evaluate new trial root
310  b=b+d
311  ELSE
312  b=b+sign(tol1,xm)
313  ENDIF
314  fb=fnc(b)
315  1000 CONTINUE
316 C
317 C
318  9000 CONTINUE
319  write(6 ,*) 'Exceeding maximum number of iterations.'
320 C WRITE(21,*) 'Exceeding maximum number of iterations.'
321  zbrent=b
322  return
323  end
324  FUNCTION fnc(ETA)
325  include 'SIZE'
326  include 'INPUT'
327  COMMON /srfcei/ iel,ifce
328  COMMON /srfcer/ x0(3),dx(3)
329  COMMON /srfcel/ succes
330  dimension x3(3)
331  integer icalld
332  save icalld
333  data icalld /0/
334 C
335  IF (ccurve(ifce,iel).EQ.'s') THEN
336  xctr = curve(1,ifce,iel)
337  yctr = curve(2,ifce,iel)
338  zctr = curve(3,ifce,iel)
339  radius = curve(4,ifce,iel)
340  x = x0(1) + dx(1)*eta - xctr
341  y = x0(2) + dx(2)*eta - yctr
342  z = x0(3) + dx(3)*eta - zctr
343 C
344  fnc = (radius**2 - (x**2+y**2+z**2))/(radius**2)
345 C
346  ENDIF
347 C
348  return
349  end
350 c-----------------------------------------------------------------------
351  subroutine setdef
352 C-------------------------------------------------------------------
353 C
354 C Set up deformed element logical switches
355 C
356 C-------------------------------------------------------------------
357  include 'SIZE'
358  include 'INPUT'
359  dimension xcc(8),ycc(8),zcc(8)
360  dimension indx(8)
361  REAL VEC(3,12)
362  LOGICAL IFVCHK
363 C
364  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
365  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
366 C
367 C Corner notation:
368 C
369 C 4+-----+3 ^ Y
370 C / /| |
371 C / / | |
372 C 8+-----+7 +2 +----> X
373 C | | / /
374 C | |/ /
375 C 5+-----+6 Z
376 C
377 C
378  DO 10 ie=1,nelt
379  ifdfrm(ie)=.true.
380  10 CONTINUE
381 C
382  IF (ifmvbd) return
383 c
384 c Force IFDFRM=.true. for all elements (for timing purposes only)
385 c
386  IF (param(59).ne.0.and.nio.eq.0)
387  $ write(6,*) 'NOTE: All elements deformed , param(59) ^=0'
388  IF (param(59).ne.0) return
389 C
390 C Check against cases which won't allow for savings in HMHOLTZ
391 C
392  indx(1)=1
393  indx(2)=2
394  indx(3)=4
395  indx(4)=3
396  indx(5)=5
397  indx(6)=6
398  indx(7)=8
399  indx(8)=7
400 C
401 C Check for deformation (rotation is acceptable).
402 C
403  DO 500 ie=1,nelt
404 C
405  call rzero(vec,36)
406  IF (if3d) THEN
407  DO 100 iedg=1,8
408  IF(ccurve(iedg,ie).NE.' ') THEN
409  ifdfrm(ie)=.true.
410  GOTO 500
411  ENDIF
412  100 CONTINUE
413 C
414  DO 105 i=1,8
415  xcc(i)=xc(indx(i),ie)
416  ycc(i)=yc(indx(i),ie)
417  zcc(i)=zc(indx(i),ie)
418  105 CONTINUE
419 C
420  DO 110 i=1,4
421  vec(1,i)=xcc(2*i)-xcc(2*i-1)
422  vec(2,i)=ycc(2*i)-ycc(2*i-1)
423  vec(3,i)=zcc(2*i)-zcc(2*i-1)
424  110 CONTINUE
425 C
426  i1=4
427  DO 120 i=0,1
428  DO 120 j=0,1
429  i1=i1+1
430  i2=4*i+j+3
431  vec(1,i1)=xcc(i2)-xcc(i2-2)
432  vec(2,i1)=ycc(i2)-ycc(i2-2)
433  vec(3,i1)=zcc(i2)-zcc(i2-2)
434  120 CONTINUE
435 C
436  i1=8
437  DO 130 i=5,8
438  i1=i1+1
439  vec(1,i1)=xcc(i)-xcc(i-4)
440  vec(2,i1)=ycc(i)-ycc(i-4)
441  vec(3,i1)=zcc(i)-zcc(i-4)
442  130 CONTINUE
443 C
444  DO 140 i=1,12
445  veclen = vec(1,i)**2 + vec(2,i)**2 + vec(3,i)**2
446  veclen = sqrt(veclen)
447  vec(1,i)=vec(1,i)/veclen
448  vec(2,i)=vec(2,i)/veclen
449  vec(3,i)=vec(3,i)/veclen
450  140 CONTINUE
451 C
452 C Check the dot product of the adjacent edges to see that it is zero.
453 C
454  ifdfrm(ie)=.false.
455  IF ( ifvchk(vec,1,5, 9) ) ifdfrm(ie)=.true.
456  IF ( ifvchk(vec,1,6,10) ) ifdfrm(ie)=.true.
457  IF ( ifvchk(vec,2,5,11) ) ifdfrm(ie)=.true.
458  IF ( ifvchk(vec,2,6,12) ) ifdfrm(ie)=.true.
459  IF ( ifvchk(vec,3,7, 9) ) ifdfrm(ie)=.true.
460  IF ( ifvchk(vec,3,8,10) ) ifdfrm(ie)=.true.
461  IF ( ifvchk(vec,4,7,11) ) ifdfrm(ie)=.true.
462  IF ( ifvchk(vec,4,8,12) ) ifdfrm(ie)=.true.
463 C
464 C Check the 2D case....
465 C
466  ELSE
467 C
468  DO 200 iedg=1,4
469  IF(ccurve(iedg,ie).NE.' ') THEN
470  ifdfrm(ie)=.true.
471  GOTO 500
472  ENDIF
473  200 CONTINUE
474 C
475  DO 205 i=1,4
476  xcc(i)=xc(indx(i),ie)
477  ycc(i)=yc(indx(i),ie)
478  205 CONTINUE
479 C
480  vec(1,1)=xcc(2)-xcc(1)
481  vec(1,2)=xcc(4)-xcc(3)
482  vec(1,3)=xcc(3)-xcc(1)
483  vec(1,4)=xcc(4)-xcc(2)
484  vec(1,5)=0.0
485  vec(2,1)=ycc(2)-ycc(1)
486  vec(2,2)=ycc(4)-ycc(3)
487  vec(2,3)=ycc(3)-ycc(1)
488  vec(2,4)=ycc(4)-ycc(2)
489  vec(2,5)=0.0
490 C
491  DO 220 i=1,4
492  veclen = vec(1,i)**2 + vec(2,i)**2
493  veclen = sqrt(veclen)
494  vec(1,i)=vec(1,i)/veclen
495  vec(2,i)=vec(2,i)/veclen
496  220 CONTINUE
497 C
498 C Check the dot product of the adjacent edges to see that it is zero.
499 C
500  ifdfrm(ie)=.false.
501  IF ( ifvchk(vec,1,3,5) ) ifdfrm(ie)=.true.
502  IF ( ifvchk(vec,1,4,5) ) ifdfrm(ie)=.true.
503  IF ( ifvchk(vec,2,3,5) ) ifdfrm(ie)=.true.
504  IF ( ifvchk(vec,2,4,5) ) ifdfrm(ie)=.true.
505  ENDIF
506  500 CONTINUE
507  return
508  end
509  LOGICAL FUNCTION ifvchk(VEC,I1,I2,I3)
510 C
511 C Take the dot product of the three components of VEC to see if it's zero.
512 C
513  dimension vec(3,12)
514  LOGICAL iftmp
515 C
516  iftmp=.false.
517  epsm=1.0e-06
518 C
519  dot1=vec(1,i1)*vec(1,i2)+vec(2,i1)*vec(2,i2)+vec(3,i1)*vec(3,i2)
520  dot2=vec(1,i2)*vec(1,i3)+vec(2,i2)*vec(2,i3)+vec(3,i2)*vec(3,i3)
521  dot3=vec(1,i1)*vec(1,i3)+vec(2,i1)*vec(2,i3)+vec(3,i1)*vec(3,i3)
522 C
523  dot1=abs(dot1)
524  dot2=abs(dot2)
525  dot3=abs(dot3)
526  dot=dot1+dot2+dot3
527  IF (dot.GT.epsm) iftmp=.true.
528 C
529  ifvchk=iftmp
530  return
531  end
532 c-----------------------------------------------------------------------
533  subroutine gencoor (xm3,ym3,zm3)
534 C-----------------------------------------------------------------------
535 C
536 C Generate xyz coordinates for all elements.
537 C Velocity formulation : mesh 3 is used
538 C Stress formulation : mesh 1 is used
539 C
540 C-----------------------------------------------------------------------
541  include 'SIZE'
542  include 'GEOM'
543  include 'INPUT'
544  dimension xm3(lx3,ly3,lz3,1),ym3(lx3,ly3,lz3,1),zm3(lx3,ly3,lz3,1)
545 C
546 C Select appropriate mesh
547 C
548  IF ( ifgmsh3 ) THEN
549  CALL genxyz (xm3,ym3,zm3,lx3,ly3,lz3)
550  ELSE
551  CALL genxyz (xm1,ym1,zm1,lx1,ly1,lz1)
552  ENDIF
553 C
554  return
555  end
556 c-----------------------------------------------------------------------
557  subroutine genxyz (xml,yml,zml,nxl,nyl,nzl)
558 C
559  include 'SIZE'
560  include 'WZ'
561  include 'GEOM'
562  include 'TOPOL'
563  include 'INPUT'
564  include 'PARALLEL'
565 
566  real xml(nxl,nyl,nzl,1),yml(nxl,nyl,nzl,1),zml(nxl,nyl,nzl,1)
567 
568 C Note : CTMP1 is used in this format in several subsequent routines
569  common /ctmp1/ h(lx1,3,2),xcrved(lx1),ycrved(ly1),zcrved(lz1)
570  $ , zgml(lx1,3),work(3,lx1,lz1)
571 
572  parameter(ldw=2*lx1*ly1*lz1)
573  common /ctmp0/ w(ldw)
574 
575  character*1 ccv
576 
577 c Initialize geometry arrays with bi- triquadratic deformations
578  call linquad(xml,yml,zml,nxl,nyl,nzl)
579 
580 
581  do ie=1,nelt
582 
583  call setzgml (zgml,ie,nxl,nyl,nzl,ifaxis)
584  call sethmat (h,zgml,nxl,nyl,nzl)
585 
586 c Deform surfaces - general 3D deformations
587 c - extruded geometry deformations
588  nfaces = 2*ldim
589  do iface=1,nfaces
590  ccv = ccurve(iface,ie)
591  if (ccv.eq.'s')
592  $ call sphsrf(xml,yml,zml,iface,ie,nxl,nyl,nzl,work)
593  if (ccv.eq.'e')
594  $ call gensrf(xml,yml,zml,iface,ie,nxl,nyl,nzl,zgml)
595  enddo
596 
597  do isid=1,8
598  ccv = ccurve(isid,ie)
599  if (ccv.eq.'C') call arcsrf(xml,yml,zml,nxl,nyl,nzl,ie,isid)
600  enddo
601 
602  enddo
603 
604 c call user_srf(xml,yml,zml,nxl,nyl,nzl)
605 c call opcopy(xm1,ym1,zm1,xml,yml,zml)
606 c call outpost(xml,yml,zml,xml,yml,' ')
607 c call exitt
608 C
609  return
610  end
611 c-----------------------------------------------------------------------
612  subroutine sethmat(h,zgml,nxl,nyl,nzl)
613 
614  include 'SIZE'
615  include 'INPUT' ! if3d
616 
617  real h(lx1,3,2),zgml(lx1,3)
618 
619  do 10 ix=1,nxl
620  h(ix,1,1)=(1.0-zgml(ix,1))*0.5
621  h(ix,1,2)=(1.0+zgml(ix,1))*0.5
622  10 continue
623  do 20 iy=1,nyl
624  h(iy,2,1)=(1.0-zgml(iy,2))*0.5
625  h(iy,2,2)=(1.0+zgml(iy,2))*0.5
626  20 continue
627  if (if3d) then
628  do 30 iz=1,nzl
629  h(iz,3,1)=(1.0-zgml(iz,3))*0.5
630  h(iz,3,2)=(1.0+zgml(iz,3))*0.5
631  30 continue
632  else
633  call rone(h(1,3,1),nzl)
634  call rone(h(1,3,2),nzl)
635  endif
636 
637  return
638  end
639 c-----------------------------------------------------------------------
640  subroutine setzgml (zgml,e,nxl,nyl,nzl,ifaxl)
641 
642  include 'SIZE'
643  include 'WZ'
644  include 'GEOM'
645 
646  real zgml(lx1,3)
647  integer e
648  logical ifaxl
649 
650  call rzero (zgml,3*lx1)
651 
652 
653  if (nxl.eq.3 .and. .not. ifaxl) then
654  do k=1,3
655  zgml(1,k) = -1
656  zgml(2,k) = 0
657  zgml(3,k) = 1
658  enddo
659  elseif (ifgmsh3.and.nxl.eq.lx3) then
660  call copy(zgml(1,1),zgm3(1,1),lx3)
661  call copy(zgml(1,2),zgm3(1,2),ly3)
662  call copy(zgml(1,3),zgm3(1,3),lz3)
663  if (ifaxl .and. ifrzer(e)) call copy(zgml(1,2),zam3,ly3)
664  elseif (nxl.eq.lx1) then
665  call copy(zgml(1,1),zgm1(1,1),lx1)
666  call copy(zgml(1,2),zgm1(1,2),ly1)
667  call copy(zgml(1,3),zgm1(1,3),lz1)
668  if (ifaxl .and. ifrzer(e)) call copy(zgml(1,2),zam1,ly1)
669  else
670  call exitti('ABORT setzgml! $',nxl)
671  endif
672 
673  return
674  end
675 c-----------------------------------------------------------------------
676  subroutine sphsrf(xml,yml,zml,ifce,ie,nx,ny,nz,xysrf)
677 C
678 C 5 Aug 1988 19:29:52
679 C
680 C Program to generate spherical shell elements for NEKTON
681 C input. Paul F. Fischer
682 C
683  include 'SIZE'
684  include 'INPUT'
685  include 'WZ'
686  include 'TOPOL'
687  dimension xml(nx,ny,nz,1),yml(nx,ny,nz,1),zml(nx,ny,nz,1)
688  dimension xysrf(3,nx,nz)
689 C
690  COMMON /ctmp1/ h(lx1,3,2),xcrved(lx1),ycrved(ly1),zcrved(lz1)
691  $ , zgml(lx1,3),work(3,lx1,lz1)
692  COMMON /ctmp0/ xcv(3,2,2),vn1(3),vn2(3)
693  $ ,x1(3),x2(3),x3(3),dx(3)
694  dimension iopp(3),nxx(3)
695 c
696 c
697 c These are representative nodes on a given face, and their opposites
698 c
699  integer cface(2,6)
700  save cface
701  data cface / 1,4 , 2,1 , 3,2 , 4,3 , 1,5 , 5,1 /
702  real vout(3),vsph(3)
703  logical ifconcv
704 c
705 C
706 C Determine geometric parameters
707 C
708  nxm1 = nx-1
709  nym1 = ny-1
710  nxy = nx*nz
711  nxy3 = 3*nx*nz
712  xctr = curve(1,ifce,ie)
713  yctr = curve(2,ifce,ie)
714  zctr = curve(3,ifce,ie)
715  radius = curve(4,ifce,ie)
716  iface = eface1(ifce)
717 C
718 C Generate (normalized) corner vectors XCV(1,i,j):
719 C
720  CALL crn3d(xcv,xc(1,ie),yc(1,ie),zc(1,ie),curve(1,ifce,ie),iface)
721 C
722 C Generate edge vectors on the sphere RR=1.0,
723 C for (r,s) = (-1,*),(1,*),(*,-1),(*,1)
724 C
725  CALL edg3d(xysrf,xcv(1,1,1),xcv(1,1,2), 1, 1, 1,ny,nx,ny)
726  CALL edg3d(xysrf,xcv(1,2,1),xcv(1,2,2),nx,nx, 1,ny,nx,ny)
727  CALL edg3d(xysrf,xcv(1,1,1),xcv(1,2,1), 1,nx, 1, 1,nx,ny)
728  CALL edg3d(xysrf,xcv(1,1,2),xcv(1,2,2), 1,nx,ny,ny,nx,ny)
729 C
730 C Generate intersection vectors for (i,j)
731 C
732 C quick check on sign of curvature: (pff , 12/08/00)
733 c
734 c
735  ivtx = cface(1,ifce)
736  ivto = cface(2,ifce)
737  vout(1) = xc(ivtx,ie)-xc(ivto,ie)
738  vout(2) = yc(ivtx,ie)-yc(ivto,ie)
739  vout(3) = zc(ivtx,ie)-zc(ivto,ie)
740 c
741  vsph(1) = xc(ivtx,ie)-xctr
742  vsph(2) = yc(ivtx,ie)-yctr
743  vsph(3) = zc(ivtx,ie)-zctr
744  ifconcv = .true.
745  sign = dot(vsph,vout,3)
746  if (sign.gt.0) ifconcv = .false.
747 c write(6,*) 'THIS IS SIGN:',sign
748 c
749  DO 200 j=2,nym1
750  CALL cross(vn1,xysrf(1,1,j),xysrf(1,nx,j))
751  DO 200 i=2,nxm1
752  CALL cross(vn2,xysrf(1,i,1),xysrf(1,i,ny))
753  if (ifconcv) then
754 c IF (IFACE.EQ.1.OR.IFACE.EQ.4.OR.IFACE.EQ.5) THEN
755  CALL cross(xysrf(1,i,j),vn2,vn1)
756  ELSE
757  CALL cross(xysrf(1,i,j),vn1,vn2)
758  ENDIF
759  200 CONTINUE
760 C
761 C Normalize all vectors to the unit sphere.
762 C
763  DO 300 i=1,nxy
764  CALL norm3d(xysrf(1,i,1))
765  300 CONTINUE
766 C
767 C Scale by actual radius
768 C
769  CALL cmult(xysrf,radius,nxy3)
770 C
771 C Add back the sphere center offset
772 C
773  DO 400 i=1,nxy
774  xysrf(1,i,1)=xysrf(1,i,1)+xctr
775  xysrf(2,i,1)=xysrf(2,i,1)+yctr
776  xysrf(3,i,1)=xysrf(3,i,1)+zctr
777  400 CONTINUE
778 C
779 C
780 C Transpose data, if necessary
781 C
782  IF (iface.EQ.1.OR.iface.EQ.4.OR.iface.EQ.5) THEN
783  DO 500 j=1 ,ny
784  DO 500 i=j+1,nx
785  tmp=xysrf(1,i,j)
786  xysrf(1,i,j)=xysrf(1,j,i)
787  xysrf(1,j,i)=tmp
788  tmp=xysrf(2,i,j)
789  xysrf(2,i,j)=xysrf(2,j,i)
790  xysrf(2,j,i)=tmp
791  tmp=xysrf(3,i,j)
792  xysrf(3,i,j)=xysrf(3,j,i)
793  xysrf(3,j,i)=tmp
794  500 CONTINUE
795  ENDIF
796 C
797 C Compute surface deflection and perturbation due to face IFACE
798 C
799  CALL dsset(nx,ny,nz)
800  js1 = skpdat(1,iface)
801  jf1 = skpdat(2,iface)
802  jskip1 = skpdat(3,iface)
803  js2 = skpdat(4,iface)
804  jf2 = skpdat(5,iface)
805  jskip2 = skpdat(6,iface)
806 C
807  iopp(1) = nx-1
808  iopp(2) = nx*(ny-1)
809  iopp(3) = nx*ny*(nz-1)
810  nxx(1) = nx
811  nxx(2) = ny
812  nxx(3) = nz
813  idir = 2*mod(iface,2) - 1
814  ifc2 = (iface+1)/2
815  delt = 0.0
816  i=0
817  DO 700 j2=js2,jf2,jskip2
818  DO 700 j1=js1,jf1,jskip1
819  i=i+1
820  jopp = j1 + iopp(ifc2)*idir
821  x2(1) = xml(j1,j2,1,ie)
822  x2(2) = yml(j1,j2,1,ie)
823  x2(3) = zml(j1,j2,1,ie)
824 C
825  dx(1) = xysrf(1,i,1)-x2(1)
826  dx(2) = xysrf(2,i,1)-x2(2)
827  dx(3) = xysrf(3,i,1)-x2(3)
828 C
829  nxs = nxx(ifc2)
830  joff = (j1-jopp)/(nxs-1)
831  DO 600 ix = 2,nxs
832  j = jopp + joff*(ix-1)
833  zeta = 0.5*(zgml(ix,ifc2) + 1.0)
834  xml(j,j2,1,ie) = xml(j,j2,1,ie)+dx(1)*zeta
835  yml(j,j2,1,ie) = yml(j,j2,1,ie)+dx(2)*zeta
836  zml(j,j2,1,ie) = zml(j,j2,1,ie)+dx(3)*zeta
837  600 CONTINUE
838  700 CONTINUE
839 C
840  return
841  end
842 c-----------------------------------------------------------------------
843  subroutine edg3d(xysrf,x1,x2,i1,i2,j1,j2,nx,ny)
844 C
845 C Generate XYZ vector along an edge of a surface.
846 C
847  include 'SIZE'
848  COMMON /ctmp1/ h(lx1,3,2),xcrved(lx1),ycrved(ly1),zcrved(lz1)
849  $ , zgml(lx1,3),work(3,lx1,lz1)
850 C
851  dimension xysrf(3,nx,ny)
852  dimension x1(3),x2(3)
853  REAL U1(3),U2(3),VN(3),B(3)
854 C
855 C Normalize incoming vectors
856 C
857  CALL copy (u1,x1,3)
858  CALL copy (u2,x2,3)
859  CALL norm3d (u1)
860  CALL norm3d (u2)
861 C
862 C Find normal to the plane and tangent to the curve.
863 C
864  CALL cross(vn,x1,x2)
865  CALL cross( b,vn,x1)
866  CALL norm3d (vn)
867  CALL norm3d (b)
868 C
869  ctheta = dot(u1,u2,3)
870  theta = acos(ctheta)
871 C
872  ij = 0
873  DO 200 j=j1,j2
874  DO 200 i=i1,i2
875  ij = ij + 1
876  thetap = 0.5*theta*(zgml(ij,1)+1.0)
877  ctp = cos(thetap)
878  stp = sin(thetap)
879  DO 200 iv = 1,3
880  xysrf(iv,i,j) = ctp*u1(iv) + stp*b(iv)
881  200 CONTINUE
882  return
883  end
884  REAL function dot(v1,v2,n)
885 C
886 C Compute Cartesian vector dot product.
887 C
888  dimension v1(n),v2(n)
889 C
890  sum = 0
891  DO 100 i=1,n
892  sum = sum + v1(i)*v2(i)
893  100 CONTINUE
894  dot = sum
895  return
896  end
897 c-----------------------------------------------------------------------
898  subroutine cross(v1,v2,v3)
899 C
900 C Compute Cartesian vector dot product.
901 C
902  dimension v1(3),v2(3),v3(3)
903 C
904  v1(1) = v2(2)*v3(3) - v2(3)*v3(2)
905  v1(2) = v2(3)*v3(1) - v2(1)*v3(3)
906  v1(3) = v2(1)*v3(2) - v2(2)*v3(1)
907 C
908  return
909  end
910 c-----------------------------------------------------------------------
911  subroutine norm3d(v1)
912 C
913 C Compute Cartesian vector dot product.
914 C
915  dimension v1(3)
916 C
917  vlngth = dot(v1,v1,3)
918  vlngth = sqrt(vlngth)
919  if (vlngth.gt.0) then
920  v1(1) = v1(1) / vlngth
921  v1(2) = v1(2) / vlngth
922  v1(3) = v1(3) / vlngth
923  endif
924 C
925  return
926  end
927 c-----------------------------------------------------------------------
928  subroutine crn3d(xcv,xc,yc,zc,curve,iface)
929  include 'SIZE'
930  include 'TOPOL'
931  dimension xcv(3,2,2),xc(8),yc(8),zc(8),curve(4)
932  dimension indvtx(4,6)
933  SAVE indvtx
934  DATA indvtx / 1,5,3,7 , 2,4,6,8 , 1,2,5,6
935  $ , 3,7,4,8 , 1,3,2,4 , 5,6,7,8 /
936 C
937  eps = 1.0e-4
938  xctr = curve(1)
939  yctr = curve(2)
940  zctr = curve(3)
941  radius = curve(4)
942 C
943  DO 10 i=1,4
944  j=indvtx(i,iface)
945  k=indx(j)
946  xcv(1,i,1)=xc(k)-xctr
947  xcv(2,i,1)=yc(k)-yctr
948  xcv(3,i,1)=zc(k)-zctr
949  10 CONTINUE
950 C
951 C Check to ensure that these points are indeed on the sphere.
952 C
953  IF (radius.LE.0.0) THEN
954  write(6,20) nid,xctr,yctr,zctr,iface
955  20 FORMAT(i5,'ERROR: Sphere of radius zero requested.'
956  $ ,/,5x,'EXITING in CRN3D',3e12.4,i3)
957  call exitt
958  ELSE
959  DO 40 i=1,4
960  radt=xcv(1,i,1)**2+xcv(2,i,1)**2+xcv(3,i,1)**2
961  radt=sqrt(radt)
962  test=abs(radt-radius)/radius
963  IF (test.GT.eps) THEN
964  write(6,30) nid
965  $ ,radt,radius,xcv(1,i,1),xcv(2,i,1),xcv(3,i,1)
966  30 FORMAT(i5,'ERROR: Element vertex not on requested sphere.'
967  $ ,/,5x,'EXITING in CRN3D',5e12.4)
968  call exitt
969  ENDIF
970  40 CONTINUE
971  ENDIF
972 C
973  return
974  end
975 c-----------------------------------------------------------------------
976  subroutine rotxyz
977 C-----------------------------------------------------------------------
978 C
979 C Establish rotation of undeformed elements.
980 C Used for fast evaluation of D*x and DT*x.
981 C Currently used for 2-d problems.
982 C
983 C-----------------------------------------------------------------------
984  include 'SIZE'
985  include 'INPUT'
986  include 'GEOM'
987  include 'ESOLV'
988  COMMON /fastmd/ ifdfrm(lelt), iffast(lelt), ifh2, ifsolv
989  LOGICAL IFDFRM, IFFAST, IFH2, IFSOLV
990 C
991  IF (ifmvbd) return
992  IF (ldim.EQ.3) return
993 C
994  eps = 1.e-6
995  epsinv = 1./eps
996  nxyz1 = lx1*ly1*lz1
997  DO 100 iel=1,nelv
998  IF (.NOT.ifdfrm(iel)) THEN
999  rxaver = vlsum(rxm1(1,1,1,iel),nxyz1)/(nxyz1)
1000  sxaver = vlsum(sxm1(1,1,1,iel),nxyz1)/(nxyz1)
1001  IF (sxaver.NE.0.) THEN
1002  rotxy = abs(sxaver/rxaver)
1003  IF (rotxy.LT.eps) THEN
1004  ifalgn(iel) = .true.
1005  ifrsxy(iel) = .true.
1006  ELSEIF (rotxy.GT.epsinv) THEN
1007  ifalgn(iel) = .true.
1008  ifrsxy(iel) = .false.
1009  ELSE
1010  ifalgn(iel) = .false.
1011  ifrsxy(iel) = .false.
1012  ENDIF
1013  ELSE
1014  ifalgn(iel) = .true.
1015  ifrsxy(iel) = .true.
1016  ENDIF
1017  ENDIF
1018  100 CONTINUE
1019  return
1020  end
1021 c-----------------------------------------------------------------------
1022  subroutine gensrf(XML,YML,ZML,IFCE,IE,MX,MY,MZ,zgml)
1023 C
1024 C 9 Mar 1994
1025 C
1026 C Program to generate surface deformations for NEKTON
1027 C input. Paul F. Fischer
1028 C
1029 c include 'basics.inc'
1030  include 'SIZE'
1031  include 'INPUT'
1032  include 'WZ'
1033  include 'TOPOL'
1034 C
1035  dimension xml(mx,my,mz,1),yml(mx,my,mz,1),zml(mx,my,mz,1)
1036  $ ,zgml(mx,3)
1037 C
1038  real IOPP(3),MXX(3),X0(3),DX(3)
1039 C
1040 C
1041 C Algorithm: .Project original point onto surface S
1042 C .Apply Gordon Hall to vector of points between x_s and
1043 C opposite face
1044 C
1045 C
1046  CALL dsset(mx,my,mz)
1047 C
1048  iface = eface1(ifce)
1049 c
1050 c Beware!! SKPDAT different from preprocessor/postprocessor!
1051 c
1052  js1 = skpdat(1,iface)
1053  jf1 = skpdat(2,iface)
1054  jskip1 = skpdat(3,iface)
1055  js2 = skpdat(4,iface)
1056  jf2 = skpdat(5,iface)
1057  jskip2 = skpdat(6,iface)
1058 c
1059  iopp(1) = mx-1
1060  iopp(2) = mx*(my-1)
1061  iopp(3) = mx*my*(mz-1)
1062  mxx(1) = mx
1063  mxx(2) = my
1064  mxx(3) = mz
1065  idir = 2*mod(iface,2) - 1
1066  ifc2 = (iface+1)/2
1067  i=0
1068 C
1069 C Find a characteristic length scale for initializing secant method
1070 C
1071  x0(1) = xml(js1,js2,1,ie)
1072  x0(2) = yml(js1,js2,1,ie)
1073  x0(3) = zml(js1,js2,1,ie)
1074  rmin = 1.0e16
1075 c
1076 c
1077 c
1078  DO 100 j2=js2,jf2,jskip2
1079  DO 100 j1=js1,jf1,jskip1
1080  if (j1.ne.js1.or.j2.ne.js2) then
1081  r2 = (x0(1) - xml(j1,j2,1,ie))**2
1082  $ + (x0(2) - yml(j1,j2,1,ie))**2
1083  $ + (x0(3) - zml(j1,j2,1,ie))**2
1084  rmin = min(r2,rmin)
1085  endif
1086  100 CONTINUE
1087  dxc = 0.05*sqrt(rmin)
1088 C
1089 C Project each point on this surface onto curved surface
1090 C
1091  DO 300 j2=js2,jf2,jskip2
1092  DO 300 j1=js1,jf1,jskip1
1093  i=i+1
1094  jopp = j1 + iopp(ifc2)*idir
1095  x0(1) = xml(j1,j2,1,ie)
1096  x0(2) = yml(j1,j2,1,ie)
1097  x0(3) = zml(j1,j2,1,ie)
1098 C
1099  call prjects(x0,dxc,curve(1,ifce,ie),ccurve(ifce,ie))
1100  dx(1) = x0(1)-xml(j1,j2,1,ie)
1101  dx(2) = x0(2)-yml(j1,j2,1,ie)
1102  dx(3) = x0(3)-zml(j1,j2,1,ie)
1103  mxs = mxx(ifc2)
1104  joff = (j1-jopp)/(mxs-1)
1105  DO 200 ix = 2,mxs
1106  j = jopp + joff*(ix-1)
1107  zeta = 0.5*(zgml(ix,1) + 1.0)
1108  xml(j,j2,1,ie) = xml(j,j2,1,ie)+dx(1)*zeta
1109  yml(j,j2,1,ie) = yml(j,j2,1,ie)+dx(2)*zeta
1110  zml(j,j2,1,ie) = zml(j,j2,1,ie)+dx(3)*zeta
1111  200 CONTINUE
1112  300 CONTINUE
1113 C
1114  return
1115  end
1116 c-----------------------------------------------------------------------
1117  subroutine prjects(x0,dxc,c,cc)
1118 c
1119 c Project the point x0 onto surface described by characteristics
1120 c given in the array c and cc.
1121 c
1122 c dxc - characteristic length scale used to estimate gradient.
1123 c
1124  real x0(3)
1125  real c(5)
1126  character*1 cc
1127  real x1(3)
1128  logical if3d
1129 c
1130  if3d = .true.
1131  if (dxc.le.0.0) then
1132  write(6,*) 'invalid dxc',dxc,x0
1133  write(6,*) 'Abandoning prjects'
1134  return
1135  endif
1136 c
1137  call copy(x1,x0,3)
1138  r0 = ressrf(x0,c,cc)
1139  if (r0.eq.0) return
1140 c
1141 c Must at least use ctr differencing to capture symmetry!
1142 c
1143  x1(1) = x0(1) - dxc
1144  r1 = ressrf(x1,c,cc)
1145  x1(1) = x0(1) + dxc
1146  r2 = ressrf(x1,c,cc)
1147  x1(1) = x0(1)
1148  rx = 0.5*(r2-r1)/dxc
1149 c
1150  x1(2) = x0(2) - dxc
1151  r1 = ressrf(x1,c,cc)/dxc
1152  x1(2) = x0(2) + dxc
1153  r2 = ressrf(x1,c,cc)/dxc
1154  x1(2) = x0(2)
1155  ry = 0.5*(r2-r1)/dxc
1156 c
1157  if (if3d) then
1158  x1(3) = x0(3) - dxc
1159  r1 = ressrf(x1,c,cc)/dxc
1160  x1(3) = x0(3) + dxc
1161  r2 = ressrf(x1,c,cc)/dxc
1162  rz = 0.5*(r2-r1)/dxc
1163  endif
1164  rnorm2 = rx**2 + ry**2 + rz**2
1165  alpha = - r0/rnorm2
1166 c
1167 c Apply secant method: Use an initial segment twice expected length
1168 c
1169  x1(1) = x0(1) + 2.0*rx * alpha
1170  x1(2) = x0(2) + 2.0*ry * alpha
1171  x1(3) = x0(3) + 2.0*rz * alpha
1172  call srfind(x1,x0,c,cc)
1173 c
1174 c write(6,6) cc,c(2),c(3),x0,x1
1175 c 6 format(1x,a1,1x,2f5.2,3f9.4,3x,3f9.4)
1176 c
1177  call copy(x0,x1,3)
1178 c
1179  return
1180  end
1181 c-----------------------------------------------------------------------
1182  subroutine srfind(x1,x0,c,cc)
1183  real x1(3),x0(3)
1184  real c(5)
1185  character*1 cc
1186 c
1187 c Find point on line segment that intersects the ellipsoid
1188 c specified by:
1189 c (x/a)**2 + (y/b)**2 + (z/b)**2 = 1
1190 c
1191 c
1192 c Algorithm: 4 rounds of secant x_k+1 = x_k - f/f'
1193 c
1194  a0 = 0.0
1195  a1 = 1.0
1196  r0 = ressrf(x0,c,cc)
1197  dx = x1(1) - x0(1)
1198  dy = x1(2) - x0(2)
1199  dz = x1(3) - x0(3)
1200 c write(6,*) 'dxyz',dx,dy,dz
1201 c write(6,*) 'cc ',x0,cc,c(2),c(3)
1202  do 10 i=1,9
1203  r1 = ressrf(x1,c,cc)
1204  if (r1.ne.r0) then
1205  da = r1*(a1-a0)/(r1-r0)
1206  r0 = r1
1207  a0 = a1
1208  a1 = a1 - da
1209  endif
1210  x1(1) = x0(1) + a1*dx
1211  x1(2) = x0(2) + a1*dy
1212  x1(3) = x0(3) + a1*dz
1213  10 continue
1214 c write(6,*) ' r1',r1,r0,a1
1215  return
1216  end
1217 c-----------------------------------------------------------------------
1218  function ressrf(x,c,cc)
1219  real x(3)
1220  real c(5)
1221  character*1 cc
1222 c
1223  ressrf = 0.0
1224  if (cc.eq.'e') then
1225  a = c(2)
1226  b = c(3)
1227  ressrf = 1.0 - (x(1)/a)**2 - (x(2)/b)**2 - (x(3)/b)**2
1228  return
1229  endif
1230 c
1231  return
1232  end
1233 c-----------------------------------------------------------------------
1234  subroutine linquad(xl,yl,zl,nxl,nyl,nzl)
1235 
1236  include 'SIZE'
1237  include 'WZ'
1238  include 'GEOM'
1239  include 'TOPOL'
1240  include 'INPUT'
1241  include 'PARALLEL'
1242 
1243  real xl(nxl*nyl*nzl,1),yl(nxl*nyl*nzl,1),zl(nxl*nyl*nzl,1)
1244 
1245  integer e
1246  logical ifmid
1247 
1248  nedge = 4 + 8*(ldim-2)
1249 
1250  do e=1,nelt ! Loop over all elements
1251 
1252  ifmid = .false.
1253  do k=1,nedge
1254  if (ccurve(k,e).eq.'m') ifmid = .true.
1255  enddo
1256 
1257  if (lx1.eq.2) ifmid = .false.
1258  if (ifmid) then
1259  call xyzquad(xl(1,e),yl(1,e),zl(1,e),nxl,nyl,nzl,e)
1260  else
1261  call xyzlin (xl(1,e),yl(1,e),zl(1,e),nxl,nyl,nzl,e,ifaxis)
1262  endif
1263  enddo
1264 
1265  return
1266  end
1267 c-----------------------------------------------------------------------
1268  subroutine xyzlin(xl,yl,zl,nxl,nyl,nzl,e,ifaxl)
1269 c Generate bi- or trilinear mesh
1270 
1271  include 'SIZE'
1272  include 'INPUT'
1273 
1274  real xl(nxl,nyl,nzl),yl(nxl,nyl,nzl),zl(nxl,nyl,nzl)
1275  integer e
1276  logical ifaxl ! local ifaxis specification
1277 
1278 c Preprocessor Corner notation: Symmetric Corner notation:
1279 c
1280 c 4+-----+3 ^ s 3+-----+4 ^ s
1281 c / /| | / /| |
1282 c / / | | / / | |
1283 c 8+-----+7 +2 +----> r 7+-----+8 +2 +----> r
1284 c | | / / | | / /
1285 c | |/ / | |/ /
1286 c 5+-----+6 t 5+-----+6 t
1287 
1288  integer indx(8)
1289  save indx
1290  data indx / 1,2,4,3,5,6,8,7 /
1291 
1292  parameter(ldw=4*lx1*ly1*lz1)
1293  common /ctmp0/ xcb(2,2,2),ycb(2,2,2),zcb(2,2,2),w(ldw)
1294 
1295  common /cxyzl/ zgml(lx1,3),jx(lx1*2),jy(lx1*2),jz(lx1*2)
1296  $ ,jxt(lx1*2),jyt(lx1*2),jzt(lx1*2)
1297  $ ,zlin(2)
1298  real jx,jy,jz,jxt,jyt,jzt
1299 
1300  call setzgml (zgml,e,nxl,nyl,nzl,ifaxl)
1301 
1302  zlin(1) = -1
1303  zlin(2) = 1
1304 
1305  k = 1
1306  do i=1,nxl
1307  call fd_weights_full(zgml(i,1),zlin,1,0,jxt(k))
1308  call fd_weights_full(zgml(i,2),zlin,1,0,jyt(k))
1309  call fd_weights_full(zgml(i,3),zlin,1,0,jzt(k))
1310  k=k+2
1311  enddo
1312  call transpose(jx,nxl,jxt,2)
1313 
1314  ldim2 = 2**ldim
1315  do ix=1,ldim2 ! Convert prex notation to lexicographical
1316  i=indx(ix)
1317  xcb(ix,1,1)=xc(i,e)
1318  ycb(ix,1,1)=yc(i,e)
1319  zcb(ix,1,1)=zc(i,e)
1320  enddo
1321 
1322 c Map R-S-T space into physical X-Y-Z space.
1323 
1324  ! NOTE: Assumes nxl=nyl=nzl !
1325 
1326  call tensr3(xl,nxl,xcb,2,jx,jyt,jzt,w)
1327  call tensr3(yl,nxl,ycb,2,jx,jyt,jzt,w)
1328  call tensr3(zl,nxl,zcb,2,jx,jyt,jzt,w)
1329 
1330  return
1331  end
1332 c-----------------------------------------------------------------------
1333  subroutine xyzquad(xl,yl,zl,nxl,nyl,nzl,e)
1334 c Generate bi- or trilinear mesh
1335 
1336  include 'SIZE'
1337  include 'INPUT'
1338 
1339  real xl(nxl,nyl,nzl),yl(nxl,nyl,nzl),zl(nxl,nyl,nzl)
1340  real xq(27),yq(27),zq(27)
1341  integer e
1342 
1343  parameter(ldw=4*lx1*ly1*lz1)
1344  common /ctmp0/ w(ldw,2),zg(3)
1345 
1346 c Note : CTMP1 is used in this format in several subsequent routines
1347 
1348  integer eindx(12) ! index of 12 edges into 3x3x3 tensor
1349  save eindx ! Follows preprocessor notation..
1350  data eindx / 2 , 6 , 8 , 4
1351  $ , 20 , 24 , 26 , 22
1352  $ , 10 , 12 , 18 , 16 / ! preproc. vtx notation
1353 
1354  common /cxyzl/ zgml(lx1,3),jx(lx1*3),jy(lx1*3),jz(lx1*3)
1355  $ ,jxt(lx1*3),jyt(lx1*3),jzt(lx1*3)
1356  $ ,zquad(3)
1357  real jx,jy,jz,jxt,jyt,jzt
1358 
1359  call xyzlin(xq,yq,zq,3,3,3,e,.false.) ! map bilin to 3x3x3
1360 
1361  nedge = 4 + 8*(ldim-2)
1362 
1363  do k=1,nedge
1364  if (ccurve(k,e).eq.'m') then
1365  j = eindx(k)
1366  xq(j) = curve(1,k,e)
1367  yq(j) = curve(2,k,e)
1368  zq(j) = curve(3,k,e)
1369  endif
1370  enddo
1371 
1372  zg(1) = -1
1373  zg(2) = 0
1374  zg(3) = 1
1375 
1376  if (if3d) then
1377  call gh_face_extend(xq,zg,3,2,w(1,1),w(1,2)) ! 2 --> edge extend
1378  call gh_face_extend(yq,zg,3,2,w(1,1),w(1,2))
1379  call gh_face_extend(zq,zg,3,2,w(1,1),w(1,2))
1380  else
1381  call gh_face_extend_2d(xq,zg,3,2,w(1,1),w(1,2)) ! 2 --> edge extend
1382  call gh_face_extend_2d(yq,zg,3,2,w(1,1),w(1,2))
1383  endif
1384  call clean_xyzq(xq,yq,zq,if3d) ! verify that midside node is in "middle"
1385 
1386 
1387 c Map R-S-T space into physical X-Y-Z space.
1388  ! NOTE: Assumes nxl=nyl=nzl !
1389 
1390  zquad(1) = -1
1391  zquad(2) = 0
1392  zquad(3) = 1
1393 
1394  call setzgml (zgml,e,nxl,nyl,nzl,ifaxis) ! Here we address axisymm.
1395 
1396  k = 1
1397  do i=1,nxl
1398  call fd_weights_full(zgml(i,1),zquad,2,0,jxt(k))
1399  call fd_weights_full(zgml(i,2),zquad,2,0,jyt(k))
1400  call fd_weights_full(zgml(i,3),zquad,2,0,jzt(k))
1401  k=k+3
1402  enddo
1403  call transpose(jx,nxl,jxt,3)
1404 
1405  call tensr3(xl,nxl,xq,3,jx,jyt,jzt,w)
1406  call tensr3(yl,nxl,yq,3,jx,jyt,jzt,w)
1407  call tensr3(zl,nxl,zq,3,jx,jyt,jzt,w)
1408 
1409  return
1410  end
1411 c-----------------------------------------------------------------------
1412  subroutine clean_xyzq(x,y,z,if3d) ! verify that midside node is in "middle"
1413 
1414  real x(3,3,3),y(3,3,3),z(3,3,3)
1415  logical if3d
1416 
1417  integer eindx(12) ! index of 12 edges into 3x3x3 tensor
1418  save eindx ! Follows preprocessor notation..
1419  data eindx / 2 , 6 , 8 , 4
1420  $ , 20 , 24 , 26 , 22
1421  $ , 10 , 12 , 18 , 16 / ! preproc. vtx notation
1422 
1423 
1424 c if (if3d) then ! 12 edges
1425 c else
1426 c endif
1427 
1428 c Here - see routine "fix_m_curve" in prenek for indexing strategy
1429 c
1430 c Note that "fix_m_curve" does not yet perform the actual fix, and
1431 c it too should be updated.
1432 
1433  return
1434  end
1435 c-----------------------------------------------------------------------
1436  subroutine mesh_metrics
1437  include 'SIZE'
1438  include 'TOTAL'
1439 
1440  parameter(nedge = 4 + 8*(ldim-2))
1441  real ledg(nedge)
1442 
1443  nxyz = nx1*ny1*nz1
1444  ntot = nxyz*nelt
1445 
1446  ! aspect ratio
1447  ddmin = 1e20
1448  ddmax = -1
1449  ddavg = 0
1450  do ie = 1,nelt
1451  ledg(1) = dist_xyzc(1,2,ie)
1452  ledg(2) = dist_xyzc(1,4,ie)
1453  ledg(3) = dist_xyzc(2,3,ie)
1454  ledg(4) = dist_xyzc(4,3,ie)
1455  if (ndim.eq.3) then
1456  ledg(5) = dist_xyzc(1,5,ie)
1457  ledg(6) = dist_xyzc(2,6,ie)
1458  ledg(7) = dist_xyzc(4,8,ie)
1459  ledg(8) = dist_xyzc(3,7,ie)
1460 
1461  ledg(9) = dist_xyzc(5,6,ie)
1462  ledg(10) = dist_xyzc(5,8,ie)
1463  ledg(11) = dist_xyzc(8,7,ie)
1464  ledg(12) = dist_xyzc(6,7,ie)
1465  endif
1466 
1467  dratio = vlmax(ledg,nedge)/vlmin(ledg,nedge)
1468  ddmin = min(ddmin,dratio)
1469  ddmax = max(ddmax,dratio)
1470  ddavg = ddavg + dratio
1471  enddo
1472  darmin = glmin(ddmin,1)
1473  darmax = glmax(ddmax,1)
1474  daravg = glsum(ddavg,1)/nelgt
1475 
1476  ! scaled Jac
1477  ddmin = 1e20
1478  ddmax = -1
1479  ddavg = 0
1480  do ie = 1,nelt
1481  dratio = vlmin(jacm1(1,1,1,ie),nxyz)/
1482  $ vlmax(jacm1(1,1,1,ie),nxyz)
1483  ddmin = min(ddmin,dratio)
1484  ddmax = max(ddmax,dratio)
1485  ddavg = ddavg + dratio
1486  enddo
1487  dsjmin = glmin(ddmin,1)
1488  dsjmax = glmax(ddmax,1)
1489  dsjavg = glsum(ddavg,1)/nelgt
1490 
1491  ! dx
1492  ddmin = 1e20
1493  ddmax = -1
1494  ddavg = 0
1495  do ie = 1,nelt
1496  ddmin = min(ddmin,dxmin_e(ie))
1497  ddmax = max(ddmax,dxmax_e(ie))
1498  dtmp = vlsum(jacm1(1,1,1,ie),nxyz)**1./3
1499  ddavg = ddavg + dtmp/(lx1-1)
1500  enddo
1501  dxmin = glmin(ddmin,1)
1502  dxmax = glmax(ddmax,1)
1503  dxavg = glsum(ddavg,1)/nelgt
1504 
1505  if (nid.eq.0) then
1506  write(6,*) 'mesh metrics:'
1507  write(6,'(A,1p2E9.2)') ' GLL grid spacing min/max :',
1508  $ dxmin,dxmax
1509  write(6,'(A,1p3E9.2)') ' scaled Jacobian min/max/avg:',
1510  $ dsjmin,dsjmax,dsjavg
1511  write(6,'(A,1p3E9.2)') ' aspect ratio min/max/avg:',
1512  $ darmin,darmax,daravg
1513  write(6,*)
1514  endif
1515 
1516  return
1517  end
1518 c-----------------------------------------------------------------------
1519  real function dist_xyzc(i,j,ie)
1520 c
1521 c distance between two element corner points
1522 c
1523  include 'SIZE'
1524  include 'INPUT'
1525 
1526  dist_xyzc = (xc(i,ie) - xc(j,ie))**2
1527  dist_xyzc = dist_xyzc + (yc(i,ie) - yc(j,ie))**2
1528  if(ndim.eq.3) dist_xyzc = dist_xyzc + (zc(i,ie) - zc(j,ie))**2
1529  dist_xyzc = sqrt(dist_xyzc)
1530 
1531  return
1532  end
1533 c-----------------------------------------------------------------------
1534  function dxmin_e(e)
1535 
1536  include 'SIZE'
1537  include 'GEOM'
1538  include 'INPUT'
1539 
1540  integer e
1541 
1542  d2m = 1.e20
1543 
1544  if (ldim.eq.3) then
1545  do k=1,nz1-1
1546  do j=1,ny1-1
1547  do i=1,nx1-1
1548  dx = xm1(i+1,j,k,e) - xm1(i,j,k,e)
1549  dy = ym1(i+1,j,k,e) - ym1(i,j,k,e)
1550  dz = zm1(i+1,j,k,e) - zm1(i,j,k,e)
1551  d2 = dx*dx + dy*dy + dz*dz
1552  d2m = min(d2m,d2)
1553 
1554  dx = xm1(i,j+1,k,e) - xm1(i,j,k,e)
1555  dy = ym1(i,j+1,k,e) - ym1(i,j,k,e)
1556  dz = zm1(i,j+1,k,e) - zm1(i,j,k,e)
1557  d2 = dx*dx + dy*dy + dz*dz
1558  d2m = min(d2m,d2)
1559 
1560  dx = xm1(i,j,k+1,e) - xm1(i,j,k,e)
1561  dy = ym1(i,j,k+1,e) - ym1(i,j,k,e)
1562  dz = zm1(i,j,k+1,e) - zm1(i,j,k,e)
1563  d2 = dx*dx + dy*dy + dz*dz
1564  d2m = min(d2m,d2)
1565  enddo
1566  enddo
1567  enddo
1568  else ! 2D
1569  do j=1,ny1-1
1570  do i=1,nx1-1
1571  dx = xm1(i+1,j,1,e) - xm1(i,j,1,e)
1572  dy = ym1(i+1,j,1,e) - ym1(i,j,1,e)
1573  d2 = dx*dx + dy*dy
1574  d2m = min(d2m,d2)
1575 
1576  dx = xm1(i,j+1,1,e) - xm1(i,j,1,e)
1577  dy = ym1(i,j+1,1,e) - ym1(i,j,1,e)
1578  d2 = dx*dx + dy*dy
1579  d2m = min(d2m,d2)
1580  enddo
1581  enddo
1582  endif
1583 
1584  dxmin_e = sqrt(d2m)
1585 
1586  return
1587  end
1588 c-----------------------------------------------------------------------
1589  function dxmax_e(e)
1590 
1591  include 'SIZE'
1592  include 'GEOM'
1593  include 'INPUT'
1594 
1595  integer e
1596 
1597  d2m = -1.e20
1598 
1599  if (ldim.eq.3) then
1600  do k=1,nz1-1
1601  do j=1,ny1-1
1602  do i=1,nx1-1
1603  dx = xm1(i+1,j,k,e) - xm1(i,j,k,e)
1604  dy = ym1(i+1,j,k,e) - ym1(i,j,k,e)
1605  dz = zm1(i+1,j,k,e) - zm1(i,j,k,e)
1606  d2 = dx*dx + dy*dy + dz*dz
1607  d2m = max(d2m,d2)
1608 
1609  dx = xm1(i,j+1,k,e) - xm1(i,j,k,e)
1610  dy = ym1(i,j+1,k,e) - ym1(i,j,k,e)
1611  dz = zm1(i,j+1,k,e) - zm1(i,j,k,e)
1612  d2 = dx*dx + dy*dy + dz*dz
1613  d2m = max(d2m,d2)
1614 
1615  dx = xm1(i,j,k+1,e) - xm1(i,j,k,e)
1616  dy = ym1(i,j,k+1,e) - ym1(i,j,k,e)
1617  dz = zm1(i,j,k+1,e) - zm1(i,j,k,e)
1618  d2 = dx*dx + dy*dy + dz*dz
1619  d2m = max(d2m,d2)
1620  enddo
1621  enddo
1622  enddo
1623  else ! 2D
1624  do j=1,ny1-1
1625  do i=1,nx1-1
1626  dx = xm1(i+1,j,1,e) - xm1(i,j,1,e)
1627  dy = ym1(i+1,j,1,e) - ym1(i,j,1,e)
1628  d2 = dx*dx + dy*dy
1629  d2m = max(d2m,d2)
1630 
1631  dx = xm1(i,j+1,1,e) - xm1(i,j,1,e)
1632  dy = ym1(i,j+1,1,e) - ym1(i,j,1,e)
1633  d2 = dx*dx + dy*dy
1634  d2m = max(d2m,d2)
1635  enddo
1636  enddo
1637  endif
1638 
1639  dxmax_e = sqrt(d2m)
1640 
1641  return
1642  end
1643 c-----------------------------------------------------------------------
exitt
void exitt()
Definition: comm_mpi.f:604
exitti
subroutine exitti(stringi, idata)
Definition: comm_mpi.f:535
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
fd_weights_full
subroutine fd_weights_full(xx, x, n, m, c)
Definition: fast3d.f:1293
tensr3
subroutine tensr3(v, nv, u, nu, A, Bt, Ct, w)
Definition: fasts.f:126
dxmin_e
function dxmin_e(e)
Definition: genxyz.f:1535
dot
real function dot(V1, V2, N)
Definition: genxyz.f:885
dist_xyzc
real function dist_xyzc(i, j, ie)
Definition: genxyz.f:1520
genxyz
subroutine genxyz(xml, yml, zml, nxl, nyl, nzl)
Definition: genxyz.f:558
xyzlin
subroutine xyzlin(xl, yl, zl, nxl, nyl, nzl, e, ifaxl)
Definition: genxyz.f:1269
ressrf
function ressrf(x, c, cc)
Definition: genxyz.f:1219
defsrf
subroutine defsrf(xml, yml, zml, nxl, nyl, nzl, ie, iface1, ccv)
Definition: genxyz.f:104
setdef
subroutine setdef
Definition: genxyz.f:352
zbrent
function zbrent(X1, X2, TOL)
Definition: genxyz.f:238
arcsrf
subroutine arcsrf(xml, yml, zml, nxl, nyl, nzl, ie, isid)
Definition: genxyz.f:3
crn3d
subroutine crn3d(xcv, xc, yc, zc, curve, iface)
Definition: genxyz.f:929
intrsc
subroutine intrsc(x3, x2, x1, delt, ie, iface)
Definition: genxyz.f:168
gensrf
subroutine gensrf(XML, YML, ZML, IFCE, IE, MX, MY, MZ, zgml)
Definition: genxyz.f:1023
gencoor
subroutine gencoor(xm3, ym3, zm3)
Definition: genxyz.f:534
linquad
subroutine linquad(xl, yl, zl, nxl, nyl, nzl)
Definition: genxyz.f:1235
prjects
subroutine prjects(x0, dxc, c, cc)
Definition: genxyz.f:1118
edg3d
subroutine edg3d(xysrf, x1, x2, i1, i2, j1, j2, nx, ny)
Definition: genxyz.f:844
sphsrf
subroutine sphsrf(xml, yml, zml, ifce, ie, nx, ny, nz, xysrf)
Definition: genxyz.f:677
xyzquad
subroutine xyzquad(xl, yl, zl, nxl, nyl, nzl, e)
Definition: genxyz.f:1334
setzgml
subroutine setzgml(zgml, e, nxl, nyl, nzl, ifaxl)
Definition: genxyz.f:641
ifvchk
logical function ifvchk(VEC, I1, I2, I3)
Definition: genxyz.f:510
srfind
subroutine srfind(x1, x0, c, cc)
Definition: genxyz.f:1183
rotxyz
subroutine rotxyz
Definition: genxyz.f:977
zbrac
subroutine zbrac(x1, x2, succes)
Definition: genxyz.f:207
norm3d
subroutine norm3d(v1)
Definition: genxyz.f:912
cross
subroutine cross(v1, v2, v3)
Definition: genxyz.f:899
dxmax_e
function dxmax_e(e)
Definition: genxyz.f:1590
clean_xyzq
subroutine clean_xyzq(x, y, z, if3d)
Definition: genxyz.f:1413
fnc
function fnc(ETA)
Definition: genxyz.f:325
sethmat
subroutine sethmat(h, zgml, nxl, nyl, nzl)
Definition: genxyz.f:613
mesh_metrics
subroutine mesh_metrics
Definition: genxyz.f:1437
glmin
function glmin(a, n)
Definition: math.f:973
transpose
subroutine transpose(a, lda, b, ldb)
Definition: math.f:1594
mod1
function mod1(i, n)
Definition: math.f:509
vlmax
real function vlmax(vec, n)
Definition: math.f:396
copy
subroutine copy(a, b, n)
Definition: math.f:260
addtnsr
subroutine addtnsr(s, h1, h2, h3, nx, ny, nz)
Definition: math.f:477
vlsum
real function vlsum(vec, n)
Definition: math.f:417
glsum
function glsum(x, n)
Definition: math.f:861
cmult
subroutine cmult(a, const, n)
Definition: math.f:315
rone
subroutine rone(a, n)
Definition: math.f:230
vlmin
real function vlmin(vec, n)
Definition: math.f:357
rzero
subroutine rzero(a, n)
Definition: math.f:208
glmax
function glmax(a, n)
Definition: math.f:960
gh_face_extend_2d
subroutine gh_face_extend_2d(x, zg, n, gh_type, e, v)
Definition: navier5.f:2437
gh_face_extend
subroutine gh_face_extend(x, zg, n, gh_type, e, v)
Definition: navier5.f:2422