KTH framework for Nek5000 toolboxes; testing version  0.0.1
coef.f
Go to the documentation of this file.
1  subroutine genwz
2 C-----------------------------------------------------------------
3 C
4 C GENERATE
5 C
6 C - DERIVATIVE OPERATORS
7 C - INTERPOLATION OPERATORS
8 C - WEIGHTS
9 C - COLLOCATION POINTS
10 C
11 C ASSOCIATED WITH THE
12 C
13 C - GAUSS-LOBATTO LEGENDRE MESH (SUFFIX M1/M2/M3)
14 C - GAUSS LEGENDRE MESH (SUFFIX M2)
15 C - GAUSS-LOBATTO JACOBI MESH (SUFFIX M1/M2/M3)
16 C
17 C-----------------------------------------------------------------
18 C
19  include 'SIZE'
20  include 'WZ'
21  include 'DXYZ'
22  include 'IXYZ'
23  include 'INPUT'
24 
25  REAL TMP(LY1,LY1),TMPT(LY1,LY1)
26 C
27  IF (ldim.EQ.2) THEN
28 C
29 C*** Two-dimensional case **********************
30 C
31 C
32 C Gauss-Lobatto Legendre mesh (suffix M1)
33 C Generate collocation points and weights
34 C
35  CALL zwgll (zgm1(1,1),wxm1,lx1)
36  CALL zwgll (zgm1(1,2),wym1,ly1)
37  zgm1(lz1,3) = 0.
38  wzm1(lz1) = 1.
39  DO 100 iy=1,ly1
40  DO 100 ix=1,lx1
41  w3m1(ix,iy,1)=wxm1(ix)*wym1(iy)
42  100 CONTINUE
43 C
44 C Compute derivative matrices
45 C
46  CALL dgll (dxm1,dxtm1,zgm1(1,1),lx1,lx1)
47  CALL dgll (dym1,dytm1,zgm1(1,2),ly1,ly1)
48  CALL rzero (dzm1 ,lz1*lz1)
49  CALL rzero (dztm1,lz1*lz1)
50 C
51 C Gauss Legendre mesh (suffix M2)
52 C Generate collocation points and weights
53 C
54  IF(ifsplit)THEN
55  CALL zwgll (zgm2(1,1),wxm2,lx2)
56  CALL zwgll (zgm2(1,2),wym2,ly2)
57  ELSE
58  CALL zwgl (zgm2(1,1),wxm2,lx2)
59  CALL zwgl (zgm2(1,2),wym2,ly2)
60  ENDIF
61  zgm2(lz2,3) = 0.
62  wzm2(lz2) = 1.
63  DO 200 iy=1,ly2
64  DO 200 ix=1,lx2
65  w3m2(ix,iy,1)=wxm2(ix)*wym2(iy)
66  200 CONTINUE
67 C
68 C Gauss-Lobatto Legendre mesh (suffix M3).
69 C Generate collocation points and weights.
70 C
71  CALL zwgll (zgm3(1,1),wxm3,lx3)
72  CALL zwgll (zgm3(1,2),wym3,ly3)
73  zgm3(lz3,3) = 0.
74  wzm3(lz3) = 1.
75  DO 300 iy=1,ly3
76  DO 300 ix=1,lx3
77  w3m3(ix,iy,1)=wxm3(ix)*wym3(iy)
78  300 CONTINUE
79 C
80 C Compute derivative matrices
81 C
82  CALL dgll (dxm3,dxtm3,zgm3(1,1),lx3,lx3)
83  CALL dgll (dym3,dytm3,zgm3(1,2),ly3,ly3)
84  CALL rzero (dzm3 ,lz3*lz3)
85  CALL rzero (dztm3,lz3*lz3)
86 C
87 C Generate interpolation operators for the staggered mesh
88 C
89  CALL igllm (ixm12,ixtm12,zgm1(1,1),zgm2(1,1),lx1,lx2,lx1,lx2)
90  CALL igllm (iym12,iytm12,zgm1(1,2),zgm2(1,2),ly1,ly2,ly1,ly2)
91  izm12(lz2,lz1) = 1.
92  iztm12(lz1,lz2) = 1.
93 C
94 C NOTE: The splitting scheme has only one mesh!!!!!
95 C
96  IF (ifsplit) THEN
97  CALL igllm (ixm21,ixtm21,zgm1(1,1),zgm2(1,1),lx1,lx2,lx1,lx2)
98  CALL igllm (iym21,iytm21,zgm1(1,2),zgm2(1,2),ly1,ly2,ly1,ly2)
99  ELSE
100  CALL iglm (ixm21,ixtm21,zgm2(1,1),zgm1(1,1),lx2,lx1,lx2,lx1)
101  CALL iglm (iym21,iytm21,zgm2(1,2),zgm1(1,2),ly2,ly1,ly2,ly1)
102  ENDIF
103  izm21(lz1,lz2) = 1.
104  iztm21(lz2,lz1) = 1.
105 C
106 C Compute derivative operators for the staggered mesh
107 C
108  IF(ifsplit)THEN
109  CALL copy (dxm12, dxm1, lx1*lx2)
110  CALL copy (dxtm12,dxtm1,lx1*lx2)
111  CALL copy (dym12, dym1, ly1*ly2)
112  CALL copy (dytm12,dytm1,ly1*ly2)
113  CALL copy (dzm12, dzm1, lz1*lz2)
114  CALL copy (dztm12,dztm1,lz1*lz2)
115  ELSE
116  CALL dgllgl (dxm12,dxtm12,zgm1(1,1),zgm2(1,1),ixm12,
117  $ lx1,lx2,lx1,lx2)
118  CALL dgllgl (dym12,dytm12,zgm1(1,2),zgm2(1,2),iym12,
119  $ ly1,ly2,ly1,ly2)
120  dzm12(lz2,lz1) = 0.
121  dztm12(lz2,lz1) = 0.
122  ENDIF
123 C
124 C Compute interpolation operators for the geometry mesh M3.
125 C
126  CALL igllm (ixm13,ixtm13,zgm1(1,1),zgm3(1,1),lx1,lx3,lx1,lx3)
127  CALL igllm (iym13,iytm13,zgm1(1,2),zgm3(1,2),ly1,ly3,ly1,ly3)
128  CALL igllm (ixm31,ixtm31,zgm3(1,1),zgm1(1,1),lx3,lx1,lx3,lx1)
129  CALL igllm (iym31,iytm31,zgm3(1,2),zgm1(1,2),ly3,ly1,ly3,ly1)
130  izm13(lz3,lz1) = 1.
131  iztm13(lz1,lz3) = 1.
132  izm31(lz1,lz3) = 1.
133  iztm31(lz3,lz1) = 1.
134 C
135 C
136  IF (ifaxis) THEN
137 C
138 C Special treatment for the axisymmetric case
139 C Generate additional points, weights, derivative operators and
140 C interpolation operators required for elements close to the axis.
141 C
142 C
143 C Gauss-Lobatto Jacobi mesh (suffix M1).
144 C Generate collocation points and weights (alpha=0, beta=1).
145 C
146  alpha = 0.
147  beta = 1.
148  CALL zwglj (zam1,wam1,ly1,alpha,beta)
149  DO 400 iy=1,ly1
150  DO 400 ix=1,lx1
151  w2am1(ix,iy)=wxm1(ix)*wam1(iy)
152  w2cm1(ix,iy)=wxm1(ix)*wym1(iy)
153  400 CONTINUE
154 C
155 C Compute derivative matrices
156 C
157  CALL copy (dcm1,dym1,ly1*ly1)
158  CALL copy (dctm1,dytm1,ly1*ly1)
159  CALL dglj (dam1,datm1,zam1,ly1,ly1,alpha,beta)
160 C
161 C Gauss Jacobi mesh (suffix M2)
162 C Generate collocation points and weights
163 C
164  IF(ifsplit)THEN
165  CALL zwglj (zam2,wam2,ly2,alpha,beta)
166  ELSE
167  CALL zwgj (zam2,wam2,ly2,alpha,beta)
168  ENDIF
169  DO 500 iy=1,ly2
170  DO 500 ix=1,lx2
171  w2cm2(ix,iy)=wxm2(ix)*wym2(iy)
172  w2am2(ix,iy)=wxm2(ix)*wam2(iy)
173  500 CONTINUE
174 C
175 C Gauss-Lobatto Jacobi mesh (suffix M3).
176 C Generate collocation points and weights.
177 C
178  CALL zwglj (zam3,wam3,ly3,alpha,beta)
179  DO 600 iy=1,ly3
180  DO 600 ix=1,lx3
181  w2cm3(ix,iy)=wxm3(ix)*wym3(iy)
182  w2am3(ix,iy)=wxm3(ix)*wam3(iy)
183  600 CONTINUE
184 C
185 C Compute derivative matrices
186 C
187  CALL copy (dcm3,dym3,ly3*ly3)
188  CALL copy (dctm3,dytm3,ly3*ly3)
189  CALL dglj (dam3,datm3,zam3,ly3,ly3,alpha,beta)
190 C
191 C Generate interpolation operators for the staggered mesh
192 C
193  CALL copy (icm12,iym12,ly2*ly1)
194  CALL copy (ictm12,iytm12,ly1*ly2)
195  CALL igljm (iam12,iatm12,zam1,zam2,ly1,ly2,ly1,ly2,alpha,beta)
196  CALL copy (icm21,iym21,ly1*ly2)
197  CALL copy (ictm21,iytm21,ly2*ly1)
198  IF (ifsplit) THEN
199  CALL igljm (iam21,iatm21,zam2,zam1,ly1,ly2,ly1,ly2,alpha,beta)
200  ELSE
201  CALL igjm (iam21,iatm21,zam2,zam1,ly2,ly1,ly2,ly1,alpha,beta)
202  ENDIF
203 C
204 C Compute derivative operators for the staggered mesh
205 C
206  CALL copy (dcm12,dym12,ly2*ly1)
207  CALL copy (dctm12,dytm12,ly1*ly2)
208  IF(ifsplit)THEN
209  CALL copy (dam12, dam1, ly1*ly2)
210  CALL copy (datm12,datm1,ly1*ly2)
211  ELSE
212  CALL dgljgj (dam12,datm12,zam1,zam2,iam12,
213  $ ly1,ly2,ly1,ly2,alpha,beta)
214  ENDIF
215 C
216 C Compute interpolation operators for the geometry mesh M3.
217 C
218  CALL copy (icm13,iym13,ly3*ly1)
219  CALL copy (ictm13,iytm13,ly1*ly3)
220  CALL igljm (iam13,iatm13,zam1,zam3,ly1,ly3,ly1,ly3,alpha,beta)
221  CALL copy (icm31,iym31,ly1*ly3)
222  CALL copy (ictm31,iytm31,ly3*ly1)
223  CALL igljm (iam31,iatm31,zam3,zam1,ly3,ly1,ly3,ly1,alpha,beta)
224 C
225 C Compute interpolation operators between Gauss-Lobatto Jacobi
226 C and Gauss-Lobatto Legendre (to be used in PREPOST).
227 C
228  CALL igljm(iajl1,iatjl1,zam1,zgm1(1,2),ly1,ly1,ly1,ly1,alpha,beta)
229  IF (ifsplit) THEN
230  CALL igljm(iajl2,iatjl2,zam2,zgm2(1,2),ly2,ly2,ly2,ly2,alpha,beta)
231  ELSE
232  CALL igjm (iajl2,iatjl2,zam2,zgm2(1,2),ly2,ly2,ly2,ly2,alpha,beta)
233  ENDIF
234 
235  CALL invmt(iajl1 ,ialj1 ,tmp ,ly1)
236  CALL invmt(iatjl1,iatlj1,tmpt,ly1)
237  CALL mxm (iatjl1,ly1,iatlj1,ly1,tmpt,ly1)
238  CALL mxm (iajl1 ,ly1,ialj1 ,ly1,tmp ,ly1)
239 
240 C
241 C Compute interpolation operators between Gauss-Lobatto Legendre
242 C and Gauss-Lobatto Jacobi (to be used in subr. genxyz IN postpre).
243 C
244 c
245 c This call is not right, and these arrays are not used. 3/27/02. pff
246 c CALL IGLLM(IALJ3,IATLJ3,ZGM3(1,2),ZAM3,ly3,ly3,ly3,ly3,ALPHA,BETA)
247  CALL igljm(ialj3,iatlj3,zgm3(1,2),zam3,ly3,ly3,ly3,ly3,alpha,beta)
248 C
249  ENDIF
250 C
251 C
252  ELSE
253 C
254 C*** Three-dimensional case ************************************
255 C
256 C
257 C Gauss-Lobatto Legendre mesh (suffix M1)
258 C Generate collocation points and weights
259 C
260  CALL zwgll (zgm1(1,1),wxm1,lx1)
261  CALL zwgll (zgm1(1,2),wym1,ly1)
262  CALL zwgll (zgm1(1,3),wzm1,lz1)
263  DO 700 iz=1,lz1
264  DO 700 iy=1,ly1
265  DO 700 ix=1,lx1
266  w3m1(ix,iy,iz)=wxm1(ix)*wym1(iy)*wzm1(iz)
267  700 CONTINUE
268 C
269 C Compute derivative matrices
270 C
271  CALL dgll (dxm1,dxtm1,zgm1(1,1),lx1,lx1)
272  CALL dgll (dym1,dytm1,zgm1(1,2),ly1,ly1)
273  CALL dgll (dzm1,dztm1,zgm1(1,3),lz1,lz1)
274 C
275 C Gauss Legendre mesh (suffix M2)
276 C Generate collocation points and weights
277 C
278  IF(ifsplit)THEN
279  CALL zwgll (zgm2(1,1),wxm2,lx2)
280  CALL zwgll (zgm2(1,2),wym2,ly2)
281  CALL zwgll (zgm2(1,3),wzm2,lz2)
282  ELSE
283  CALL zwgl (zgm2(1,1),wxm2,lx2)
284  CALL zwgl (zgm2(1,2),wym2,ly2)
285  CALL zwgl (zgm2(1,3),wzm2,lz2)
286  ENDIF
287  DO 800 iz=1,lz2
288  DO 800 iy=1,ly2
289  DO 800 ix=1,lx2
290  w3m2(ix,iy,iz)=wxm2(ix)*wym2(iy)*wzm2(iz)
291  800 CONTINUE
292 C
293 C Gauss-Loabtto Legendre mesh (suffix M3).
294 C Generate collocation points and weights.
295 C
296  CALL zwgll (zgm3(1,1),wxm3,lx3)
297  CALL zwgll (zgm3(1,2),wym3,ly3)
298  CALL zwgll (zgm3(1,3),wzm3,lz3)
299  DO 900 iz=1,lz3
300  DO 900 iy=1,ly3
301  DO 900 ix=1,lx3
302  w3m3(ix,iy,iz)=wxm3(ix)*wym3(iy)*wzm3(iz)
303  900 CONTINUE
304 C
305 C Compute derivative matrices
306 C
307  CALL dgll (dxm3,dxtm3,zgm3(1,1),lx3,lx3)
308  CALL dgll (dym3,dytm3,zgm3(1,2),ly3,ly3)
309  CALL dgll (dzm3,dztm3,zgm3(1,3),lz3,lz3)
310 C
311 C Generate interpolation operators for the staggered mesh
312 C
313  CALL igllm (ixm12,ixtm12,zgm1(1,1),zgm2(1,1),lx1,lx2,lx1,lx2)
314  CALL igllm (iym12,iytm12,zgm1(1,2),zgm2(1,2),ly1,ly2,ly1,ly2)
315  CALL igllm (izm12,iztm12,zgm1(1,3),zgm2(1,3),lz1,lz2,lz1,lz2)
316 C
317 C NOTE: The splitting scheme has only one mesh!!!!!
318 C
319  IF (ifsplit) THEN
320  CALL igllm (ixm21,ixtm21,zgm1(1,1),zgm2(1,1),lx1,lx2,lx1,lx2)
321  CALL igllm (iym21,iytm21,zgm1(1,2),zgm2(1,2),ly1,ly2,ly1,ly2)
322  CALL igllm (izm21,iztm21,zgm1(1,3),zgm2(1,3),lz1,lz2,lz1,lz2)
323  ELSE
324  CALL iglm (ixm21,ixtm21,zgm2(1,1),zgm1(1,1),lx2,lx1,lx2,lx1)
325  CALL iglm (iym21,iytm21,zgm2(1,2),zgm1(1,2),ly2,ly1,ly2,ly1)
326  CALL iglm (izm21,iztm21,zgm2(1,3),zgm1(1,3),lz2,lz1,lz2,lz1)
327  ENDIF
328 C
329 C Compute derivative operators for the staggered mesh
330 C
331  IF(ifsplit)THEN
332  CALL copy (dxm12, dxm1, lx1*lx2)
333  CALL copy (dxtm12,dxtm1,lx1*lx2)
334  CALL copy (dym12, dym1, ly1*ly2)
335  CALL copy (dytm12,dytm1,ly1*ly2)
336  CALL copy (dzm12, dzm1, lz1*lz2)
337  CALL copy (dztm12,dztm1,lz1*lz2)
338  ELSE
339  CALL dgllgl (dxm12,dxtm12,zgm1(1,1),zgm2(1,1),ixm12,
340  $ lx1,lx2,lx1,lx2)
341  CALL dgllgl (dym12,dytm12,zgm1(1,2),zgm2(1,2),iym12,
342  $ ly1,ly2,ly1,ly2)
343  CALL dgllgl (dzm12,dztm12,zgm1(1,3),zgm2(1,3),izm12,
344  $ lz1,lz2,lz1,lz2)
345  ENDIF
346 C
347 C Compute interpolation operators for the geometry mesh M3.
348 C
349  CALL igllm (ixm13,ixtm13,zgm1(1,1),zgm3(1,1),lx1,lx3,lx1,lx3)
350  CALL igllm (iym13,iytm13,zgm1(1,2),zgm3(1,2),ly1,ly3,ly1,ly3)
351  CALL igllm (izm13,iztm13,zgm1(1,3),zgm3(1,3),lz1,lz3,lz1,lz3)
352  CALL igllm (ixm31,ixtm31,zgm3(1,1),zgm1(1,1),lx3,lx1,lx3,lx1)
353  CALL igllm (iym31,iytm31,zgm3(1,2),zgm1(1,2),ly3,ly1,ly3,ly1)
354  CALL igllm (izm31,iztm31,zgm3(1,3),zgm1(1,3),lz3,lz1,lz3,lz1)
355 C
356  ENDIF
357 C
358  RETURN
359  END
360  subroutine geom1 (xm3,ym3,zm3)
361 C-----------------------------------------------------------------------
362 C
363 C Routine to generate all elemental geometric data for mesh 1.
364 C
365 C Velocity formulation : global-to-local mapping based on mesh 3
366 C Stress formulation : global-to-local mapping based on mesh 1
367 C
368 C-----------------------------------------------------------------------
369  include 'SIZE'
370  include 'GEOM'
371  include 'INPUT'
372  include 'TSTEP'
373 C
374 C Note : XM3,YM3,ZM3 should come from COMMON /SCRUZ/.
375 C
376  dimension xm3(lx3,ly3,lz3,1)
377  $ , ym3(lx3,ly3,lz3,1)
378  $ , zm3(lx3,ly3,lz3,1)
379 C
380  IF (ifgmsh3 .AND. istep.EQ.0) THEN
381  CALL glmapm3 (xm3,ym3,zm3)
382  ELSE
383  CALL glmapm1
384  ENDIF
385 C
386  CALL geodat1
387 C
388  RETURN
389  END
390  subroutine glmapm3 (xm3,ym3,zm3)
391 C-------------------------------------------------------------------
392 C
393 C Routine to generate mapping data based on mesh 3
394 C (Gauss-Legendre Lobatto meshes).
395 C
396 C XRM3, YRM3, ZRM3 - dx/dr, dy/dr, dz/dr
397 C XSM3, YSM3, ZSM3 - dx/ds, dy/ds, dz/ds
398 C XTM3, YTM3, ZTM3 - dx/dt, dy/dt, dz/dt
399 C RXM3, RYM3, RZM3 - dr/dx, dr/dy, dr/dz
400 C SXM3, SYM3, SZM3 - ds/dx, ds/dy, ds/dz
401 C TXM3, TYM3, TZM3 - dt/dx, dt/dy, dt/dz
402 C JACM3 - Jacobian
403 C
404 C------------------------------------------------------------------
405  include 'SIZE'
406  include 'TOTAL'
407 C
408 C Note : work arrays for mesh 3 in scratch commons will be
409 C changed after exit of routine.
410 C
411  COMMON /scrns/ xrm3(lx3,ly3,lz3,lelt)
412  $ , xsm3(lx3,ly3,lz3,lelt)
413  $ , xtm3(lx3,ly3,lz3,lelt)
414  $ , yrm3(lx3,ly3,lz3,lelt)
415  $ , ysm3(lx3,ly3,lz3,lelt)
416  $ , ytm3(lx3,ly3,lz3,lelt)
417  $ , zrm3(lx3,ly3,lz3,lelt)
418  COMMON /ctmp0/ zsm3(lx3,ly3,lz3,lelt)
419  $ , ztm3(lx3,ly3,lz3,lelt)
420  COMMON /ctmp1/ rxm3(lx3,ly3,lz3,lelt)
421  $ , rym3(lx3,ly3,lz3,lelt)
422  $ , rzm3(lx3,ly3,lz3,lelt)
423  $ , sxm3(lx3,ly3,lz3,lelt)
424  COMMON /scrmg/ sym3(lx3,ly3,lz3,lelt)
425  $ , szm3(lx3,ly3,lz3,lelt)
426  $ , txm3(lx3,ly3,lz3,lelt)
427  $ , tym3(lx3,ly3,lz3,lelt)
428  COMMON /screv/ tzm3(lx3,ly3,lz3,lelt)
429  $ , jacm3(lx3,ly3,lz3,lelt)
430  REAL JACM3
431  dimension xm3(lx3,ly3,lz3,1)
432  $ , ym3(lx3,ly3,lz3,1)
433  $ , zm3(lx3,ly3,lz3,1)
434 C
435 C
436  nxy3 = lx3*ly3
437  nyz3 = ly3*lz3
438  nxyz3 = lx3*ly3*lz3
439  ntot3 = nxyz3*nelt
440  nxyz1 = lx1*ly1*lz1
441  ntot1 = nxyz1*nelt
442 C
443 C
444 C Compute isoparametric partials.
445 C
446 
447  IF (ldim.EQ.2) THEN
448 C
449 C Two-dimensional case
450 C
451  DO 200 iel=1,nelt
452 C
453 C Use the appropriate derivative- and interpolation operator in
454 C the y-direction (= radial direction if axisymmetric).
455 C
456  IF (ifaxis) THEN
457  ly33 = ly3*ly3
458  IF (ifrzer(iel)) THEN
459  CALL copy (dytm3,datm3,ly33)
460  ELSE
461  CALL copy (dytm3,dctm3,ly33)
462  ENDIF
463  ENDIF
464 C
465  CALL mxm(dxm3,lx3,xm3(1,1,1,iel),lx3,xrm3(1,1,1,iel),ly3)
466  CALL mxm(dxm3,lx3,ym3(1,1,1,iel),lx3,yrm3(1,1,1,iel),ly3)
467  CALL mxm(xm3(1,1,1,iel),lx3,dytm3,ly3,xsm3(1,1,1,iel),ly3)
468  CALL mxm(ym3(1,1,1,iel),lx3,dytm3,ly3,ysm3(1,1,1,iel),ly3)
469 C
470  200 CONTINUE
471 C
472  CALL rzero (jacm3,ntot3)
473  CALL addcol3 (jacm3,xrm3,ysm3,ntot3)
474  CALL subcol3 (jacm3,xsm3,yrm3,ntot3)
475 C
476  CALL copy (rxm3,ysm3,ntot3)
477  CALL copy (rym3,xsm3,ntot3)
478  CALL chsign (rym3,ntot3)
479  CALL copy (sxm3,yrm3,ntot3)
480  CALL chsign (sxm3,ntot3)
481  CALL copy (sym3,xrm3,ntot3)
482 C
483  ELSE
484 C
485 C Three-dimensional case
486 C
487  DO 300 iel=1,nelt
488 C
489  CALL mxm(dxm3,lx3,xm3(1,1,1,iel),lx3,xrm3(1,1,1,iel),nyz3)
490  CALL mxm(dxm3,lx3,ym3(1,1,1,iel),lx3,yrm3(1,1,1,iel),nyz3)
491  CALL mxm(dxm3,lx3,zm3(1,1,1,iel),lx3,zrm3(1,1,1,iel),nyz3)
492 C
493  DO 310 iz=1,lz3
494  CALL mxm(xm3(1,1,iz,iel),lx3,dytm3,ly3,xsm3(1,1,iz,iel),ly3)
495  CALL mxm(ym3(1,1,iz,iel),lx3,dytm3,ly3,ysm3(1,1,iz,iel),ly3)
496  CALL mxm(zm3(1,1,iz,iel),lx3,dytm3,ly3,zsm3(1,1,iz,iel),ly3)
497  310 CONTINUE
498 C
499  CALL mxm(xm3(1,1,1,iel),nxy3,dztm3,lz3,xtm3(1,1,1,iel),lz3)
500  CALL mxm(ym3(1,1,1,iel),nxy3,dztm3,lz3,ytm3(1,1,1,iel),lz3)
501  CALL mxm(zm3(1,1,1,iel),nxy3,dztm3,lz3,ztm3(1,1,1,iel),lz3)
502 C
503  300 CONTINUE
504 C
505  CALL rzero (jacm3,ntot3)
506  CALL addcol4 (jacm3,xrm3,ysm3,ztm3,ntot3)
507  CALL addcol4 (jacm3,xtm3,yrm3,zsm3,ntot3)
508  CALL addcol4 (jacm3,xsm3,ytm3,zrm3,ntot3)
509  CALL subcol4 (jacm3,xrm3,ytm3,zsm3,ntot3)
510  CALL subcol4 (jacm3,xsm3,yrm3,ztm3,ntot3)
511  CALL subcol4 (jacm3,xtm3,ysm3,zrm3,ntot3)
512 C
513  CALL ascol5 (rxm3,ysm3,ztm3,ytm3,zsm3,ntot3)
514  CALL ascol5 (rym3,xtm3,zsm3,xsm3,ztm3,ntot3)
515  CALL ascol5 (rzm3,xsm3,ytm3,xtm3,ysm3,ntot3)
516  CALL ascol5 (sxm3,ytm3,zrm3,yrm3,ztm3,ntot3)
517  CALL ascol5 (sym3,xrm3,ztm3,xtm3,zrm3,ntot3)
518  CALL ascol5 (szm3,xtm3,yrm3,xrm3,ytm3,ntot3)
519  CALL ascol5 (txm3,yrm3,zsm3,ysm3,zrm3,ntot3)
520  CALL ascol5 (tym3,xsm3,zrm3,xrm3,zsm3,ntot3)
521  CALL ascol5 (tzm3,xrm3,ysm3,xsm3,yrm3,ntot3)
522 C
523  ENDIF
524 C
525 C Mapping from space P(n-2) to space P(n) (mesh M3 to mesh M1).
526 C
527  IF (ldim.EQ.2) THEN
528  CALL rzero (rzm1,ntot1)
529  CALL rzero (szm1,ntot1)
530  CALL rone (tzm1,ntot1)
531  ENDIF
532 C
533  kerr = 0
534  DO 400 ie=1,nelt
535 
536 c write(6,*) 'chkj1'
537 c call outxm3j(xm3,ym3,jacm3)
538 
539  CALL chkjac(jacm3(1,1,1,ie),nxyz3,ie,xm3(1,1,1,ie),
540  $ ym3(1,1,1,ie),zm3(1,1,1,ie),ldim,ierr)
541  if (ierr.eq.1) kerr = kerr+1
542  CALL map31 (rxm1(1,1,1,ie),rxm3(1,1,1,ie),ie)
543  CALL map31 (rym1(1,1,1,ie),rym3(1,1,1,ie),ie)
544  CALL map31 (sxm1(1,1,1,ie),sxm3(1,1,1,ie),ie)
545  CALL map31 (sym1(1,1,1,ie),sym3(1,1,1,ie),ie)
546  IF (ldim.EQ.3) THEN
547  CALL map31 (rzm1(1,1,1,ie),rzm3(1,1,1,ie),ie)
548  CALL map31 (szm1(1,1,1,ie),szm3(1,1,1,ie),ie)
549  CALL map31 (txm1(1,1,1,ie),txm3(1,1,1,ie),ie)
550  CALL map31 (tym1(1,1,1,ie),tym3(1,1,1,ie),ie)
551  CALL map31 (tzm1(1,1,1,ie),tzm3(1,1,1,ie),ie)
552  ENDIF
553  CALL map31 (jacm1(1,1,1,ie),jacm3(1,1,1,ie),ie)
554  CALL map31 (xm1(1,1,1,ie),xm3(1,1,1,ie),ie)
555  CALL map31 (ym1(1,1,1,ie),ym3(1,1,1,ie),ie)
556  CALL map31 (zm1(1,1,1,ie),zm3(1,1,1,ie),ie)
557  400 CONTINUE
558  kerr = iglsum(kerr,1)
559  if (kerr.gt.0) then
560  ifxyo = .true.
561  ifvo = .false.
562  ifpo = .false.
563  ifto = .false.
564  param(66) = 4
565  call outpost(vx,vy,vz,pr,t,'xyz')
566  if (nid.eq.0) write(6,*) 'Jac error 3, setting p66=4, ifxyo=t'
567  call exitt
568  endif
569 
570  call invers2(jacmi,jacm1,ntot1)
571 
572  RETURN
573  END
574  subroutine glmapm1
575 C-----------------------------------------------------------------------
576 C
577 C Routine to generate mapping data based on mesh 1
578 C (Gauss-Legendre Lobatto meshes).
579 C
580 C XRM1, YRM1, ZRM1 - dx/dr, dy/dr, dz/dr
581 C XSM1, YSM1, ZSM1 - dx/ds, dy/ds, dz/ds
582 C XTM1, YTM1, ZTM1 - dx/dt, dy/dt, dz/dt
583 C RXM1, RYM1, RZM1 - dr/dx, dr/dy, dr/dz
584 C SXM1, SYM1, SZM1 - ds/dx, ds/dy, ds/dz
585 C TXM1, TYM1, TZM1 - dt/dx, dt/dy, dt/dz
586 C JACM1 - Jacobian
587 C
588 C-----------------------------------------------------------------------
589  include 'SIZE'
590  include 'GEOM'
591  include 'INPUT'
592  include 'SOLN'
593 C
594 C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
595 C share the same array structure in Scratch Common /SCRNS/.
596 C
597  COMMON /scrns/ xrm1(lx1,ly1,lz1,lelt)
598  $ , yrm1(lx1,ly1,lz1,lelt)
599  $ , xsm1(lx1,ly1,lz1,lelt)
600  $ , ysm1(lx1,ly1,lz1,lelt)
601  $ , xtm1(lx1,ly1,lz1,lelt)
602  $ , ytm1(lx1,ly1,lz1,lelt)
603  $ , zrm1(lx1,ly1,lz1,lelt)
604  COMMON /ctmp1/ zsm1(lx1,ly1,lz1,lelt)
605  $ , ztm1(lx1,ly1,lz1,lelt)
606 C
607  nxy1 = lx1*ly1
608  nyz1 = ly1*lz1
609  nxyz1 = lx1*ly1*lz1
610  ntot1 = nxyz1*nelt
611 C
612  CALL xyzrst (xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,xtm1,ytm1,ztm1,
613  $ ifaxis)
614 C
615  IF (ldim.EQ.2) THEN
616  CALL rzero (jacm1,ntot1)
617  CALL addcol3 (jacm1,xrm1,ysm1,ntot1)
618  CALL subcol3 (jacm1,xsm1,yrm1,ntot1)
619  CALL copy (rxm1,ysm1,ntot1)
620  CALL copy (rym1,xsm1,ntot1)
621  CALL chsign (rym1,ntot1)
622  CALL copy (sxm1,yrm1,ntot1)
623  CALL chsign (sxm1,ntot1)
624  CALL copy (sym1,xrm1,ntot1)
625  CALL rzero (rzm1,ntot1)
626  CALL rzero (szm1,ntot1)
627  CALL rone (tzm1,ntot1)
628  ELSE
629  CALL rzero (jacm1,ntot1)
630  CALL addcol4 (jacm1,xrm1,ysm1,ztm1,ntot1)
631  CALL addcol4 (jacm1,xtm1,yrm1,zsm1,ntot1)
632  CALL addcol4 (jacm1,xsm1,ytm1,zrm1,ntot1)
633  CALL subcol4 (jacm1,xrm1,ytm1,zsm1,ntot1)
634  CALL subcol4 (jacm1,xsm1,yrm1,ztm1,ntot1)
635  CALL subcol4 (jacm1,xtm1,ysm1,zrm1,ntot1)
636  CALL ascol5 (rxm1,ysm1,ztm1,ytm1,zsm1,ntot1)
637  CALL ascol5 (rym1,xtm1,zsm1,xsm1,ztm1,ntot1)
638  CALL ascol5 (rzm1,xsm1,ytm1,xtm1,ysm1,ntot1)
639  CALL ascol5 (sxm1,ytm1,zrm1,yrm1,ztm1,ntot1)
640  CALL ascol5 (sym1,xrm1,ztm1,xtm1,zrm1,ntot1)
641  CALL ascol5 (szm1,xtm1,yrm1,xrm1,ytm1,ntot1)
642  CALL ascol5 (txm1,yrm1,zsm1,ysm1,zrm1,ntot1)
643  CALL ascol5 (tym1,xsm1,zrm1,xrm1,zsm1,ntot1)
644  CALL ascol5 (tzm1,xrm1,ysm1,xsm1,yrm1,ntot1)
645  ENDIF
646 C
647  kerr = 0
648  DO 500 ie=1,nelt
649  CALL chkjac(jacm1(1,1,1,ie),nxyz1,ie,xm1(1,1,1,ie),
650  $ ym1(1,1,1,ie),zm1(1,1,1,ie),ldim,ierr)
651  if (ierr.ne.0) kerr = kerr+1
652  500 CONTINUE
653  kerr = iglsum(kerr,1)
654  if (kerr.gt.0) then
655  ifxyo = .true.
656  ifvo = .false.
657  ifpo = .false.
658  ifto = .false.
659  param(66) = 4
660  call outpost(vx,vy,vz,pr,t,'xyz')
661  if (nid.eq.0) write(6,*) 'Jac error 1, setting p66=4, ifxyo=t'
662  call exitt
663  endif
664 
665  call invers2(jacmi,jacm1,ntot1)
666 
667  RETURN
668  END
669  subroutine geodat1
670 C-----------------------------------------------------------------------
671 C
672 C Routine to generate elemental geometric matrices on mesh 1
673 C (Gauss-Legendre Lobatto mesh).
674 C
675 C-----------------------------------------------------------------------
676  include 'SIZE'
677  include 'GEOM'
678  include 'INPUT'
679  include 'MASS'
680  include 'TSTEP'
681  include 'WZ'
682 C
683 C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
684 C share the same array structure in Scratch Common /SCRNS/.
685 C
686  COMMON /scrns/ xrm1(lx1,ly1,lz1,lelt)
687  $ , yrm1(lx1,ly1,lz1,lelt)
688  $ , xsm1(lx1,ly1,lz1,lelt)
689  $ , ysm1(lx1,ly1,lz1,lelt)
690  $ , xtm1(lx1,ly1,lz1,lelt)
691  $ , ytm1(lx1,ly1,lz1,lelt)
692  $ , zrm1(lx1,ly1,lz1,lelt)
693  COMMON /ctmp1/ zsm1(lx1,ly1,lz1,lelt)
694  $ , ztm1(lx1,ly1,lz1,lelt)
695  $ , wj(lx1,ly1,lz1,lelt)
696 C
697  nxyz1 = lx1*ly1*lz1
698  ntot1 = nxyz1*nelt
699 C
700  IF (ifgmsh3 .AND. istep.EQ.0)
701  $ CALL xyzrst (xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,xtm1,ytm1,ztm1,
702  $ ifaxis)
703 C
704  IF (.NOT.ifaxis) THEN
705  CALL invers2 (wj,jacm1,ntot1)
706  ELSE
707  DO 500 iel=1,nelt
708  IF (ifrzer(iel)) THEN
709  DO 510 j=1,ly1
710  DO 510 i=1,lx1
711  IF (j.GT.1) THEN
712  wj(i,j,1,iel) = ym1(i,j,1,iel)/
713  $ (jacm1(i,j,1,iel)*(1.+zam1(j)))
714  ELSE
715  wj(i,j,1,iel) = ysm1(i,j,1,iel)/jacm1(i,j,1,iel)
716  ENDIF
717  510 CONTINUE
718  ELSE
719  CALL invcol3 (wj(1,1,1,iel),ym1(1,1,1,iel),
720  $ jacm1(1,1,1,iel),nxyz1)
721  ENDIF
722  500 CONTINUE
723  ENDIF
724 C
725 C Compute geometric factors for integrated del-squared operator.
726 C
727  IF (ldim.EQ.2) THEN
728  CALL vdot2 (g1m1,rxm1,rym1,rxm1,rym1,ntot1)
729  CALL vdot2 (g2m1,sxm1,sym1,sxm1,sym1,ntot1)
730  CALL vdot2 (g4m1,rxm1,rym1,sxm1,sym1,ntot1)
731  CALL col2 (g1m1,wj,ntot1)
732  CALL col2 (g2m1,wj,ntot1)
733  CALL col2 (g4m1,wj,ntot1)
734  CALL rzero (g3m1,ntot1)
735  CALL rzero (g5m1,ntot1)
736  CALL rzero (g6m1,ntot1)
737  ELSE
738  CALL vdot3 (g1m1,rxm1,rym1,rzm1,rxm1,rym1,rzm1,ntot1)
739  CALL vdot3 (g2m1,sxm1,sym1,szm1,sxm1,sym1,szm1,ntot1)
740  CALL vdot3 (g3m1,txm1,tym1,tzm1,txm1,tym1,tzm1,ntot1)
741  CALL vdot3 (g4m1,rxm1,rym1,rzm1,sxm1,sym1,szm1,ntot1)
742  CALL vdot3 (g5m1,rxm1,rym1,rzm1,txm1,tym1,tzm1,ntot1)
743  CALL vdot3 (g6m1,sxm1,sym1,szm1,txm1,tym1,tzm1,ntot1)
744  CALL col2 (g1m1,wj,ntot1)
745  CALL col2 (g2m1,wj,ntot1)
746  CALL col2 (g3m1,wj,ntot1)
747  CALL col2 (g4m1,wj,ntot1)
748  CALL col2 (g5m1,wj,ntot1)
749  CALL col2 (g6m1,wj,ntot1)
750  ENDIF
751 C
752 C Multiply the geometric factors GiM1,i=1,5 with the
753 C weights on mesh M1.
754 C
755  DO 580 iel=1,nelt
756  IF (ifaxis) CALL setaxw1 ( ifrzer(iel) )
757  CALL col2 (g1m1(1,1,1,iel),w3m1,nxyz1)
758  CALL col2 (g2m1(1,1,1,iel),w3m1,nxyz1)
759  CALL col2 (g4m1(1,1,1,iel),w3m1,nxyz1)
760  IF (ldim.EQ.3) THEN
761  CALL col2 (g3m1(1,1,1,iel),w3m1,nxyz1)
762  CALL col2 (g5m1(1,1,1,iel),w3m1,nxyz1)
763  CALL col2 (g6m1(1,1,1,iel),w3m1,nxyz1)
764  ENDIF
765  580 CONTINUE
766 C
767 C Compute the mass matrix on mesh M1.
768 C
769  DO 700 iel=1,nelt
770  IF (ifaxis) CALL setaxw1 ( ifrzer(iel) )
771  CALL col3 (bm1(1,1,1,iel),jacm1(1,1,1,iel),w3m1,nxyz1)
772  IF (ifaxis) THEN
773  CALL col3(baxm1(1,1,1,iel),jacm1(1,1,1,iel),w3m1,nxyz1)
774  IF (ifrzer(iel)) THEN
775  DO 600 j=1,ly1
776  IF (j.GT.1) THEN
777  DO 610 i=1,lx1
778  bm1(i,j,1,iel) = bm1(i,j,1,iel)*ym1(i,j,1,iel)
779  $ /(1.+zam1(j))
780  baxm1(i,j,1,iel)=baxm1(i,j,1,iel)/(1.+zam1(j))
781  610 CONTINUE
782  ELSE
783  DO 620 i=1,lx1
784  bm1(i,j,1,iel) = bm1(i,j,1,iel)*ysm1(i,j,1,iel)
785  baxm1(i,j,1,iel)=baxm1(i,j,1,iel)
786  620 CONTINUE
787  ENDIF
788  600 CONTINUE
789  ELSE
790  CALL col2 (bm1(1,1,1,iel),ym1(1,1,1,iel),nxyz1)
791  ENDIF
792  ENDIF
793 C
794  700 CONTINUE
795 
796  IF(ifaxis) THEN
797  DO iel=1,nelt
798  IF(ifrzer(iel)) THEN
799  DO j=1,ly1
800  DO i=1,lx1
801  IF(j.EQ.1) THEN
802  yinvm1(i,j,1,iel)=1.0d0/ysm1(i,j,1,iel)
803  ELSE
804  yinvm1(i,j,1,iel)=1.0d0/ym1(i,j,1,iel)
805  ENDIF
806  ENDDO
807  ENDDO
808  ELSE
809  CALL invers2(yinvm1(1,1,1,iel),ym1(1,1,1,iel),nxyz1)
810  ENDIF
811  ENDDO
812  ELSE
813  CALL cfill(yinvm1,1.0d0,nxyz1*nelt)
814  ENDIF
815 C
816 C Compute normals, tangents, and areas on elemental surfaces
817 C
818  CALL setarea
819 C
820  RETURN
821  END
822  subroutine geom2
823 C-------------------------------------------------------------------
824 C
825 C Routine to generate all elemental geometric data for mesh 2
826 C (Gauss-Legendre mesh).
827 C
828 C RXM2, RYM2, RZM2 - dr/dx, dr/dy, dr/dz
829 C SXM2, SYM2, SZM2 - ds/dx, ds/dy, ds/dz
830 C TXM2, TYM2, TZM2 - dt/dx, dt/dy, dt/dz
831 C JACM2 - Jacobian
832 C BM2 - Mass matrix
833 C
834 C------------------------------------------------------------------
835  include 'SIZE'
836  include 'TOTAL'
837 C
838  nxyz2 = lx2*ly2*lz2
839  ntot2 = nxyz2*nelv
840 C
841  IF (ifsplit) THEN
842 C
843 C Mesh 1 and 2 are identical
844 C
845  CALL copy (rxm2,rxm1,ntot2)
846  CALL copy (rym2,rym1,ntot2)
847  CALL copy (rzm2,rzm1,ntot2)
848  CALL copy (sxm2,sxm1,ntot2)
849  CALL copy (sym2,sym1,ntot2)
850  CALL copy (szm2,szm1,ntot2)
851  CALL copy (txm2,txm1,ntot2)
852  CALL copy (tym2,tym1,ntot2)
853  CALL copy (tzm2,tzm1,ntot2)
854  CALL copy (jacm2,jacm1,ntot2)
855  CALL copy (bm2,bm1,ntot2)
856 
857  CALL copy (xm2,xm1,ntot2)
858  CALL copy (ym2,ym1,ntot2)
859  CALL copy (zm2,zm1,ntot2)
860 
861  ELSE
862 C
863 C Consistent approximation spaces (UZAWA)
864 C
865  IF (ldim.EQ.2) THEN
866  CALL rzero (rzm2,ntot2)
867  CALL rzero (szm2,ntot2)
868  CALL rone (tzm2,ntot2)
869  ENDIF
870 C
871  DO 1000 iel=1,nelv
872 C
873 C Mapping from mesh M1 to mesh M2
874 C
875  CALL map12 (rxm2(1,1,1,iel),rxm1(1,1,1,iel),iel)
876  CALL map12 (rym2(1,1,1,iel),rym1(1,1,1,iel),iel)
877  CALL map12 (sxm2(1,1,1,iel),sxm1(1,1,1,iel),iel)
878  CALL map12 (sym2(1,1,1,iel),sym1(1,1,1,iel),iel)
879  IF (ldim.EQ.3) THEN
880  CALL map12 (rzm2(1,1,1,iel),rzm1(1,1,1,iel),iel)
881  CALL map12 (szm2(1,1,1,iel),szm1(1,1,1,iel),iel)
882  CALL map12 (txm2(1,1,1,iel),txm1(1,1,1,iel),iel)
883  CALL map12 (tym2(1,1,1,iel),tym1(1,1,1,iel),iel)
884  CALL map12 (tzm2(1,1,1,iel),tzm1(1,1,1,iel),iel)
885  ENDIF
886  CALL map12 (jacm2(1,1,1,iel),jacm1(1,1,1,iel),iel)
887 C
888  CALL map12 (xm2(1,1,1,iel),xm1(1,1,1,iel),iel)
889  CALL map12 (ym2(1,1,1,iel),ym1(1,1,1,iel),iel)
890  CALL map12 (zm2(1,1,1,iel),zm1(1,1,1,iel),iel)
891 C
892 C Compute the mass matrix on mesh M2.
893 C
894  IF (ifaxis) CALL setaxw2 ( ifrzer(iel) )
895  CALL col3 (bm2(1,1,1,iel),w3m2,jacm2(1,1,1,iel),nxyz2)
896 C
897  IF (ifaxis.AND.ifrzer(iel)) THEN
898  DO 300 j=1,ly2
899  DO 300 i=1,lx2
900  bm2(i,j,1,iel) = bm2(i,j,1,iel)*ym2(i,j,1,iel)
901  $ /(1.+zam2(j))
902  300 CONTINUE
903  ELSEIF (ifaxis.AND.(.NOT.ifrzer(iel))) THEN
904  CALL col2 (bm2(1,1,1,iel),ym2(1,1,1,iel),nxyz2)
905  ENDIF
906  1000 CONTINUE
907 C
908  ENDIF
909 C
910 C Compute inverse of mesh 2 mass matrix, pff 3/5/92
911  CALL invers2(bm2inv,bm2,ntot2)
912 C
913  RETURN
914  END
915  subroutine xyzrst (xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,
916  $ XTM1,YTM1,ZTM1,IFAXIS)
917 C-----------------------------------------------------------------------
918 C
919 C Compute global-to-local derivatives on mesh 1.
920 C
921 C-----------------------------------------------------------------------
922  include 'SIZE'
923  include 'GEOM'
924  include 'DXYZ'
925 C
926  dimension xrm1(lx1,ly1,lz1,1),yrm1(lx1,ly1,lz1,1)
927  $ , zrm1(lx1,ly1,lz1,1),xsm1(lx1,ly1,lz1,1)
928  $ , ysm1(lx1,ly1,lz1,1),zsm1(lx1,ly1,lz1,1)
929  $ , xtm1(lx1,ly1,lz1,1),ytm1(lx1,ly1,lz1,1)
930  $ , ztm1(lx1,ly1,lz1,1)
931  LOGICAL IFAXIS
932 C
933  nxy1=lx1*ly1
934  nyz1=ly1*lz1
935 C
936  DO 100 iel=1,nelt
937 C
938  IF (ifaxis) CALL setaxdy ( ifrzer(iel) )
939 C
940  CALL mxm (dxm1,lx1,xm1(1,1,1,iel),lx1,xrm1(1,1,1,iel),nyz1)
941  CALL mxm (dxm1,lx1,ym1(1,1,1,iel),lx1,yrm1(1,1,1,iel),nyz1)
942  CALL mxm (dxm1,lx1,zm1(1,1,1,iel),lx1,zrm1(1,1,1,iel),nyz1)
943 C
944  DO 10 iz=1,lz1
945  CALL mxm (xm1(1,1,iz,iel),lx1,dytm1,ly1,xsm1(1,1,iz,iel),ly1)
946  CALL mxm (ym1(1,1,iz,iel),lx1,dytm1,ly1,ysm1(1,1,iz,iel),ly1)
947  CALL mxm (zm1(1,1,iz,iel),lx1,dytm1,ly1,zsm1(1,1,iz,iel),ly1)
948  10 CONTINUE
949 C
950  IF (ldim.EQ.3) THEN
951  CALL mxm (xm1(1,1,1,iel),nxy1,dztm1,lz1,xtm1(1,1,1,iel),lz1)
952  CALL mxm (ym1(1,1,1,iel),nxy1,dztm1,lz1,ytm1(1,1,1,iel),lz1)
953  CALL mxm (zm1(1,1,1,iel),nxy1,dztm1,lz1,ztm1(1,1,1,iel),lz1)
954  ELSE
955  CALL rzero (xtm1(1,1,1,iel),nxy1)
956  CALL rzero (ytm1(1,1,1,iel),nxy1)
957  CALL rone (ztm1(1,1,1,iel),nxy1)
958  ENDIF
959 C
960  100 CONTINUE
961 C
962  RETURN
963  END
964  subroutine chkjac(jac,n,iel,X,Y,Z,ND,IERR)
965 c
966  include 'SIZE'
967  include 'PARALLEL'
968 C
969 C Check the array JAC for a change in sign.
970 C
971  REAL JAC(N),x(1),y(1),z(1)
972 c
973  ierr = 1
974  sign = jac(1)
975  DO 100 i=2,n
976  IF (sign*jac(i).LE.0.0) THEN
977  ieg = lglel(iel)
978  WRITE(6,101) nid,i,ieg
979  write(6,*) jac(i-1),jac(i)
980  if (ldim.eq.3) then
981  write(6,7) nid,x(i-1),y(i-1),z(i-1)
982  write(6,7) nid,x(i),y(i),z(i)
983  else
984  write(6,7) nid,x(i-1),y(i-1)
985  write(6,7) nid,x(i),y(i)
986  endif
987  7 format(i5,' xyz:',1p3e14.5)
988 c if (np.eq.1) call out_xyz_el(x,y,z,iel)
989 c ierr=0
990  return
991  ENDIF
992  100 CONTINUE
993  101 FORMAT(//,i5,2x
994  $ ,'ERROR: Vanishing Jacobian near',i7,'th node of element'
995  $ ,i10,'.')
996 c
997 c
998  ierr = 0
999  RETURN
1000  END
1001 c-----------------------------------------------------------------------
1002  subroutine volume
1003 C
1004 C Compute the volume based on mesh M1 and mesh M2
1005 C
1006 
1007  include 'SIZE'
1008  include 'ESOLV'
1009  include 'INPUT'
1010  include 'MASS'
1011  include 'TSTEP'
1012  integer e
1013 C
1014  volvm1=glsum(bm1,lx1*ly1*lz1*nelv)
1015  volvm2=glsum(bm2,lx2*ly2*lz2*nelv)
1016  voltm1=glsum(bm1,lx1*ly1*lz1*nelt)
1017  voltm2=glsum(bm2,lx2*ly2*lz2*nelt)
1018  mfield=1
1019  if (ifmvbd) mfield=0
1020  nfldt = nfield
1021  if (ifmhd) nfldt = nfield+1
1022 
1023  do ifld=mfield,nfldt
1024  if (iftmsh(ifld)) then
1025  volfld(ifld) = voltm1
1026  else
1027  volfld(ifld) = volvm1
1028  endif
1029  enddo
1030 
1031 c if (nio.eq.0) write(6,*) 'vol_t,vol_v:',voltm1,volvm1
1032 
1033 
1034  nxyz = lx1*ly1*lz1
1035  do e=1,nelt
1036  volel(e) = vlsum(bm1(1,1,1,e),nxyz)
1037  enddo
1038 
1039  return
1040  end
1041 c-----------------------------------------------------------------------
1042  subroutine setarea
1043 C
1044 C Compute surface data: areas, normals and tangents
1045 C
1046  include 'SIZE'
1047  include 'GEOM'
1048  include 'INPUT'
1049 C
1050  nsrf = 6*lx1*lz1*nelt
1051 C
1052  CALL rzero (area,nsrf)
1053  CALL rzero3 (unx,uny,unz,nsrf)
1054  CALL rzero3 (t1x,t1y,t1z,nsrf)
1055  CALL rzero3 (t2x,t2y,t2z,nsrf)
1056 C
1057  IF (ldim.EQ.2) THEN
1058  CALL area2
1059  ELSE
1060  CALL area3
1061  ENDIF
1062 C
1063  RETURN
1064  END
1065  subroutine area2
1066 C--------------------------------------------------------------------
1067 C
1068 C Compute areas, normals and tangents (2D and Axisymmetric geom.)
1069 C
1070 C--------------------------------------------------------------------
1071  include 'SIZE'
1072  include 'GEOM'
1073 C
1074 C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
1075 C share the same array structure in Scratch Common /SCRNS/.
1076 C
1077  COMMON /scrns/ xrm1(lx1,ly1,lz1,lelt)
1078  $ , yrm1(lx1,ly1,lz1,lelt)
1079  $ , xsm1(lx1,ly1,lz1,lelt)
1080  $ , ysm1(lx1,ly1,lz1,lelt)
1081  COMMON /ctmp0/ wgtr1(lx1,lelt)
1082  $ , wgtr2(ly1,lelt)
1083  $ , wgtr3(lx1,lelt)
1084  $ , wgtr4(ly1,lelt)
1085 C
1086  CALL setwgtr (wgtr1,wgtr2,wgtr3,wgtr4)
1087 C
1088 C "R"
1089 C
1090  DO 100 iel=1,nelt
1091  DO 100 iy=1,ly1
1092  xs2 = xsm1(lx1,iy,1,iel)
1093  ys2 = ysm1(lx1,iy,1,iel)
1094  xs4 = xsm1( 1,iy,1,iel)
1095  ys4 = ysm1( 1,iy,1,iel)
1096  ss2 = sqrt( xs2**2 + ys2**2 )
1097  ss4 = sqrt( xs4**2 + ys4**2 )
1098  t1x(iy,1,2,iel) = xs2 / ss2
1099  t1y(iy,1,2,iel) = ys2 / ss2
1100  t1x(iy,1,4,iel) = -xs4 / ss4
1101  t1y(iy,1,4,iel) = -ys4 / ss4
1102  unx(iy,1,2,iel) = t1y(iy,1,2,iel)
1103  uny(iy,1,2,iel) = -t1x(iy,1,2,iel)
1104  unx(iy,1,4,iel) = t1y(iy,1,4,iel)
1105  uny(iy,1,4,iel) = -t1x(iy,1,4,iel)
1106  area(iy,1,2,iel) = ss2 * wgtr2(iy,iel)
1107  area(iy,1,4,iel) = ss4 * wgtr4(iy,iel)
1108  100 CONTINUE
1109 C
1110 C "S"
1111 C
1112  DO 200 iel=1,nelt
1113  DO 200 ix=1,lx1
1114  xr1 = xrm1(ix, 1,1,iel)
1115  yr1 = yrm1(ix, 1,1,iel)
1116  xr3 = xrm1(ix,ly1,1,iel)
1117  yr3 = yrm1(ix,ly1,1,iel)
1118  rr1 = sqrt( xr1**2 + yr1**2 )
1119  rr3 = sqrt( xr3**2 + yr3**2 )
1120  t1x(ix,1,1,iel) = xr1 / rr1
1121  t1y(ix,1,1,iel) = yr1 / rr1
1122  t1x(ix,1,3,iel) = -xr3 / rr3
1123  t1y(ix,1,3,iel) = -yr3 / rr3
1124  unx(ix,1,1,iel) = t1y(ix,1,1,iel)
1125  uny(ix,1,1,iel) = -t1x(ix,1,1,iel)
1126  unx(ix,1,3,iel) = t1y(ix,1,3,iel)
1127  uny(ix,1,3,iel) = -t1x(ix,1,3,iel)
1128  area(ix,1,1,iel) = rr1 * wgtr1(ix,iel)
1129  area(ix,1,3,iel) = rr3 * wgtr3(ix,iel)
1130  200 CONTINUE
1131 C
1132  RETURN
1133  END
1134  subroutine setwgtr (wgtr1,wgtr2,wgtr3,wgtr4)
1135 C
1136  include 'SIZE'
1137  include 'GEOM'
1138  include 'INPUT'
1139  include 'WZ'
1140 C
1141 C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
1142 C share the same array structure in Scratch Common /SCRNS/.
1143 C
1144  COMMON /scrns/ xrm1(lx1,ly1,lz1,lelt)
1145  $ , yrm1(lx1,ly1,lz1,lelt)
1146  $ , xsm1(lx1,ly1,lz1,lelt)
1147  $ , ysm1(lx1,ly1,lz1,lelt)
1148 C
1149  dimension wgtr1(lx1,1)
1150  $ , wgtr2(ly1,1)
1151  $ , wgtr3(lx1,1)
1152  $ , wgtr4(ly1,1)
1153 C
1154  IF (ifaxis) THEN
1155  DO 100 iel=1,nelt
1156  DO 120 ix=1,lx1
1157  wgtr1(ix,iel) = ym1(ix, 1,1,iel) * wxm1(ix)
1158  wgtr3(ix,iel) = ym1(ix,ly1,1,iel) * wxm1(ix)
1159  120 CONTINUE
1160  IF ( ifrzer(iel) ) THEN
1161  iy = 1
1162  wgtr2(iy,iel) = ysm1(lx1,iy,1,iel) * wam1(iy)
1163  wgtr4(iy,iel) = ysm1( 1,iy,1,iel) * wam1(iy)
1164  DO 160 iy=2,ly1
1165  dnr = 1. + zam1(iy)
1166  wgtr2(iy,iel) = ym1(lx1,iy,1,iel) * wam1(iy) / dnr
1167  wgtr4(iy,iel) = ym1( 1,iy,1,iel) * wam1(iy) / dnr
1168  160 CONTINUE
1169  ELSE
1170  DO 180 iy=1,ly1
1171  wgtr2(iy,iel) = ym1(lx1,iy,1,iel) * wym1(iy)
1172  wgtr4(iy,iel) = ym1( 1,iy,1,iel) * wym1(iy)
1173  180 CONTINUE
1174  ENDIF
1175  100 CONTINUE
1176  ELSE
1177  DO 200 iel=1,nelt
1178  CALL copy (wgtr1(1,iel),wxm1,lx1)
1179  CALL copy (wgtr2(1,iel),wym1,ly1)
1180  CALL copy (wgtr3(1,iel),wxm1,lx1)
1181  CALL copy (wgtr4(1,iel),wym1,ly1)
1182  200 CONTINUE
1183  ENDIF
1184 C
1185  RETURN
1186  END
1187  subroutine area3
1188 C--------------------------------------------------------------------
1189 C
1190 C Compute areas, normals and tangents (3D geom.)
1191 C
1192 C--------------------------------------------------------------------
1193  include 'SIZE'
1194  include 'WZ'
1195  include 'GEOM'
1196 C
1197 C Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
1198 C share the same array structure in Scratch Common /SCRNS/.
1199 C
1200  COMMON /scrns/ xrm1(lx1,ly1,lz1,lelt)
1201  $ , yrm1(lx1,ly1,lz1,lelt)
1202  $ , xsm1(lx1,ly1,lz1,lelt)
1203  $ , ysm1(lx1,ly1,lz1,lelt)
1204  $ , xtm1(lx1,ly1,lz1,lelt)
1205  $ , ytm1(lx1,ly1,lz1,lelt)
1206  $ , zrm1(lx1,ly1,lz1,lelt)
1207  COMMON /ctmp1/ zsm1(lx1,ly1,lz1,lelt)
1208  $ , ztm1(lx1,ly1,lz1,lelt)
1209  $ , a(lx1,ly1,lz1,lelt)
1210  $ , b(lx1,ly1,lz1,lelt)
1211  COMMON /ctmp0/ c(lx1,ly1,lz1,lelt)
1212  $ , dot(lx1,ly1,lz1,lelt)
1213 C
1214  nxy1 = lx1*ly1
1215  nface = 2*ldim
1216  ntot = lx1*ly1*lz1*nelt
1217  nsrf = 6*lx1*ly1*nelt
1218 C
1219 C "R"
1220 C
1221  CALL vcross(a,b,c,xsm1,ysm1,zsm1,xtm1,ytm1,ztm1,ntot)
1222  CALL vdot3 (dot,a,b,c,a,b,c,ntot)
1223 C
1224  DO 100 iel=1,nelt
1225  DO 100 iz=1,lz1
1226  DO 100 iy=1,ly1
1227  weight = wym1(iy)*wzm1(iz)
1228  area(iy,iz,2,iel) = sqrt(dot(lx1,iy,iz,iel))*weight
1229  area(iy,iz,4,iel) = sqrt(dot( 1,iy,iz,iel))*weight
1230  unx(iy,iz,4,iel) = -a( 1,iy,iz,iel)
1231  unx(iy,iz,2,iel) = a(lx1,iy,iz,iel)
1232  uny(iy,iz,4,iel) = -b( 1,iy,iz,iel)
1233  uny(iy,iz,2,iel) = b(lx1,iy,iz,iel)
1234  unz(iy,iz,4,iel) = -c( 1,iy,iz,iel)
1235  unz(iy,iz,2,iel) = c(lx1,iy,iz,iel)
1236  100 CONTINUE
1237 C
1238 C "S"
1239 C
1240  CALL vcross(a,b,c,xrm1,yrm1,zrm1,xtm1,ytm1,ztm1,ntot)
1241  CALL vdot3 (dot,a,b,c,a,b,c,ntot)
1242  DO 200 iel=1,nelt
1243  DO 200 iz=1,lz1
1244  DO 200 ix=1,lx1
1245  weight=wxm1(ix)*wzm1(iz)
1246  area(ix,iz,1,iel) = sqrt(dot(ix, 1,iz,iel))*weight
1247  area(ix,iz,3,iel) = sqrt(dot(ix,ly1,iz,iel))*weight
1248  unx(ix,iz,1,iel) = a(ix, 1,iz,iel)
1249  unx(ix,iz,3,iel) = -a(ix,ly1,iz,iel)
1250  uny(ix,iz,1,iel) = b(ix, 1,iz,iel)
1251  uny(ix,iz,3,iel) = -b(ix,ly1,iz,iel)
1252  unz(ix,iz,1,iel) = c(ix, 1,iz,iel)
1253  unz(ix,iz,3,iel) = -c(ix,ly1,iz,iel)
1254  200 CONTINUE
1255 C
1256 C "T"
1257 C
1258  CALL vcross(a,b,c,xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,ntot)
1259  CALL vdot3 (dot,a,b,c,a,b,c,ntot)
1260  DO 300 iel=1,nelt
1261  DO 300 ix=1,lx1
1262  DO 300 iy=1,ly1
1263  weight=wxm1(ix)*wym1(iy)
1264  area(ix,iy,5,iel) = sqrt(dot(ix,iy, 1,iel))*weight
1265  area(ix,iy,6,iel) = sqrt(dot(ix,iy,lz1,iel))*weight
1266  unx(ix,iy,5,iel) = -a(ix,iy, 1,iel)
1267  unx(ix,iy,6,iel) = a(ix,iy,lz1,iel)
1268  uny(ix,iy,5,iel) = -b(ix,iy, 1,iel)
1269  uny(ix,iy,6,iel) = b(ix,iy,lz1,iel)
1270  unz(ix,iy,5,iel) = -c(ix,iy, 1,iel)
1271  unz(ix,iy,6,iel) = c(ix,iy,lz1,iel)
1272  300 CONTINUE
1273 C
1274  CALL unitvec (unx,uny,unz,nsrf)
1275 C
1276 C COMPUTE UNIT TANGENT T1
1277 C
1278  DO 600 iel=1,nelt
1279  DO 600 ifc=1,nface
1280  IF (ifc.EQ.1 .OR. ifc.EQ.6) THEN
1281  CALL facexv (t1x(1,1,ifc,iel),t1y(1,1,ifc,iel),
1282  $ t1z(1,1,ifc,iel),
1283  $ xrm1(1,1,1,iel),yrm1(1,1,1,iel),
1284  $ zrm1(1,1,1,iel),ifc,0)
1285  ELSEIF (ifc.EQ.2 .OR. ifc.EQ.5) THEN
1286  CALL facexv (t1x(1,1,ifc,iel),t1y(1,1,ifc,iel),
1287  $ t1z(1,1,ifc,iel),
1288  $ xsm1(1,1,1,iel),ysm1(1,1,1,iel),
1289  $ zsm1(1,1,1,iel),ifc,0)
1290  ELSE
1291  CALL facexv (t1x(1,1,ifc,iel),t1y(1,1,ifc,iel),
1292  $ t1z(1,1,ifc,iel),
1293  $ xtm1(1,1,1,iel),ytm1(1,1,1,iel),
1294  $ ztm1(1,1,1,iel),ifc,0)
1295  ENDIF
1296  600 CONTINUE
1297 C
1298  CALL unitvec (t1x,t1y,t1z,nsrf)
1299 C
1300 C COMPUTE UNIT TANGENT T2 ( T2 = Normal X T1 )
1301 C
1302  DO 700 iel=1,nelt
1303  DO 700 ifc=1,nface
1304  CALL vcross (t2x(1,1,ifc,iel),t2y(1,1,ifc,iel),
1305  $ t2z(1,1,ifc,iel),
1306  $ unx(1,1,ifc,iel),uny(1,1,ifc,iel),
1307  $ unz(1,1,ifc,iel),
1308  $ t1x(1,1,ifc,iel),t1y(1,1,ifc,iel),
1309  $ t1z(1,1,ifc,iel),nxy1)
1310  700 CONTINUE
1311 C
1312  RETURN
1313  END
1314  subroutine lagmass
1315 C--------------------------------------------------------------------
1316 C
1317 C Lag the mass matrix (matrices)
1318 C
1319 C--------------------------------------------------------------------
1320  include 'SIZE'
1321  include 'MASS'
1322  include 'TSTEP'
1323 C
1324  ntot1 = lx1*ly1*lz1*nelt
1325  DO 100 ilag=nbdinp-1,2,-1
1326  CALL copy (bm1lag(1,1,1,1,ilag),bm1lag(1,1,1,1,ilag-1),ntot1)
1327  100 CONTINUE
1328  CALL copy (bm1lag(1,1,1,1,1),bm1,ntot1)
1329 C
1330  RETURN
1331  END
1332 C
1333  subroutine setinvm
1334 C--------------------------------------------------------------------
1335 C
1336 C Invert the mass matrix.
1337 C
1338 C 1) Copy BM1 to BINVM1
1339 C 2) Perform direct stiffness summation on BINVM1
1340 C 3) Compute BINVM1 = 1/BINVM1
1341 C 4) Two inverse mass matrices required because of difference
1342 C in DSSUM routine for IMESH=1 and IMESH=2.
1343 C
1344 C--------------------------------------------------------------------
1345  include 'SIZE'
1346  include 'MASS'
1347  include 'GEOM'
1348  include 'INPUT'
1349  include 'TSTEP'
1350  include 'WZ'
1351 
1352  nxyz1 = lx1*ly1*lz1
1353 
1354  ifld = ifield
1355 
1356 csk IF (IFFLOW) THEN ! Velocity mass matrix
1357  ifield = 1
1358  ntot = nxyz1*nelv
1359  CALL copy (binvm1,bm1,ntot)
1360  CALL dssum (binvm1,lx1,ly1,lz1)
1361  CALL invcol1 (binvm1,ntot)
1362 csk ENDIF
1363 
1364 
1365  IF (ifheat) THEN ! Temperature mass matrix
1366  ifield = 2
1367  ntot = nxyz1*nelt
1368  CALL copy (bintm1,bm1,ntot)
1369  CALL dssum (bintm1,lx1,ly1,lz1)
1370  CALL invcol1 (bintm1,ntot)
1371  ENDIF
1372 
1373  ifield = ifld
1374 
1375  return
1376  end
1377 
1378 c-----------------------------------------------------------------------
1379 
1380  subroutine maprs(y,x,xa,nrest,iel)
1381 C
1382 C Map the elemental array X from Restart mesh to Y on mesh M1
1383 C Conforming elements, i.e. lx1=ly1=lz1.
1384 C
1385 C---------------------------------------------------------------
1386 C
1387  include 'SIZE'
1388  include 'GEOM'
1389  include 'IXYZ'
1390  include 'WZ'
1391  include 'INPUT'
1392 C
1393  REAL X(NREST,NREST,NREST)
1394  REAL Y(LX1,LY1,LZ1)
1395 C
1396  REAL XA(lx1,NREST,NREST)
1397  COMMON /ctmp0/ xb(lx1,ly1,lz1)
1398 C
1399  REAL IXRES(LX1,LX1),IXTRES(LX1,LX1)
1400  REAL IYRES(LY1,LY1),IYTRES(LY1,LY1)
1401  REAL IZRES(LZ1,LZ1),IZTRES(LZ1,LZ1)
1402  REAL ZCRES(20),WCRES(20)
1403  REAL ZARES(20),WARES(20)
1404 C
1405  nzrest = nrest
1406  IF(lz1.EQ.1) nzrest=1
1407  nyzres = nrest*nzrest
1408  nxy1 = lx1 *ly1
1409 C
1410  CALL zwgll (zcres,wcres,nrest)
1411  CALL igllm (ixres,ixtres,zcres,zgm1,nrest,lx1,nrest,lx1)
1412  IF (.NOT.ifaxis) THEN
1413  CALL copy (iyres,ixres,lx1*nrest)
1414  CALL copy (iytres,ixtres,lx1*nrest)
1415  CALL copy (izres,ixres,lx1*nrest)
1416  CALL copy (iztres,ixtres,lx1*nrest)
1417  ELSE
1418 C
1419 C Use the appropriate derivative- and interpolation operator in
1420 C the y-direction (= radial direction if axisymmetric).
1421 C
1422  IF (ifrzer(iel)) THEN
1423  alpha = 0.
1424  beta = 1.
1425  CALL zwglj (zares,wares,nrest,alpha,beta)
1426  CALL igljm (iyres,iytres,zares,zgm1,nrest,ly1,nrest,ly1,
1427  $ alpha,beta)
1428  ly1r = ly1*nrest
1429  ELSE
1430  CALL copy (iyres,ixres,lx1*nrest)
1431  CALL copy (iytres,ixtres,lx1*nrest)
1432  ENDIF
1433  ENDIF
1434 C
1435  IF (ldim.EQ.2) THEN
1436  CALL mxm (ixres,lx1,x,nrest,xa,nrest)
1437  CALL mxm (xa,lx1,iytres,nrest,y,ly1)
1438  ELSE
1439  CALL mxm (ixres,lx1,x,nrest,xa,nyzres)
1440  DO 100 iz=1,nzrest
1441  CALL mxm (xa(1,1,iz),lx1,iytres,nrest,xb(1,1,iz),ly1)
1442  100 CONTINUE
1443  CALL mxm (xb,nxy1,iztres,nzrest,y,lz1)
1444  ENDIF
1445 C
1446  RETURN
1447  END
1448 C
1449  subroutine map31 (y,x,iel)
1450 C---------------------------------------------------------------
1451 C
1452 C Map the elemental array X from mesh M3 to mesh M1
1453 C
1454 C---------------------------------------------------------------
1455 C
1456  include 'SIZE'
1457  include 'GEOM'
1458  include 'IXYZ'
1459  include 'INPUT'
1460 C
1461  REAL X(LX3,LY3,LZ3)
1462  REAL Y(LX1,LY1,LZ1)
1463 C
1464  COMMON /ctmp0/ xa(lx1,ly3,lz3), xb(lx1,ly1,lz3)
1465 C
1466  nyz3 = ly3*lz3
1467  nxy1 = lx1*ly1
1468 C
1469 C Use the appropriate derivative- and interpolation operator in
1470 C the y-direction (= radial direction if axisymmetric).
1471 C
1472  IF (ifaxis) THEN
1473  ly31 = ly1*ly3
1474  IF (ifrzer(iel)) CALL copy (iytm31,iatm31,ly31)
1475  IF (.NOT.ifrzer(iel)) CALL copy (iytm31,ictm31,ly31)
1476  ENDIF
1477 C
1478  IF (if3d) THEN
1479  CALL mxm (ixm31,lx1,x,lx3,xa,nyz3)
1480  DO 100 iz=1,lz3
1481  CALL mxm (xa(1,1,iz),lx1,iytm31,ly3,xb(1,1,iz),ly1)
1482  100 CONTINUE
1483  CALL mxm (xb,nxy1,iztm31,lz3,y,lz1)
1484  ELSE
1485  CALL mxm (ixm31,lx1,x,lx3,xa,nyz3)
1486  CALL mxm (xa,lx1,iytm31,ly3,y,ly1)
1487  ENDIF
1488 C
1489  RETURN
1490  END
1491 C
1492  subroutine map13 (y,x,iel)
1493 C---------------------------------------------------------------
1494 C
1495 C Map the elemental array X from mesh M1 to mesh M3
1496 C
1497 C---------------------------------------------------------------
1498 C
1499  include 'SIZE'
1500  include 'GEOM'
1501  include 'IXYZ'
1502  include 'INPUT'
1503 C
1504  REAL X(LX1,LY1,LZ1)
1505  REAL Y(LX3,LY3,LZ3)
1506 C
1507  COMMON /ctmp0/ xa(lx3,ly1,lz1), xb(lx3,ly3,lz1)
1508 C
1509  nyz1 = ly1*lz1
1510  nxy3 = lx3*ly3
1511 C
1512 C Use the appropriate derivative- and interpolation operator in
1513 C the y-direction (= radial direction if axisymmetric).
1514 C
1515  IF (ifaxis) THEN
1516  ly13 = ly1*ly3
1517  IF (ifrzer(iel)) CALL copy (iytm13,iatm13,ly13)
1518  IF (.NOT.ifrzer(iel)) CALL copy (iytm13,ictm13,ly13)
1519  ENDIF
1520 C
1521  CALL mxm (ixm13,lx3,x,lx1,xa,nyz1)
1522  DO 100 iz=1,lz1
1523  CALL mxm (xa(1,1,iz),lx3,iytm13,ly1,xb(1,1,iz),ly3)
1524  100 CONTINUE
1525  CALL mxm (xb,nxy3,iztm13,lz1,y,lz3)
1526 C
1527  RETURN
1528  END
1529  subroutine map12 (y,x,iel)
1530 C---------------------------------------------------------------
1531 C
1532 C Map the elemental array X from mesh M1 to mesh M2
1533 C
1534 C---------------------------------------------------------------
1535 C
1536  include 'SIZE'
1537  include 'GEOM'
1538  include 'IXYZ'
1539  include 'INPUT'
1540 C
1541  REAL X(LX1,LY1,LZ1)
1542  REAL Y(LX2,LY2,LZ2)
1543 C
1544  COMMON /ctmp00/ xa(lx2,ly1,lz1), xb(lx2,ly2,lz1)
1545 C
1546  nyz1 = ly1*lz1
1547  nxy2 = lx2*ly2
1548 C
1549 C Use the appropriate derivative- and interpolation operator in
1550 C the y-direction (= radial direction if axisymmetric).
1551 C
1552  IF (ifaxis) THEN
1553  ly12 = ly1*ly2
1554  IF (ifrzer(iel)) CALL copy (iytm12,iatm12,ly12)
1555  IF (.NOT.ifrzer(iel)) CALL copy (iytm12,ictm12,ly12)
1556  ENDIF
1557 C
1558  CALL mxm (ixm12,lx2,x,lx1,xa,nyz1)
1559  DO 100 iz=1,lz1
1560  CALL mxm (xa(1,1,iz),lx2,iytm12,ly1,xb(1,1,iz),ly2)
1561  100 CONTINUE
1562  CALL mxm (xb,nxy2,iztm12,lz1,y,lz2)
1563 C
1564  RETURN
1565  END
1566 C
1567  subroutine map21t (y,x,iel)
1568 C---------------------------------------------------------------
1569 C
1570 C Map the elemental array X from mesh M2 to mesh M1 (Y)
1571 C
1572 C---------------------------------------------------------------
1573 C
1574  include 'SIZE'
1575  include 'GEOM'
1576  include 'IXYZ'
1577  include 'INPUT'
1578 C
1579  REAL X(LX2,LY2,LZ2)
1580  REAL Y(LX1,LY1,LZ1)
1581 C
1582  COMMON /ctmp0/ xa(lx1,ly2,lz2), xb(lx1,ly1,lz2)
1583 C
1584  nyz2 = ly2*lz2
1585  nxy1 = lx1*ly1
1586  nxyz = lx1*ly1*lz1
1587 C
1588 C Use the appropriate derivative- and interpolation operator in
1589 C the y-direction (= radial direction if axisymmetric).
1590 C
1591  IF (ifsplit) THEN
1592  CALL copy(y,x,nxyz)
1593  RETURN
1594  ENDIF
1595 C
1596  IF (if3d) THEN
1597  CALL mxm (ixm21,lx1,x,lx2,xa,nyz2)
1598  DO 100 iz=1,lz2
1599  CALL mxm (xa(1,1,iz),lx1,iytm21,ly2,xb(1,1,iz),ly1)
1600  100 CONTINUE
1601  CALL mxm (xb,nxy1,iztm21,lz2,y,lz1)
1602  ELSE
1603  CALL mxm (ixm21,lx1,x,lx2,xa,nyz2)
1604  CALL mxm (xa,lx1,iytm21,ly2,y,ly1)
1605  ENDIF
1606  RETURN
1607  END
1608  subroutine map21e (y,x,iel)
1609 C---------------------------------------------------------------
1610 C
1611 C Map the elemental array X from mesh M2 to mesh M1
1612 C
1613 C---------------------------------------------------------------
1614 C
1615  include 'SIZE'
1616  include 'GEOM'
1617  include 'IXYZ'
1618  include 'INPUT'
1619 C
1620  REAL X(LX2,LY2,LZ2)
1621  REAL Y(LX1,LY1,LZ1)
1622 C
1623  COMMON /ctmp0/ xa(lx1,ly2,lz2), xb(lx1,ly1,lz2)
1624 C
1625  nyz2 = ly2*lz2
1626  nxy1 = lx1*ly1
1627 C
1628 C Use the appropriate derivative- and interpolation operator in
1629 C the y-direction (= radial direction if axisymmetric).
1630 C
1631  IF (ifaxis) THEN
1632  ly21 = ly1*ly2
1633  IF (ifrzer(iel)) CALL copy (iym12,iam12,ly21)
1634  IF (.NOT.ifrzer(iel)) CALL copy (iym12,icm12,ly21)
1635  ENDIF
1636 C
1637  CALL mxm (ixtm12,lx1,x,lx2,xa,nyz2)
1638  DO 100 iz=1,lz2
1639  CALL mxm (xa(1,1,iz),lx1,iym12,ly2,xb(1,1,iz),ly1)
1640  100 CONTINUE
1641  CALL mxm (xb,nxy1,izm12,lz2,y,lz1)
1642 C
1643  RETURN
1644  END
1645 c-----------------------------------------------------------------------
1646  subroutine out_xyz_el(x,y,z,e)
1647  include 'SIZE'
1648  integer e
1649  real x(1),y(1),z(1)
1650 c
1651  call out_fld_el(x,e,'XQ')
1652  call out_fld_el(y,e,'YQ')
1653  call out_fld_el(z,e,'ZQ')
1654 c
1655  return
1656  end
1657 c-----------------------------------------------------------------------
1658  subroutine out_fld_el(x,e,c2)
1659  include 'SIZE'
1660  real x(lx1,ly1,lz1,lelt)
1661  integer e
1662  character*2 c2
1663 c
1664  write(6,1) c2,e
1665  nx8 = min(lx1,8)
1666  do k=1,lz1
1667  do j=1,ly1
1668  write(6,1) c2,e,(x(i,j,k,e),i=1,nx8)
1669  enddo
1670  enddo
1671  1 format(a2,i6,1p8e11.3)
1672  return
1673  end
1674 c-----------------------------------------------------------------------
1675  subroutine outxm3j(xm3,ym3,jm3)
1676  include 'SIZE'
1677  include 'TOTAL'
1678 
1679  real xm3(lx1,ly1,lz1,lelv)
1680  real ym3(lx1,ly1,lz1,lelv)
1681  real jm3(lx1,ly1,lz1,lelv)
1682 
1683  integer e
1684 
1685  do e=1,nelt
1686  write(6,*) e,nelfld(e),iftmsh(e),' iftmsh'
1687  call outmat(xm3(1,1,1,e),lx3,ly3,' xm3 ',e)
1688  call outmat(ym3(1,1,1,e),lx3,ly3,' ym3 ',e)
1689  call outmat(jm3(1,1,1,e),lx3,ly3,' jm3 ',e)
1690  enddo
1691 
1692  return
1693  end
1694 c-----------------------------------------------------------------------
1695  SUBROUTINE invmt(A,B,AA,N)
1696 C
1697  REAL A(N,N),AA(N,N),B(N,N)
1698  INTEGER INDX(100)
1699 C
1700  nn = n*n
1701  DO 12 i=1,n
1702  DO 11 j=1,n
1703  b(i,j) = 0.0
1704  11 CONTINUE
1705  b(i,i) = 1.0
1706  12 CONTINUE
1707 C
1708  CALL copy (aa,a,nn)
1709  CALL ludcmp(aa,n,n,indx,d)
1710  DO 13 j=1,n
1711  CALL lubksb(aa,n,n,indx,b(1,j))
1712  13 CONTINUE
1713 C
1714  RETURN
1715  END
1716 
1717  SUBROUTINE lubksb(A,N,NP,INDX,B)
1718  REAL A(NP,NP),B(N)
1719  INTEGER INDX(N)
1720  ii=0
1721  DO 12 i=1,n
1722  ll=indx(i)
1723  sum=b(ll)
1724  b(ll)=b(i)
1725  IF (ii.NE.0)THEN
1726  DO 11 j=ii,i-1
1727  sum=sum-a(i,j)*b(j)
1728 11 CONTINUE
1729  ELSE IF (sum.NE.0.0) THEN
1730  ii=i
1731  ENDIF
1732  b(i)=sum
1733 12 CONTINUE
1734  DO 14 i=n,1,-1
1735  sum=b(i)
1736  IF(i.LT.n)THEN
1737  DO 13 j=i+1,n
1738  sum=sum-a(i,j)*b(j)
1739 13 CONTINUE
1740  ENDIF
1741  b(i)=sum/a(i,i)
1742 14 CONTINUE
1743  RETURN
1744  END
1745  SUBROUTINE ludcmp(A,N,NP,INDX,D)
1746  parameter(nmax=100,tiny=1.0e-20)
1747  REAL A(NP,NP),VV(NMAX)
1748  INTEGER INDX(N)
1749  d=1.0
1750  DO 12 i=1,n
1751  aamax=0.0
1752  DO 11 j=1,n
1753  IF (abs(a(i,j)).GT.aamax) aamax=abs(a(i,j))
1754 11 CONTINUE
1755  IF (aamax.EQ.0.0) THEN
1756  write(6,*) 'Singular matrix.'
1757  call exitt
1758  ENDIF
1759  vv(i)=1.0/aamax
1760 12 CONTINUE
1761  DO 19 j=1,n
1762  IF (j.GT.1) THEN
1763  DO 14 i=1,j-1
1764  sum=a(i,j)
1765  IF (i.GT.1)THEN
1766  DO 13 k=1,i-1
1767  sum=sum-a(i,k)*a(k,j)
1768 13 CONTINUE
1769  a(i,j)=sum
1770  ENDIF
1771 14 CONTINUE
1772  ENDIF
1773  aamax=0.0
1774  DO 16 i=j,n
1775  sum=a(i,j)
1776  IF (j.GT.1)THEN
1777  DO 15 k=1,j-1
1778  sum=sum-a(i,k)*a(k,j)
1779 15 CONTINUE
1780  a(i,j)=sum
1781  ENDIF
1782  dum=vv(i)*abs(sum)
1783  IF (dum.GE.aamax) THEN
1784  imax=i
1785  aamax=dum
1786  ENDIF
1787 16 CONTINUE
1788  IF (j.NE.imax)THEN
1789  DO 17 k=1,n
1790  dum=a(imax,k)
1791  a(imax,k)=a(j,k)
1792  a(j,k)=dum
1793 17 CONTINUE
1794  d=-d
1795  vv(imax)=vv(j)
1796  ENDIF
1797  indx(j)=imax
1798  IF(j.NE.n)THEN
1799  IF(a(j,j).EQ.0.)a(j,j)=tiny
1800  dum=1.0/a(j,j)
1801  DO 18 i=j+1,n
1802  a(i,j)=a(i,j)*dum
1803 18 CONTINUE
1804  ENDIF
1805 19 CONTINUE
1806  IF(a(n,n).EQ.0.0)a(n,n)=tiny
1807  RETURN
1808  END
1809 c-----------------------------------------------------------------------
1810  subroutine set_unr
1811  include 'SIZE'
1812  include 'INPUT'
1813  include 'GEOM'
1814  include 'TOPOL'
1815 
1816  integer e,f,pf
1817 
1818  nface = 2*ldim
1819  call dsset(lx1,ly1,lz1)
1820 
1821  do e=1,nelt
1822  do f=1,nface
1823 
1824  pf = eface1(f)
1825  js1 = skpdat(1,pf)
1826  jf1 = skpdat(2,pf)
1827  jskip1 = skpdat(3,pf)
1828  js2 = skpdat(4,pf)
1829  jf2 = skpdat(5,pf)
1830  jskip2 = skpdat(6,pf)
1831 
1832  i = 0
1833  if (if3d) then
1834  do j2=js2,jf2,jskip2
1835  do j1=js1,jf1,jskip1
1836  i = i+1
1837  a = area(i,1,f,e)/jacm1(j1,j2,1,e)
1838  unr(i,f,e) = a * ( rxm1(j1,j2,1,e)*unx(i,1,f,e)
1839  $ + rym1(j1,j2,1,e)*uny(i,1,f,e)
1840  $ + rzm1(j1,j2,1,e)*unz(i,1,f,e) )
1841  uns(i,f,e) = a * ( sxm1(j1,j2,1,e)*unx(i,1,f,e)
1842  $ + sym1(j1,j2,1,e)*uny(i,1,f,e)
1843  $ + szm1(j1,j2,1,e)*unz(i,1,f,e) )
1844  unt(i,f,e) = a * ( txm1(j1,j2,1,e)*unx(i,1,f,e)
1845  $ + tym1(j1,j2,1,e)*uny(i,1,f,e)
1846  $ + tzm1(j1,j2,1,e)*unz(i,1,f,e) )
1847  enddo
1848  enddo
1849  else
1850  do j2=js2,jf2,jskip2
1851  do j1=js1,jf1,jskip1
1852  i = i+1
1853  a = area(i,1,f,e)/jacm1(j1,j2,1,e)
1854  unr(i,f,e) = a * ( rxm1(j1,j2,1,e)*unx(i,1,f,e)
1855  $ + rym1(j1,j2,1,e)*uny(i,1,f,e) )
1856  uns(i,f,e) = a * ( sxm1(j1,j2,1,e)*unx(i,1,f,e)
1857  $ + sym1(j1,j2,1,e)*uny(i,1,f,e) )
1858  unt(i,f,e) = 0.
1859  enddo
1860  enddo
1861  endif
1862  enddo
1863  enddo
1864 
1865  return
1866  end
1867 c-----------------------------------------------------------------------
unitvec
subroutine unitvec(X, Y, Z, N)
Definition: bdry.f:1831
rzero3
subroutine rzero3(A, B, C, N)
Definition: bdry.f:1821
exitt
void exitt()
Definition: comm_mpi.f:604
invmt
subroutine invmt(A, B, AA, N)
Definition: coef.f:1696
chkjac
subroutine chkjac(jac, n, iel, X, Y, Z, ND, IERR)
Definition: coef.f:965
genwz
subroutine genwz
Definition: coef.f:2
map31
subroutine map31(y, x, iel)
Definition: coef.f:1450
setarea
subroutine setarea
Definition: coef.f:1043
setinvm
subroutine setinvm
Definition: coef.f:1334
lubksb
subroutine lubksb(A, N, NP, INDX, B)
Definition: coef.f:1718
setwgtr
subroutine setwgtr(wgtr1, wgtr2, wgtr3, wgtr4)
Definition: coef.f:1135
map21e
subroutine map21e(y, x, iel)
Definition: coef.f:1609
glmapm1
subroutine glmapm1
Definition: coef.f:575
outxm3j
subroutine outxm3j(xm3, ym3, jm3)
Definition: coef.f:1676
area3
subroutine area3
Definition: coef.f:1188
out_fld_el
subroutine out_fld_el(x, e, c2)
Definition: coef.f:1659
glmapm3
subroutine glmapm3(xm3, ym3, zm3)
Definition: coef.f:391
xyzrst
subroutine xyzrst(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, XTM1, YTM1, ZTM1, IFAXIS)
Definition: coef.f:917
geom1
subroutine geom1(xm3, ym3, zm3)
Definition: coef.f:361
map21t
subroutine map21t(y, x, iel)
Definition: coef.f:1568
ludcmp
subroutine ludcmp(A, N, NP, INDX, D)
Definition: coef.f:1746
volume
subroutine volume
Definition: coef.f:1003
geom2
subroutine geom2
Definition: coef.f:823
lagmass
subroutine lagmass
Definition: coef.f:1315
set_unr
subroutine set_unr
Definition: coef.f:1811
area2
subroutine area2
Definition: coef.f:1066
map12
subroutine map12(y, x, iel)
Definition: coef.f:1530
geodat1
subroutine geodat1
Definition: coef.f:670
map13
subroutine map13(y, x, iel)
Definition: coef.f:1493
maprs
subroutine maprs(y, x, xa, nrest, iel)
Definition: coef.f:1381
out_xyz_el
subroutine out_xyz_el(x, y, z, e)
Definition: coef.f:1647
dsset
subroutine dsset(nx, ny, nz)
Definition: connect1.f:553
dssum
subroutine dssum(u, nx, ny, nz)
Definition: dssum.f:34
outmat
subroutine outmat(a, m, n, name6, ie)
Definition: fast3d.f:891
dot
real function dot(V1, V2, N)
Definition: genxyz.f:885
col3
subroutine col3(a, b, c, n)
Definition: math.f:598
ascol5
subroutine ascol5(a, b, c, d, e, n)
Definition: math.f:153
invers2
subroutine invers2(a, b, n)
Definition: math.f:51
addcol3
subroutine addcol3(a, b, c, n)
Definition: math.f:654
iglsum
function iglsum(a, n)
Definition: math.f:926
col2
subroutine col2(a, b, n)
Definition: math.f:564
copy
subroutine copy(a, b, n)
Definition: math.f:260
addcol4
subroutine addcol4(a, b, c, d, n)
Definition: math.f:142
vdot2
subroutine vdot2(dot, u1, u2, v1, v2, n)
Definition: math.f:447
vlsum
real function vlsum(vec, n)
Definition: math.f:417
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
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
invcol3
subroutine invcol3(a, b, c, n)
Definition: math.f:98
vdot3
subroutine vdot3(dot, u1, u2, u3, v1, v2, v3, n)
Definition: math.f:462
vcross
subroutine vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
Definition: math.f:430
rzero
subroutine rzero(a, n)
Definition: math.f:208
chsign
subroutine chsign(a, n)
Definition: math.f:305
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Definition: mxm_wrapper.f:2
outpost
subroutine outpost(v1, v2, v3, vp, vt, name3)
Definition: prepost.f:1378
zwglj
subroutine zwglj(Z, W, NP, ALPHA, BETA)
Definition: speclib.f:208
zwgl
subroutine zwgl(Z, W, NP)
Definition: speclib.f:91
dglj
subroutine dglj(D, DT, Z, NZ, lzd, ALPHA, BETA)
Definition: speclib.f:712
zwgll
subroutine zwgll(Z, W, NP)
Definition: speclib.f:108
iglm
subroutine iglm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2)
Definition: speclib.f:1077
igljm
subroutine igljm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2, ALPHA, BETA)
Definition: speclib.f:1153
dgllgl
subroutine dgllgl(D, DT, ZM1, ZM2, IM12, NZM1, NZM2, ND1, ND2)
Definition: speclib.f:935
igllm
subroutine igllm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2)
Definition: speclib.f:1102
igjm
subroutine igjm(I12, IT12, Z1, Z2, lz1, lz2, ND1, ND2, ALPHA, BETA)
Definition: speclib.f:1127
dgll
subroutine dgll(D, DT, Z, NZ, lzd)
Definition: speclib.f:801
zwgj
subroutine zwgj(Z, W, NP, ALPHA, BETA)
Definition: speclib.f:125
dgljgj
subroutine dgljgj(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
Definition: speclib.f:970
setaxdy
subroutine setaxdy(ifaxdy)
Definition: subs1.f:2342
setaxw2
subroutine setaxw2(IFAXWG)
Definition: subs2.f:17
facexv
subroutine facexv(A1, A2, A3, B1, B2, B3, IFACE1, IOP)
Definition: subs2.f:170
setaxw1
subroutine setaxw1(IFAXWG)
Definition: subs2.f:2