Nek5000
SEM for Incompressible NS
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
coef.F90
Go to the documentation of this file.
1 !-----------------------------------------------------------------
12 !-----------------------------------------------------------------
13 subroutine genwz
14  use size_m, only : ndim, nx1, ny1, nz1, nx2, ny2, nz2, nx3, ny3, nz3
15  use dxyz, only : dxm1, dxtm1, dym1, dytm1, dzm1, dztm1, dxm3, dxtm3
16  use dxyz, only : dym3, dytm3, dzm3, dztm3, dxm12, dxtm12, dym12, dytm12
17  use dxyz, only : dzm12, dztm12
18  use input, only : ifsplit
19  use ixyz, only : ixm12, ixtm12, iym12, iytm12, izm12, izm12
20  use ixyz, only : iztm12, ixm21, ixtm21, iym21, iytm21, izm21, iztm21
21  use ixyz, only : ixm13, ixtm13, iym13, iytm13, izm13, iztm13
22  use ixyz, only : ixm31, ixtm31, iym31, iytm31, izm31, iztm31
23  use wz_m, only : zgm1, zgm2, zgm3
24  use wz_m, only : wxm1, wym1, wzm1, wxm2, wym2, wzm2, wxm3, wym3, wzm3
25  use wz_m, only : w3m1, w3m2, w3m3
26  use speclib, only : zwgl, zwgll, dgll, zwglj, igllm, dgllgl
27  implicit none
28 
29  integer :: ix, iy, iz
30 
31  IF (ndim == 2) THEN
32  write (*,*) "Oops: ndim == 2"
33 #if 0
34  !*** Two-dimensional case **********************
35 
36 
37  ! Gauss-Lobatto Legendre mesh (suffix M1)
38  ! Generate collocation points and weights
39 
40  CALL zwgll(zgm1(1,1),wxm1,nx1)
41  CALL zwgll(zgm1(1,2),wym1,ny1)
42  zgm1(nz1,3) = 0.
43  wzm1(nz1) = 1.
44  DO 100 iy=1,ny1
45  DO 100 ix=1,nx1
46  w3m1(ix,iy,1)=wxm1(ix)*wym1(iy)
47  100 END DO
48 
49  ! Compute derivative matrices
50 
51  CALL dgll(dxm1,dxtm1,zgm1(1,1),nx1,nx1)
52  CALL dgll(dym1,dytm1,zgm1(1,2),ny1,ny1)
53  CALL rzero(dzm1 ,nz1*nz1)
54  CALL rzero(dztm1,nz1*nz1)
55 
56  ! Gauss Legendre mesh (suffix M2)
57  ! Generate collocation points and weights
58 
59  IF(ifsplit)THEN
60  CALL zwgll(zgm2(1,1),wxm2,nx2)
61  CALL zwgll(zgm2(1,2),wym2,ny2)
62  ELSE
63  CALL zwgl(zgm2(1,1),wxm2,nx2)
64  CALL zwgl(zgm2(1,2),wym2,ny2)
65  ENDIF
66  zgm2(nz2,3) = 0.
67  wzm2(nz2) = 1.
68  DO 200 iy=1,ny2
69  DO 200 ix=1,nx2
70  w3m2(ix,iy,1)=wxm2(ix)*wym2(iy)
71  200 END DO
72 
73  ! Gauss-Lobatto Legendre mesh (suffix M3).
74  ! Generate collocation points and weights.
75 
76  CALL zwgll(zgm3(1,1),wxm3,nx3)
77  CALL zwgll(zgm3(1,2),wym3,ny3)
78  zgm3(nz3,3) = 0.
79  wzm3(nz3) = 1.
80  DO 300 iy=1,ny3
81  DO 300 ix=1,nx3
82  w3m3(ix,iy,1)=wxm3(ix)*wym3(iy)
83  300 END DO
84 
85  ! Compute derivative matrices
86 
87  CALL dgll(dxm3,dxtm3,zgm3(1,1),nx3,nx3)
88  CALL dgll(dym3,dytm3,zgm3(1,2),ny3,ny3)
89  CALL rzero(dzm3 ,nz3*nz3)
90  CALL rzero(dztm3,nz3*nz3)
91 
92  ! Generate interpolation operators for the staggered mesh
93 
94  CALL igllm(ixm12,ixtm12,zgm1(1,1),zgm2(1,1),nx1,nx2,nx1,nx2)
95  CALL igllm(iym12,iytm12,zgm1(1,2),zgm2(1,2),ny1,ny2,ny1,ny2)
96  izm12(nz2,nz1) = 1.
97  iztm12(nz1,nz2) = 1.
98 
99  ! NOTE: The splitting scheme has only one mesh!!!!!
100 
101  IF (ifsplit) THEN
102  CALL igllm(ixm21,ixtm21,zgm1(1,1),zgm2(1,1),nx1,nx2,nx1,nx2)
103  CALL igllm(iym21,iytm21,zgm1(1,2),zgm2(1,2),ny1,ny2,ny1,ny2)
104  ELSE
105  write(*,*) "Oops: ifsplit false"
106 #if 0
107  CALL iglm(ixm21,ixtm21,zgm2(1,1),zgm1(1,1),nx2,nx1,nx2,nx1)
108  CALL iglm(iym21,iytm21,zgm2(1,2),zgm1(1,2),ny2,ny1,ny2,ny1)
109 #endif
110  ENDIF
111  izm21(nz1,nz2) = 1.
112  iztm21(nz2,nz1) = 1.
113 
114  ! Compute derivative operators for the staggered mesh
115 
116  IF(ifsplit)THEN
117  CALL copy(dxm12, dxm1, nx1*nx2)
118  CALL copy(dxtm12,dxtm1,nx1*nx2)
119  CALL copy(dym12, dym1, ny1*ny2)
120  CALL copy(dytm12,dytm1,ny1*ny2)
121  CALL copy(dzm12, dzm1, nz1*nz2)
122  CALL copy(dztm12,dztm1,nz1*nz2)
123  ELSE
124  CALL dgllgl(dxm12,dxtm12,zgm1(1,1),zgm2(1,1),ixm12, &
125  nx1,nx2,nx1,nx2)
126  CALL dgllgl(dym12,dytm12,zgm1(1,2),zgm2(1,2),iym12, &
127  ny1,ny2,ny1,ny2)
128  dzm12(nz2,nz1) = 0.
129  dztm12(nz2,nz1) = 0.
130  ENDIF
131 
132  ! Compute interpolation operators for the geometry mesh M3.
133 
134  CALL igllm(ixm13,ixtm13,zgm1(1,1),zgm3(1,1),nx1,nx3,nx1,nx3)
135  CALL igllm(iym13,iytm13,zgm1(1,2),zgm3(1,2),ny1,ny3,ny1,ny3)
136  CALL igllm(ixm31,ixtm31,zgm3(1,1),zgm1(1,1),nx3,nx1,nx3,nx1)
137  CALL igllm(iym31,iytm31,zgm3(1,2),zgm1(1,2),ny3,ny1,ny3,ny1)
138  izm13(nz3,nz1) = 1.
139  iztm13(nz1,nz3) = 1.
140  izm31(nz1,nz3) = 1.
141  iztm31(nz3,nz1) = 1.
142 
143 
144  IF (ifaxis) THEN
145 
146  ! Special treatment for the axisymmetric case
147  ! Generate additional points, weights, derivative operators and
148  ! interpolation operators required for elements close to the axis.
149 
150 
151  ! Gauss-Lobatto Jacobi mesh (suffix M1).
152  ! Generate collocation points and weights (alpha=0, beta=1).
153 
154  alpha = 0.
155  beta = 1.
156  CALL zwglj(zam1,wam1,ny1,alpha,beta)
157  DO 400 iy=1,ny1
158  DO 400 ix=1,nx1
159  w2am1(ix,iy)=wxm1(ix)*wam1(iy)
160  w2cm1(ix,iy)=wxm1(ix)*wym1(iy)
161  400 END DO
162 
163  ! Compute derivative matrices
164 
165  CALL copy(dcm1,dym1,ny1*ny1)
166  CALL copy(dctm1,dytm1,ny1*ny1)
167  CALL dglj(dam1,datm1,zam1,ny1,ny1,alpha,beta)
168 
169  ! Gauss Jacobi mesh (suffix M2)
170  ! Generate collocation points and weights
171 
172  IF(ifsplit)THEN
173  CALL zwglj(zam2,wam2,ny2,alpha,beta)
174  ELSE
175  CALL zwgj(zam2,wam2,ny2,alpha,beta)
176  ENDIF
177  DO 500 iy=1,ny2
178  DO 500 ix=1,nx2
179  w2cm2(ix,iy)=wxm2(ix)*wym2(iy)
180  w2am2(ix,iy)=wxm2(ix)*wam2(iy)
181  500 END DO
182 
183  ! Gauss-Lobatto Jacobi mesh (suffix M3).
184  ! Generate collocation points and weights.
185 
186  CALL zwglj(zam3,wam3,ny3,alpha,beta)
187  DO 600 iy=1,ny3
188  DO 600 ix=1,nx3
189  w2cm3(ix,iy)=wxm3(ix)*wym3(iy)
190  w2am3(ix,iy)=wxm3(ix)*wam3(iy)
191  600 END DO
192 
193  ! Compute derivative matrices
194 
195  CALL copy(dcm3,dym3,ny3*ny3)
196  CALL copy(dctm3,dytm3,ny3*ny3)
197  CALL dglj(dam3,datm3,zam3,ny3,ny3,alpha,beta)
198 
199  ! Generate interpolation operators for the staggered mesh
200 
201  CALL copy(icm12,iym12,ny2*ny1)
202  CALL copy(ictm12,iytm12,ny1*ny2)
203  CALL igljm(iam12,iatm12,zam1,zam2,ny1,ny2,ny1,ny2,alpha,beta)
204  CALL copy(icm21,iym21,ny1*ny2)
205  CALL copy(ictm21,iytm21,ny2*ny1)
206  IF (ifsplit) THEN
207  CALL igljm(iam21,iatm21,zam2,zam1,ny1,ny2,ny1,ny2,alpha,beta)
208  ELSE
209  CALL igjm(iam21,iatm21,zam2,zam1,ny2,ny1,ny2,ny1,alpha,beta)
210  ENDIF
211 
212  ! Compute derivative operators for the staggered mesh
213 
214  CALL copy(dcm12,dym12,ny2*ny1)
215  CALL copy(dctm12,dytm12,ny1*ny2)
216  IF(ifsplit)THEN
217  CALL copy(dam12, dam1, ny1*ny2)
218  CALL copy(datm12,datm1,ny1*ny2)
219  ELSE
220  CALL dgljgj(dam12,datm12,zam1,zam2,iam12, &
221  ny1,ny2,ny1,ny2,alpha,beta)
222  ENDIF
223 
224  ! Compute interpolation operators for the geometry mesh M3.
225 
226  CALL copy(icm13,iym13,ny3*ny1)
227  CALL copy(ictm13,iytm13,ny1*ny3)
228  CALL igljm(iam13,iatm13,zam1,zam3,ny1,ny3,ny1,ny3,alpha,beta)
229  CALL copy(icm31,iym31,ny1*ny3)
230  CALL copy(ictm31,iytm31,ny3*ny1)
231  CALL igljm(iam31,iatm31,zam3,zam1,ny3,ny1,ny3,ny1,alpha,beta)
232 
233  ! Compute interpolation operators between Gauss-Lobatto Jacobi
234  ! and Gauss-Lobatto Legendre (to be used in PREPOST).
235 
236  CALL igljm(iajl1,iatjl1,zam1,zgm1(1,2),ny1,ny1,ny1,ny1,alpha,beta)
237  IF (ifsplit) THEN
238  CALL igljm(iajl2,iatjl2,zam2,zgm2(1,2),ny2,ny2,ny2,ny2,alpha,beta)
239  ELSE
240  CALL igjm(iajl2,iatjl2,zam2,zgm2(1,2),ny2,ny2,ny2,ny2,alpha,beta)
241  ENDIF
242 
243  CALL invmt(iajl1 ,ialj1 ,tmp ,ny1)
244  CALL invmt(iatjl1,iatlj1,tmpt,ny1)
245  CALL mxm(iatjl1,ny1,iatlj1,ny1,tmpt,ny1)
246  CALL mxm(iajl1 ,ny1,ialj1 ,ny1,tmp ,ny1)
247 
248 
249  ! Compute interpolation operators between Gauss-Lobatto Legendre
250  ! and Gauss-Lobatto Jacobi (to be used in subr. genxyz IN postpre).
251 
252 
253  ! This call is not right, and these arrays are not used. 3/27/02. pff
254  ! CALL IGLLM(IALJ3,IATLJ3,ZGM3(1,2),ZAM3,NY3,NY3,NY3,NY3,ALPHA,BETA)
255  CALL igljm(ialj3,iatlj3,zgm3(1,2),zam3,ny3,ny3,ny3,ny3,alpha,beta)
256 
257  ENDIF
258 #endif
259 
260  ELSE
261 
262  !*** Three-dimensional case ************************************
263 
264 
265  ! Gauss-Lobatto Legendre mesh (suffix M1)
266  ! Generate collocation points and weights
267 
268  CALL zwgll(zgm1(1,1),wxm1,nx1)
269  CALL zwgll(zgm1(1,2),wym1,ny1)
270  CALL zwgll(zgm1(1,3),wzm1,nz1)
271  DO iz=1,nz1
272  DO iy=1,ny1
273  DO ix=1,nx1
274  w3m1(ix,iy,iz)=wxm1(ix)*wym1(iy)*wzm1(iz)
275  enddo
276  enddo
277  END DO
278 
279  ! Compute derivative matrices
280 
281  CALL dgll(dxm1,dxtm1,zgm1(1,1),nx1,nx1)
282  CALL dgll(dym1,dytm1,zgm1(1,2),ny1,ny1)
283  CALL dgll(dzm1,dztm1,zgm1(1,3),nz1,nz1)
284 
285  ! Gauss Legendre mesh (suffix M2)
286  ! Generate collocation points and weights
287 
288  IF(ifsplit)THEN
289  CALL zwgll(zgm2(1,1),wxm2,nx2)
290  CALL zwgll(zgm2(1,2),wym2,ny2)
291  CALL zwgll(zgm2(1,3),wzm2,nz2)
292  ELSE
293  CALL zwgl(zgm2(1,1),wxm2,nx2)
294  CALL zwgl(zgm2(1,2),wym2,ny2)
295  CALL zwgl(zgm2(1,3),wzm2,nz2)
296  ENDIF
297  DO iz=1,nz2
298  DO iy=1,ny2
299  DO ix=1,nx2
300  w3m2(ix,iy,iz)=wxm2(ix)*wym2(iy)*wzm2(iz)
301  enddo
302  enddo
303  END DO
304 
305  ! Gauss-Loabtto Legendre mesh (suffix M3).
306  ! Generate collocation points and weights.
307 
308  CALL zwgll(zgm3(1,1),wxm3,nx3)
309  CALL zwgll(zgm3(1,2),wym3,ny3)
310  CALL zwgll(zgm3(1,3),wzm3,nz3)
311  DO iz=1,nz3
312  DO iy=1,ny3
313  DO ix=1,nx3
314  w3m3(ix,iy,iz)=wxm3(ix)*wym3(iy)*wzm3(iz)
315  enddo
316  enddo
317  END DO
318 
319  ! Compute derivative matrices
320 
321  CALL dgll(dxm3,dxtm3,zgm3(1,1),nx3,nx3)
322  CALL dgll(dym3,dytm3,zgm3(1,2),ny3,ny3)
323  CALL dgll(dzm3,dztm3,zgm3(1,3),nz3,nz3)
324 
325  ! Generate interpolation operators for the staggered mesh
326 
327  CALL igllm(ixm12,ixtm12,zgm1(1,1),zgm2(1,1),nx1,nx2,nx1,nx2)
328  CALL igllm(iym12,iytm12,zgm1(1,2),zgm2(1,2),ny1,ny2,ny1,ny2)
329  CALL igllm(izm12,iztm12,zgm1(1,3),zgm2(1,3),nz1,nz2,nz1,nz2)
330 
331  ! NOTE: The splitting scheme has only one mesh!!!!!
332 
333  IF (ifsplit) THEN
334  CALL igllm(ixm21,ixtm21,zgm1(1,1),zgm2(1,1),nx1,nx2,nx1,nx2)
335  CALL igllm(iym21,iytm21,zgm1(1,2),zgm2(1,2),ny1,ny2,ny1,ny2)
336  CALL igllm(izm21,iztm21,zgm1(1,3),zgm2(1,3),nz1,nz2,nz1,nz2)
337  ELSE
338  write(*,*) "Oops: ifsplit is false"
339 #if 0
340  CALL iglm(ixm21,ixtm21,zgm2(1,1),zgm1(1,1),nx2,nx1,nx2,nx1)
341  CALL iglm(iym21,iytm21,zgm2(1,2),zgm1(1,2),ny2,ny1,ny2,ny1)
342  CALL iglm(izm21,iztm21,zgm2(1,3),zgm1(1,3),nz2,nz1,nz2,nz1)
343 #endif
344  ENDIF
345 
346  ! Compute derivative operators for the staggered mesh
347 
348  IF(ifsplit)THEN
349  CALL copy(dxm12, dxm1, nx1*nx2)
350  CALL copy(dxtm12,dxtm1,nx1*nx2)
351  CALL copy(dym12, dym1, ny1*ny2)
352  CALL copy(dytm12,dytm1,ny1*ny2)
353  CALL copy(dzm12, dzm1, nz1*nz2)
354  CALL copy(dztm12,dztm1,nz1*nz2)
355  ELSE
356  CALL dgllgl(dxm12,dxtm12,zgm1(1,1),zgm2(1,1),ixm12, &
357  nx1,nx2,nx1,nx2)
358  CALL dgllgl(dym12,dytm12,zgm1(1,2),zgm2(1,2),iym12, &
359  ny1,ny2,ny1,ny2)
360  CALL dgllgl(dzm12,dztm12,zgm1(1,3),zgm2(1,3),izm12, &
361  nz1,nz2,nz1,nz2)
362  ENDIF
363 
364  ! Compute interpolation operators for the geometry mesh M3.
365 
366  CALL igllm(ixm13,ixtm13,zgm1(1,1),zgm3(1,1),nx1,nx3,nx1,nx3)
367  CALL igllm(iym13,iytm13,zgm1(1,2),zgm3(1,2),ny1,ny3,ny1,ny3)
368  CALL igllm(izm13,iztm13,zgm1(1,3),zgm3(1,3),nz1,nz3,nz1,nz3)
369  CALL igllm(ixm31,ixtm31,zgm3(1,1),zgm1(1,1),nx3,nx1,nx3,nx1)
370  CALL igllm(iym31,iytm31,zgm3(1,2),zgm1(1,2),ny3,ny1,ny3,ny1)
371  CALL igllm(izm31,iztm31,zgm3(1,3),zgm1(1,3),nz3,nz1,nz3,nz1)
372 
373  ENDIF
374 
375  RETURN
376 end subroutine genwz
377 
378 !-----------------------------------------------------------------------
382 !-----------------------------------------------------------------------
383 subroutine geom1 ()!xm3,ym3,zm3)
384  use kinds, only : dp
385  use size_m, only : lx1, ly1, lz1, lelt
386  use geom, only : ifgmsh3
387  use tstep, only : istep
388  use mesh, only : if_ortho
389  implicit none
390 
391  real(DP), allocatable, dimension(:,:,:,:) :: &
392  XRM1, YRM1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1
393 
394  allocate(xrm1(lx1,ly1,lz1,lelt) &
395  , ysm1(lx1,ly1,lz1,lelt) &
396  , ztm1(lx1,ly1,lz1,lelt) )
397  if (if_ortho) then
398  allocate( &
399  yrm1(lx1,ly1,lz1,1) &
400  , xsm1(lx1,ly1,lz1,1) &
401  , xtm1(lx1,ly1,lz1,1) &
402  , ytm1(lx1,ly1,lz1,1) &
403  , zrm1(lx1,ly1,lz1,1) &
404  , zsm1(lx1,ly1,lz1,1) )
405  else
406  allocate( &
407  yrm1(lx1,ly1,lz1,lelt) &
408  , xsm1(lx1,ly1,lz1,lelt) &
409  , xtm1(lx1,ly1,lz1,lelt) &
410  , ytm1(lx1,ly1,lz1,lelt) &
411  , zrm1(lx1,ly1,lz1,lelt) &
412  , zsm1(lx1,ly1,lz1,lelt) )
413  endif
414 
415 
416  IF (ifgmsh3 .AND. istep == 0) THEN
417  write(*,*) "Oops: IFGMSH3"
418 !max CALL GLMAPM3 (XM3,YM3,ZM3)
419  ELSE
420  CALL glmapm1(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
421  ENDIF
422 
423  CALL geodat1(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
424 
425  RETURN
426 end subroutine geom1
427 
428 !-----------------------------------------------------------------------
439 !-----------------------------------------------------------------------
440 subroutine glmapm1(XRM1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
441  use kinds, only : dp
442  use size_m, only : ndim, nid
443  use size_m, only : nx1, ny1, nz1, nelt
444  use size_m, only : lx1, ly1, lz1, lelt
445  use geom, only : jacm1, jacmi
446  use geom, only : rxm1, rym1, rzm1, sxm1, sym1, szm1, txm1, tym1, tzm1
447  use geom, only : xm1, ym1, zm1
448  use input, only : ifaxis, ifxyo, ifvo, ifpo, ifto, param
449  use soln, only : vx, vy, vz, pr, t
450  use mesh, only : if_ortho
451  implicit none
452 
453 ! Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
454 ! share the same array structure in Scratch Common /SCRNS/.
455 ! real(DP) :: xrm1, yrm1, zrm1
456 ! real(DP) :: xsm1, ysm1
457 ! real(DP) :: xtm1, ytm1
458 
459  real(DP) :: XRM1(lx1,ly1,lz1,lelt) &
460  , YRM1(lx1,ly1,lz1,lelt) &
461  , XSM1(lx1,ly1,lz1,lelt) &
462  , YSM1(lx1,ly1,lz1,lelt) &
463  , XTM1(lx1,ly1,lz1,lelt) &
464  , YTM1(lx1,ly1,lz1,lelt) &
465  , ZRM1(lx1,ly1,lz1,lelt)
466  real(DP) :: ZSM1(lx1,ly1,lz1,lelt), ZTM1(lx1,ly1,lz1,lelt)
467 
468  integer :: nxy1, nyz1, nxyz1, ntot1
469  integer :: ie, ierr, kerr
470  integer, external :: iglsum
471 
472  nxy1 = nx1*ny1
473  nyz1 = ny1*nz1
474  nxyz1 = nx1*ny1*nz1
475  ntot1 = nxyz1*nelt
476 
477  CALL xyzrst(xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,xtm1,ytm1,ztm1, &
478  ifaxis)
479 
480  IF (ndim == 2) THEN
481 #if 0
482  CALL rzero(jacm1,ntot1)
483  CALL addcol3(jacm1,xrm1,ysm1,ntot1)
484  CALL subcol3(jacm1,xsm1,yrm1,ntot1)
485  CALL copy(rxm1,ysm1,ntot1)
486  CALL copy(rym1,xsm1,ntot1)
487  CALL chsign(rym1,ntot1)
488  CALL copy(sxm1,yrm1,ntot1)
489  CALL chsign(sxm1,ntot1)
490  CALL copy(sym1,xrm1,ntot1)
491  CALL rzero(rzm1,ntot1)
492  CALL rzero(szm1,ntot1)
493  CALL rone(tzm1,ntot1)
494 #endif
495  ELSE
496 
497  if (if_ortho) then
498  rxm1 = ysm1*ztm1
499  sym1 = xrm1*ztm1
500  tzm1 = xrm1*ysm1
501 
502  jacm1 = xrm1*ysm1*ztm1
503  else
504  rxm1 = ysm1*ztm1 - ytm1*zsm1
505  sym1 = xrm1*ztm1 - xtm1*zrm1
506  tzm1 = xrm1*ysm1 - xsm1*yrm1
507 
508  rym1 = xtm1*zsm1 - xsm1*ztm1
509  rzm1 = xsm1*ytm1 - xtm1*ysm1
510  sxm1 = ytm1*zrm1 - yrm1*ztm1
511  szm1 = xtm1*yrm1 - xrm1*ytm1
512  txm1 = yrm1*zsm1 - ysm1*zrm1
513  tym1 = xsm1*zrm1 - xrm1*zsm1
514 
515  jacm1 = xrm1*ysm1*ztm1 + xtm1*yrm1*zsm1 + xsm1*ytm1*zrm1 &
516  - xrm1*ytm1*zsm1 - xsm1*yrm1*ztm1 - xtm1*ysm1*zrm1
517  endif
518 
519  ENDIF
520 
521  kerr = 0
522  DO 500 ie=1,nelt
523  CALL chkjac(jacm1(1,1,1,ie),nxyz1,ie,xm1,ym1,zm1,ierr)
524  if (ierr /= 0) kerr = kerr+1
525  500 END DO
526  kerr = iglsum(kerr,1)
527  if (kerr > 0) then
528  ifxyo = .true.
529  ifvo = .false.
530  ifpo = .false.
531  ifto = .false.
532  param(66) = 4
533  call outpost(vx,vy,vz,pr,t,'xyz')
534  if (nid == 0) write(6,*) 'Jac error 1, setting p66=4, ifxyo=t'
535  call exitt
536  endif
537 
538  jacmi = 1_dp / jacm1
539 
540  RETURN
541 end subroutine glmapm1
542 
543 !-----------------------------------------------------------------------
546 !-----------------------------------------------------------------------
547 subroutine geodat1(XRM1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
548  use kinds, only : dp
549  use size_m, only : lx1, ly1, lz1, lelt
550  use size_m, only : nx1, ny1, nz1, nelt, ndim
551  use geom, only : ifgmsh3, jacm1, ifrzer, ym1, bm1
552  use geom, only : g1m1, rxm1, rym1, rzm1, g2m1, sxm1, sym1, szm1
553  use geom, only : g3m1, txm1, tym1, tzm1, g4m1, g5m1, g6m1
554  use input, only : ifaxis
555  use tstep, only : istep
556  use wz_m, only : zam1, w3m1
557  use mesh, only : if_ortho
558  implicit none
559 
560 ! Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
561 ! share the same array structure in Scratch Common /SCRNS/.
562 
563  real(DP) :: XRM1(lx1,ly1,lz1,lelt) &
564  , YRM1(lx1,ly1,lz1,lelt) &
565  , XSM1(lx1,ly1,lz1,lelt) &
566  , YSM1(lx1,ly1,lz1,lelt) &
567  , XTM1(lx1,ly1,lz1,lelt) &
568  , YTM1(lx1,ly1,lz1,lelt) &
569  , ZRM1(lx1,ly1,lz1,lelt)
570  real(DP) :: ZSM1(lx1,ly1,lz1,lelt) &
571  , ZTM1(lx1,ly1,lz1,lelt)
572  real(DP), allocatable :: WJ(:,:,:,:)
573 
574  integer :: nxyz1, ntot1, iel, j, i
575 
576  nxyz1 = nx1*ny1*nz1
577  ntot1 = nxyz1*nelt
578 
579 
580  IF (ifgmsh3 .AND. istep == 0) then
581  write(*,*) "Oops: IFGMSH3"
582 !max CALL XYZRST (XRM1,YRM1,ZRM1,XSM1,YSM1,ZSM1,XTM1,YTM1,ZTM1, &
583 ! IFAXIS)
584  endif
585 
586  allocate(wj(lx1,ly1,lz1,lelt))
587  IF ( .NOT. ifaxis) THEN
588  wj = 1_dp / jacm1
589  ELSE
590  DO 500 iel=1,nelt
591  IF (ifrzer(iel)) THEN
592  DO j=1,ny1
593  DO i=1,nx1
594  IF (j > 1) THEN
595  wj(i,j,1,iel) = ym1(i,j,1,iel)/ &
596  (jacm1(i,j,1,iel)*(1.+zam1(j)))
597  ELSE
598  wj(i,j,1,iel) = ysm1(i,j,1,iel)/jacm1(i,j,1,iel)
599  ENDIF
600  enddo
601  END DO
602  ELSE
603  wj(:,:,:,iel) = ym1(:,:,:,iel) / jacm1(:,:,:,iel)
604  ENDIF
605  500 END DO
606  ENDIF
607 
608 ! Compute geometric factors for integrated del-squared operator.
609 
610  IF (ndim == 2) THEN
611  write(*,*) "Whoops! geodat1"
612 #if 0
613  CALL vdot2(g1m1,rxm1,rym1,rxm1,rym1,ntot1)
614  CALL vdot2(g2m1,sxm1,sym1,sxm1,sym1,ntot1)
615  CALL vdot2(g4m1,rxm1,rym1,sxm1,sym1,ntot1)
616  CALL col2(g1m1,wj,ntot1)
617  CALL col2(g2m1,wj,ntot1)
618  CALL col2(g4m1,wj,ntot1)
619  CALL rzero(g3m1,ntot1)
620  CALL rzero(g5m1,ntot1)
621  CALL rzero(g6m1,ntot1)
622 #endif
623  ELSE
624  if (if_ortho) then
625  g1m1 = wj * (rxm1 * rxm1)
626  g2m1 = wj * (sym1 * sym1)
627  g3m1 = wj * (tzm1 * tzm1)
628 
629  else
630  g1m1 = wj * (rxm1 * rxm1 + rym1 * rym1 + rzm1 * rzm1)
631  g2m1 = wj * (sxm1 * sxm1 + sym1 * sym1 + szm1 * szm1)
632  g3m1 = wj * (txm1 * txm1 + tym1 * tym1 + tzm1 * tzm1)
633 
634  g4m1 = wj * (rxm1 * sxm1 + rym1 * sym1 + rzm1 * szm1)
635  g5m1 = wj * (rxm1 * txm1 + rym1 * tym1 + rzm1 * tzm1)
636  g6m1 = wj * (txm1 * sxm1 + tym1 * sym1 + tzm1 * szm1)
637  endif
638  ENDIF
639  deallocate(wj)
640 
641 ! Multiply the geometric factors GiM1,i=1,5 with the
642 ! weights on mesh M1.
643 
644  DO iel=1,nelt
645 !max IF (IFAXIS) CALL SETAXW1 ( IFRZER(IEL) )
646  g1m1(:,:,:,iel) = g1m1(:,:,:,iel) * w3m1
647  g2m1(:,:,:,iel) = g2m1(:,:,:,iel) * w3m1
648  ! CALL COL2 (G4M1(1,1,1,IEL),W3M1,NXYZ1)
649  IF (ndim == 3) THEN
650  g3m1(:,:,:,iel) = g3m1(:,:,:,iel) * w3m1
651  ! CALL COL2 (G5M1(1,1,1,IEL),W3M1,NXYZ1)
652  ! CALL COL2 (G6M1(1,1,1,IEL),W3M1,NXYZ1)
653  ENDIF
654  END DO
655 
656 ! Compute the mass matrix on mesh M1.
657 
658  DO iel=1,nelt
659 !max IF (IFAXIS) CALL SETAXW1 ( IFRZER(IEL) )
660  bm1(:,:,:,iel) = jacm1(:,:,:,iel) * w3m1
661  IF (ifaxis) THEN
662  write(*,*) "Oops: ifaxis"
663 #if 0
664  CALL col3(baxm1(1,1,1,iel),jacm1(1,1,1,iel),w3m1,nxyz1)
665  IF (ifrzer(iel)) THEN
666  DO 600 j=1,ny1
667  IF (j > 1) THEN
668  DO 610 i=1,nx1
669  bm1(i,j,1,iel) = bm1(i,j,1,iel)*ym1(i,j,1,iel) &
670  /(1.+zam1(j))
671  baxm1(i,j,1,iel)=baxm1(i,j,1,iel)/(1.+zam1(j))
672  610 END DO
673  ELSE
674  DO 620 i=1,nx1
675  bm1(i,j,1,iel) = bm1(i,j,1,iel)*ysm1(i,j,1,iel)
676  baxm1(i,j,1,iel)=baxm1(i,j,1,iel)
677  620 END DO
678  ENDIF
679  600 END DO
680  ELSE
681  CALL col2(bm1(1,1,1,iel),ym1(1,1,1,iel),nxyz1)
682  ENDIF
683 #endif
684  ENDIF
685 
686  END DO
687 
688  IF(ifaxis) THEN
689  write(*,*) "Oops: ifaxis"
690 #if 0
691  DO iel=1,nelt
692  IF(ifrzer(iel)) THEN
693  DO j=1,ny1
694  DO i=1,nx1
695  IF(j == 1) THEN
696  yinvm1(i,j,1,iel)=1.0d0/ysm1(i,j,1,iel)
697  ELSE
698  yinvm1(i,j,1,iel)=1.0d0/ym1(i,j,1,iel)
699  ENDIF
700  ENDDO
701  ENDDO
702  ELSE
703  yinvm1(:,:,:,iel) = 1_dp / ym(:,:,:,iel)
704  ENDIF
705  ENDDO
706 #endif
707  ENDIF
708 
709 ! Compute normals, tangents, and areas on elemental surfaces
710 
711  CALL setarea(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
712 
713  RETURN
714 end subroutine geodat1
715 
716 !-------------------------------------------------------------------
724 !------------------------------------------------------------------
725 subroutine geom2
726  use size_m, only : nx2, ny2, nz2, nelv
727  use geom, only : rxm2, rxm1, rym2, rym1, rzm2, rzm1
728  use geom, only : sxm2, sxm1, sym2, sym1, szm2, szm1
729  use geom, only : txm2, txm1, tym2, tym1, tzm2, tzm1
730  use geom, only : jacm2, jacm1, xm2, xm1, ym2, ym1, zm2, zm1
731  use input, only : ifsplit
732  use geom, only : bm2, bm1
733  implicit none
734 
735  integer :: nxyz2, ntot2
736 
737  nxyz2 = nx2*ny2*nz2
738  ntot2 = nxyz2*nelv
739 
740  IF (ifsplit) THEN
741  ! Mesh 1 and 2 are identical
742  rxm2 => rxm1
743  rym2 => rym1
744  rzm2 => rzm1
745 
746  sxm2 => sxm1
747  sym2 => sym1
748  szm2 => szm1
749 
750  txm2 => txm1
751  tym2 => tym1
752  tzm2 => tzm1
753 
754  jacm2 => jacm1
755  bm2 => bm1
756 
757  xm2 => xm1
758  ym2 => ym1
759  zm2 => zm1
760 
761  ELSE
762  write(*,*) "Oops: ifsplit"
763 #if 0
764  ! Consistent approximation spaces (UZAWA)
765 
766  IF (ndim == 2) THEN
767  CALL rzero(rzm2,ntot2)
768  CALL rzero(szm2,ntot2)
769  CALL rone(tzm2,ntot2)
770  ENDIF
771 
772  DO 1000 iel=1,nelv
773 
774  ! Mapping from mesh M1 to mesh M2
775 
776  CALL map12(rxm2(1,1,1,iel),rxm1(1,1,1,iel),iel)
777  CALL map12(rym2(1,1,1,iel),rym1(1,1,1,iel),iel)
778  CALL map12(sxm2(1,1,1,iel),sxm1(1,1,1,iel),iel)
779  CALL map12(sym2(1,1,1,iel),sym1(1,1,1,iel),iel)
780  IF (ndim == 3) THEN
781  CALL map12(rzm2(1,1,1,iel),rzm1(1,1,1,iel),iel)
782  CALL map12(szm2(1,1,1,iel),szm1(1,1,1,iel),iel)
783  CALL map12(txm2(1,1,1,iel),txm1(1,1,1,iel),iel)
784  CALL map12(tym2(1,1,1,iel),tym1(1,1,1,iel),iel)
785  CALL map12(tzm2(1,1,1,iel),tzm1(1,1,1,iel),iel)
786  ENDIF
787  CALL map12(jacm2(1,1,1,iel),jacm1(1,1,1,iel),iel)
788 
789  CALL map12(xm2(1,1,1,iel),xm1(1,1,1,iel),iel)
790  CALL map12(ym2(1,1,1,iel),ym1(1,1,1,iel),iel)
791  CALL map12(zm2(1,1,1,iel),zm1(1,1,1,iel),iel)
792 
793  ! Compute the mass matrix on mesh M2.
794 
795  IF (ifaxis) CALL setaxw2( ifrzer(iel) )
796  CALL col3(bm2(1,1,1,iel),w3m2,jacm2(1,1,1,iel),nxyz2)
797 
798  IF (ifaxis .AND. ifrzer(iel)) THEN
799  DO 300 j=1,ny2
800  DO 300 i=1,nx2
801  bm2(i,j,1,iel) = bm2(i,j,1,iel)*ym2(i,j,1,iel) &
802  /(1.+zam2(j))
803  300 END DO
804  ELSEIF (ifaxis .AND. ( .NOT. ifrzer(iel))) THEN
805  CALL col2(bm2(1,1,1,iel),ym2(1,1,1,iel),nxyz2)
806  ENDIF
807  1000 END DO
808 #endif
809  ENDIF
810 
811 ! Compute inverse of mesh 2 mass matrix, pff 3/5/92
812 !max CALL INVERS2(BM2INV,BM2,NTOT2)
813 
814  RETURN
815 end subroutine geom2
816 
817 !-----------------------------------------------------------------------
819 !-----------------------------------------------------------------------
820 subroutine xyzrst (xrm1,yrm1,zrm1,xsm1,ysm1,zsm1, XTM1,YTM1,ZTM1,IFAXIS)
821  use kinds, only : dp
822  use size_m, only : lx1, ly1, lz1, nx1, ny1, nz1, nelt, ndim
823  use dxyz, only : dxm1, dytm1, dztm1
824  use geom, only : xm1, ym1, zm1
825  use mesh, only : if_ortho
826  implicit none
827 
828  real(DP) :: XRM1(lx1,ly1,lz1,*), YRM1(lx1,ly1,lz1,*), ZRM1(lx1,ly1,lz1,*)
829  real(DP) :: XSM1(lx1,ly1,lz1,*), YSM1(lx1,ly1,lz1,*), ZSM1(lx1,ly1,lz1,*)
830  real(DP) :: XTM1(lx1,ly1,lz1,*), YTM1(lx1,ly1,lz1,*), ZTM1(lx1,ly1,lz1,*)
831  LOGICAL :: IFAXIS
832 
833  integer :: nxy1, nyz1, iel, iz
834 
835  nxy1=nx1*ny1
836  nyz1=ny1*nz1
837 
839  DO iel=1,nelt
840 
841  IF (ifaxis) then
842  write(*,*) "Oops: ifaxis"
843  !CALL SETAXDY ( IFRZER(IEL) )
844  endif
845 
846  CALL mxm(dxm1,nx1,xm1(1,1,1,iel),nx1,xrm1(1,1,1,iel),nyz1)
847  if (.not. if_ortho) CALL mxm(dxm1,nx1,ym1(1,1,1,iel),nx1,yrm1(1,1,1,iel),nyz1)
848  if (.not. if_ortho) CALL mxm(dxm1,nx1,zm1(1,1,1,iel),nx1,zrm1(1,1,1,iel),nyz1)
849 
850  DO iz=1,nz1
851  if (.not. if_ortho) CALL mxm(xm1(1,1,iz,iel),nx1,dytm1,ny1,xsm1(1,1,iz,iel),ny1)
852  CALL mxm(ym1(1,1,iz,iel),nx1,dytm1,ny1,ysm1(1,1,iz,iel),ny1)
853  if (.not. if_ortho) CALL mxm(zm1(1,1,iz,iel),nx1,dytm1,ny1,zsm1(1,1,iz,iel),ny1)
854  END DO
855 
856  IF (ndim == 3) THEN
857  if (.not. if_ortho) CALL mxm(xm1(1,1,1,iel),nxy1,dztm1,nz1,xtm1(1,1,1,iel),nz1)
858  if (.not. if_ortho) CALL mxm(ym1(1,1,1,iel),nxy1,dztm1,nz1,ytm1(1,1,1,iel),nz1)
859  CALL mxm(zm1(1,1,1,iel),nxy1,dztm1,nz1,ztm1(1,1,1,iel),nz1)
860  ELSE
861  if (.not. if_ortho) xtm1(:,:,:,iel) = 0._dp
862  if (.not. if_ortho) ytm1(:,:,:,iel) = 0._dp
863  ztm1(:,:,:,iel) = 1._dp
864  ENDIF
865 
866  END DO
867 
868  RETURN
869 end subroutine xyzrst
870 
872 subroutine chkjac(jac,n,iel,X,Y,Z,IERR)
873  use kinds, only : dp
874  use size_m
875  use parallel
876  implicit none
877 
878  integer :: n, iel, ierr
879  REAL(DP) :: JAC(n),x(1),y(1),z(1)
880  real(DP) :: sign
881  integer :: i, ieg
882 
883  ierr = 1
884  sign = jac(1)
885  DO 100 i=2,n
886  IF (sign*jac(i) <= 0._dp) THEN
887  ieg = lglel(iel)
888  WRITE(6,101) nid,i,ieg
889  write(6,*) jac(i-1),jac(i)
890  if (ndim == 3) then
891  write(6,7) nid,x(i-1),y(i-1),z(i-1)
892  write(6,7) nid,x(i),y(i),z(i)
893  else
894  write(6,7) nid,x(i-1),y(i-1)
895  write(6,7) nid,x(i),y(i)
896  endif
897  7 format(i5,' xyz:',1p3e14.5)
898  ! if (np.eq.1) call out_xyz_el(x,y,z,iel)
899  ! ierr=0
900  return
901  ENDIF
902  100 END DO
903  101 FORMAT(//,i5,2x &
904  ,'ERROR: Vanishing Jacobian near',i7,'th node of element' &
905  ,i10,'.')
906 
907 
908  ierr = 0
909  RETURN
910 end subroutine chkjac
911 
912 !-----------------------------------------------------------------------
914 subroutine volume
915  use kinds, only : dp
916  use size_m, only : nx1, ny1, nz1, nelv, nelt, nx2, ny2, nz2, nfield, nid
917  use esolv, only : volel
918  use input, only : ifmvbd, ifmhd, iftmsh
919  use geom, only : volvm1, volvm2, voltm1, voltm2, bm1, bm2
920  use tstep, only : volfld
921  implicit none
922 
923  real(DP), external :: glsum
924  integer :: e, mfield, nfldt, ifld, nxyz
925 
926  volvm1=glsum(bm1,nx1*ny1*nz1*nelv)
927  volvm2=glsum(bm2,nx2*ny2*nz2*nelv)
928  voltm1=glsum(bm1,nx1*ny1*nz1*nelt)
929  voltm2=glsum(bm2,nx2*ny2*nz2*nelt)
930  mfield=1
931  if (ifmvbd) mfield=0
932  nfldt = nfield
933  if (ifmhd) nfldt = nfield+1
934 
935  do ifld=mfield,nfldt
936  if (iftmsh(ifld)) then
937  volfld(ifld) = voltm1
938  else
939  volfld(ifld) = volvm1
940  endif
941  enddo
942 
943  if (nid == 0) write(6,*) 'vol_t,vol_v:',voltm1,volvm1
944 
945 
946  nxyz = nx1*ny1*nz1
947  do e=1,nelt
948  volel(e) = sum(bm1(:,:,:,e))
949  enddo
950 
951  return
952 end subroutine volume
953 
954 !-----------------------------------------------------------------------
956 subroutine setarea(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
957  use kinds, only : dp
958  use size_m, only : lx1, ly1, lz1
959  use size_m, only : nx1, nz1, nelt, ndim
960  use geom, only : area, unx, uny, unz
961  implicit none
962 
963  real(DP) :: XRM1(lx1,ly1,lz1,*) &
964  , YRM1(lx1,ly1,lz1,*) &
965  , XSM1(lx1,ly1,lz1,*) &
966  , YSM1(lx1,ly1,lz1,*) &
967  , XTM1(lx1,ly1,lz1,*) &
968  , YTM1(lx1,ly1,lz1,*) &
969  , ZRM1(lx1,ly1,lz1,*) &
970  , ZSM1(lx1,ly1,lz1,*) &
971  , ZTM1(lx1,ly1,lz1,*)
972 
973  integer :: nsrf
974 
975  nsrf = 6*nx1*nz1*nelt
976 
977  area = 0_dp
978  unx = 0_dp; uny = 0_dp; unz = 0_dp
979 ! CALL RZERO3 (T1X,T1Y,T1Z,NSRF)
980 ! CALL RZERO3 (T2X,T2Y,T2Z,NSRF)
981 
982  IF (ndim == 2) THEN
983 !max CALL AREA2
984  ELSE
985  CALL area3(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
986  ENDIF
987 
988  RETURN
989 end subroutine setarea
990 
991 !--------------------------------------------------------------------
993 !--------------------------------------------------------------------
994 subroutine area3(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
995  use kinds, only : dp
996  use size_m, only : lx1, ly1, lz1, lelt
997  use size_m, only : nx1, ny1, nz1, nelt, ndim
998  use geom, only : area, unx, uny, unz
999  use wz_m, only : wxm1, wym1, wzm1
1000  use mesh, only : if_ortho
1001  implicit none
1002 
1003 ! Note: Subroutines GLMAPM1, GEODAT1, AREA2, SETWGTR and AREA3
1004 ! share the same array structure in Scratch Common /SCRNS/.
1005 
1006  real(DP) :: XRM1(lx1,ly1,lz1,*) &
1007  , YRM1(lx1,ly1,lz1,*) &
1008  , XSM1(lx1,ly1,lz1,*) &
1009  , YSM1(lx1,ly1,lz1,*) &
1010  , XTM1(lx1,ly1,lz1,*) &
1011  , YTM1(lx1,ly1,lz1,*) &
1012  , ZRM1(lx1,ly1,lz1,*)
1013  real(DP) :: ZSM1(lx1,ly1,lz1,*) &
1014  , ZTM1(lx1,ly1,lz1,*)
1015 
1016  real(DP), allocatable :: A(:,:,:,:), B(:,:,:,:), C(:,:,:,:), dot(:,:,:,:)
1017  integer :: nxy1, nface, ntot, nsrf, iel, iz, iy, ix
1018  real(DP) :: weight
1019 
1020  nxy1 = nx1*ny1
1021  nface = 2*ndim
1022  ntot = nx1*ny1*nz1*nelt
1023  nsrf = 6*nx1*ny1*nelt
1024 
1025 ! "R"
1026  allocate(a(lx1,ly1,lz1,lelt), b(lx1,ly1,lz1,lelt), c(lx1,ly1,lz1,lelt))
1027  allocate(dot(lx1,ly1,lz1,lelt))
1028  if (if_ortho) then
1029  a = ysm1(:,:,:,1:nelt) * ztm1(:,:,:,1:nelt)
1030  b = 0._dp; c = 0._dp
1031  else
1032  CALL vcross(a,b,c,xsm1,ysm1,zsm1,xtm1,ytm1,ztm1,ntot)
1033  endif
1034  dot = a*a + b*b + c*c
1035 
1036  DO iel=1,nelt
1037  DO iz=1,nz1
1038  DO iy=1,ny1
1039  weight = wym1(iy)*wzm1(iz)
1040  area(iy,iz,2,iel) = sqrt(dot(nx1,iy,iz,iel))*weight
1041  area(iy,iz,4,iel) = sqrt(dot( 1,iy,iz,iel))*weight
1042  unx(iy,iz,4,iel) = -a( 1,iy,iz,iel)
1043  unx(iy,iz,2,iel) = a(nx1,iy,iz,iel)
1044  uny(iy,iz,4,iel) = -b( 1,iy,iz,iel)
1045  uny(iy,iz,2,iel) = b(nx1,iy,iz,iel)
1046  unz(iy,iz,4,iel) = -c( 1,iy,iz,iel)
1047  unz(iy,iz,2,iel) = c(nx1,iy,iz,iel)
1048  enddo
1049  enddo
1050  END DO
1051 
1052 ! "S"
1053  if (if_ortho) then
1054  b = -xrm1(:,:,:,1:nelt) * ztm1(:,:,:,1:nelt)
1055  a = 0._dp; c = 0._dp
1056  else
1057  CALL vcross(a,b,c,xrm1,yrm1,zrm1,xtm1,ytm1,ztm1,ntot)
1058  endif
1059  dot = a*a + b*b + c*c
1060  DO iel=1,nelt
1061  DO iz=1,nz1
1062  DO ix=1,nx1
1063  weight=wxm1(ix)*wzm1(iz)
1064  area(ix,iz,1,iel) = sqrt(dot(ix, 1,iz,iel))*weight
1065  area(ix,iz,3,iel) = sqrt(dot(ix,ny1,iz,iel))*weight
1066  unx(ix,iz,1,iel) = a(ix, 1,iz,iel)
1067  unx(ix,iz,3,iel) = -a(ix,ny1,iz,iel)
1068  uny(ix,iz,1,iel) = b(ix, 1,iz,iel)
1069  uny(ix,iz,3,iel) = -b(ix,ny1,iz,iel)
1070  unz(ix,iz,1,iel) = c(ix, 1,iz,iel)
1071  unz(ix,iz,3,iel) = -c(ix,ny1,iz,iel)
1072  enddo
1073  enddo
1074  END DO
1075 
1076 ! "T"
1077  if (if_ortho) then
1078  c = xrm1(:,:,:,1:nelt) * ysm1(:,:,:,1:nelt)
1079  a = 0._dp; b = 0._dp
1080  else
1081  CALL vcross(a,b,c,xrm1,yrm1,zrm1,xsm1,ysm1,zsm1,ntot)
1082  endif
1083  dot = a*a + b*b + c*c
1084  DO iel=1,nelt
1085  DO ix=1,nx1
1086  DO iy=1,ny1
1087  weight=wxm1(ix)*wym1(iy)
1088  area(ix,iy,5,iel) = sqrt(dot(ix,iy, 1,iel))*weight
1089  area(ix,iy,6,iel) = sqrt(dot(ix,iy,nz1,iel))*weight
1090  unx(ix,iy,5,iel) = -a(ix,iy, 1,iel)
1091  unx(ix,iy,6,iel) = a(ix,iy,nz1,iel)
1092  uny(ix,iy,5,iel) = -b(ix,iy, 1,iel)
1093  uny(ix,iy,6,iel) = b(ix,iy,nz1,iel)
1094  unz(ix,iy,5,iel) = -c(ix,iy, 1,iel)
1095  unz(ix,iy,6,iel) = c(ix,iy,nz1,iel)
1096  enddo
1097  enddo
1098  END DO
1099 
1100  CALL unitvec(unx,uny,unz,nsrf)
1101 
1102 ! COMPUTE UNIT TANGENT T1
1103 #if 0
1104  DO 600 iel=1,nelt
1105  DO 600 ifc=1,nface
1106  IF (ifc == 1 .OR. ifc == 6) THEN
1107  CALL facexv(t1x(1,1,ifc,iel),t1y(1,1,ifc,iel), &
1108  t1z(1,1,ifc,iel), &
1109  xrm1(1,1,1,iel),yrm1(1,1,1,iel), &
1110  zrm1(1,1,1,iel),ifc,0)
1111  ELSEIF (ifc == 2 .OR. ifc == 5) THEN
1112  CALL facexv(t1x(1,1,ifc,iel),t1y(1,1,ifc,iel), &
1113  t1z(1,1,ifc,iel), &
1114  xsm1(1,1,1,iel),ysm1(1,1,1,iel), &
1115  zsm1(1,1,1,iel),ifc,0)
1116  ELSE
1117  CALL facexv(t1x(1,1,ifc,iel),t1y(1,1,ifc,iel), &
1118  t1z(1,1,ifc,iel), &
1119  xtm1(1,1,1,iel),ytm1(1,1,1,iel), &
1120  ztm1(1,1,1,iel),ifc,0)
1121  ENDIF
1122  600 END DO
1123 
1124  CALL unitvec(t1x,t1y,t1z,nsrf)
1125 
1126 ! COMPUTE UNIT TANGENT T2 ( T2 = Normal X T1 )
1127  DO 700 iel=1,nelt
1128  DO 700 ifc=1,nface
1129  CALL vcross(t2x(1,1,ifc,iel),t2y(1,1,ifc,iel), &
1130  t2z(1,1,ifc,iel), &
1131  unx(1,1,ifc,iel),uny(1,1,ifc,iel), &
1132  unz(1,1,ifc,iel), &
1133  t1x(1,1,ifc,iel),t1y(1,1,ifc,iel), &
1134  t1z(1,1,ifc,iel),nxy1)
1135  700 END DO
1136 #endif
1137 
1138  RETURN
1139 end subroutine area3
1140 
1141 !--------------------------------------------------------------------
1143 !--------------------------------------------------------------------
1144 subroutine lagmass()
1145  use size_m, only : nx1, ny1, nz1, nelt
1146  use geom, only : bm1, bm1lag
1147  use tstep, only : nbdinp
1148  implicit none
1149 
1150  integer :: ntot1, ilag
1151  ntot1 = nx1*ny1*nz1*nelt
1152  DO ilag=nbdinp-1,2,-1
1153  CALL copy(bm1lag(1,1,1,1,ilag),bm1lag(1,1,1,1,ilag-1),ntot1)
1154  END DO
1155  CALL copy(bm1lag(1,1,1,1,1),bm1,ntot1)
1156 
1157  RETURN
1158 end subroutine lagmass
1159 
1160 !--------------------------------------------------------------------
1167 !--------------------------------------------------------------------
1168 subroutine setinvm()
1169  use kinds, only : dp
1170  use size_m, only : nx1, ny1, nz1, nelv, nelt, nfield
1171  use input, only : ifflow, ifheat, iftmsh
1172  use geom, only : bm1, binvm1, bintm1
1173  use tstep, only : ifield
1174  implicit none
1175 
1176  logical :: any_iftmsh
1177  integer :: nxyz1, ifld, store_field, ntot
1178  nxyz1 = nx1*ny1*nz1
1179 
1180  store_field = ifield
1181 
1182  IF (ifflow) THEN ! Velocity mass matrix
1183  ifield = 1
1184  ntot = nxyz1*nelv
1185  CALL copy(binvm1,bm1,ntot)
1186  CALL dssum(binvm1)
1187  binvm1 = 1._dp / binvm1
1188  ENDIF
1189 
1190  any_iftmsh = .false.
1191  do ifld = 1,nfield
1192  if (iftmsh(ifld)) any_iftmsh = .true.
1193  enddo
1194 
1195  IF (ifheat .and. any_iftmsh) THEN ! Temperature mass matrix
1196  ifield = 2
1197  ntot = nxyz1*nelt
1198  CALL copy(bintm1,bm1,ntot)
1199  CALL dssum(bintm1)
1200  bintm1 = 1._dp / bintm1
1201  ENDIF
1202 
1203  ifield = store_field
1204 
1205  return
1206 end subroutine setinvm
1207 !-----------------------------------------------------------------------
sum
static double sum(struct xxt *data, double v, uint n, uint tag)
Definition: xxt.c:400
dssum
subroutine dssum(u)
Direct stiffness sum.
Definition: dssum.F90:54
tstep
cleaned
Definition: tstep_mod.F90:2
input
Input parameters from preprocessors.
Definition: input_mod.F90:11
soln
Definition: soln_mod.F90:1
speclib::igjm
subroutine igjm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
Definition: speclib.F90:1227
wz_m
cleaned
Definition: wz_mod.F90:2
mxm
subroutine mxm(a, n1, b, n2, c, n3)
Compute matrix-matrix product C = A*B.
Definition: mxm_wrapper.F90:7
speclib::dgllgl
subroutine, public dgllgl(D, DT, ZM1, ZM2, IM12, NZM1, NZM2, ND1, ND2)
Definition: speclib.F90:1015
speclib::igljm
subroutine igljm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2, ALPHA, BETA)
Definition: speclib.F90:1257
dxyz
Definition: dxyz_mod.F90:3
glmapm1
subroutine glmapm1(XRM1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
Routine to generate mapping data based on mesh 1 (Gauss-Legendre Lobatto meshes). ...
Definition: coef.F90:440
speclib::dgljgj
subroutine dgljgj(D, DT, ZGL, ZG, IGLG, NPGL, NPG, ND1, ND2, ALPHA, BETA)
Definition: speclib.F90:1055
vcross
subroutine vcross(u1, u2, u3, v1, v2, v3, w1, w2, w3, n)
Compute a Cartesian vector cross product.
Definition: math.F90:95
speclib::iglm
subroutine iglm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
Definition: speclib.F90:1167
genwz
subroutine genwz
Generate derivative and interpolation operators. GENERATE.
Definition: coef.F90:13
exitt
void exitt()
Definition: comm_mpi.F90:411
xyzrst
subroutine xyzrst(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, XTM1, YTM1, ZTM1, IFAXIS)
Compute global-to-local derivatives on mesh 1.
Definition: coef.F90:820
area3
subroutine area3(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
Compute areas, normals and tangents (3D geom.)
Definition: coef.F90:994
mesh
cleaned
Definition: mesh_mod.F90:2
geodat1
subroutine geodat1(XRM1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
Routine to generate elemental geometric matrices on mesh 1 (Gauss-Legendre Lobatto mesh)...
Definition: coef.F90:547
parallel::lglel
integer function lglel(iel)
Definition: parallel_mod.F90:176
copy
subroutine copy(a, b, n)
Definition: math.F90:52
chkjac
subroutine chkjac(jac, n, iel, X, Y, Z, IERR)
Check the array JAC for a change in sign.
Definition: coef.F90:872
speclib::zwgj
subroutine zwgj(Z, W, NP, ALPHA, BETA)
Definition: speclib.F90:113
outpost
subroutine outpost(v1, v2, v3, vp, vt, name3)
Definition: prepost.F90:1152
speclib::dgll
subroutine, public dgll(D, DT, Z, NZ, NZD)
Definition: speclib.F90:867
parallel
cleaned
Definition: parallel_mod.F90:2
i
for i
Definition: xxt_test.m:74
facexv
subroutine facexv(A1, A2, A3, B1, B2, B3, IFACE1, IOP)
Definition: subs1.F90:903
geom
Geometry arrays.
Definition: geom_mod.F90:2
speclib::zwgll
subroutine, public zwgll(Z, W, NP)
Definition: speclib.F90:95
lagmass
subroutine lagmass()
Lag the mass matrix (matrices)
Definition: coef.F90:1144
setinvm
subroutine setinvm()
Invert the mass matrix. 1) Copy BM1 to BINVM1 2) Perform direct stiffness summation on BINVM1 3) Comp...
Definition: coef.F90:1168
volume
subroutine volume
Compute the volume based on mesh M1 and mesh M2.
Definition: coef.F90:914
geom2
subroutine geom2
Routine to generate all elemental geometric data for mesh 2 (Gauss-Legendre mesh). RXM2, RYM2, RZM2 - dr/dx, dr/dy, dr/dz SXM2, SYM2, SZM2 - ds/dx, ds/dy, ds/dz TXM2, TYM2, TZM2 - dt/dx, dt/dy, dt/dz JACM2 - Jacobian BM2 - Mass matrix.
Definition: coef.F90:725
speclib::zwglj
subroutine, public zwglj(Z, W, NP, ALPHA, BETA)
Definition: speclib.F90:204
unitvec
subroutine unitvec(X, Y, Z, N)
Definition: bdry.F90:1061
setarea
subroutine setarea(xrm1, yrm1, zrm1, xsm1, ysm1, zsm1, xtm1, ytm1, ztm1)
Compute surface data: areas, normals and tangents.
Definition: coef.F90:956
geom1
subroutine geom1()
Routine to generate all elemental geometric data for mesh 1. Velocity formulation : global-to-local m...
Definition: coef.F90:383
speclib::zwgl
subroutine, public zwgl(Z, W, NP)
Generate NP Gauss Legendre points (Z) and weights (W) associated with Jacobi polynomial P(N)(alpha=0...
Definition: speclib.F90:77
speclib::igllm
subroutine, public igllm(I12, IT12, Z1, Z2, NZ1, NZ2, ND1, ND2)
Definition: speclib.F90:1197
ixyz
Definition: ixyz_mod.F90:2
speclib::dglj
subroutine dglj(D, DT, Z, NZ, NZD, ALPHA, BETA)
Definition: speclib.F90:772
Generated on Tue Jul 14 2015 10:52:58 for Nek5000 by   doxygen 1.8.8