source: CIVL/examples/fortran/nek5000/core/subs2.f

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c

  • Property mode set to 100755
File size: 63.3 KB
Line 
1 SUBROUTINE SETAXW1 (IFAXWG)
2C
3 INCLUDE 'SIZE'
4 INCLUDE 'WZ'
5C
6 LOGICAL IFAXWG
7C
8 IF (IFAXWG) THEN
9 CALL COPY (W3M1,W2AM1,lx1*ly1)
10 ELSE
11 CALL COPY (W3M1,W2CM1,lx1*ly1)
12 ENDIF
13C
14 RETURN
15 END
16 SUBROUTINE SETAXW2 (IFAXWG)
17C-----------------------------------------------------------------------
18 INCLUDE 'SIZE'
19 INCLUDE 'WZ'
20C
21 LOGICAL IFAXWG
22C
23 IF (IFAXWG) THEN
24 CALL COPY (W3M2,W2AM2,lx2*ly2)
25 ELSE
26 CALL COPY (W3M2,W2CM2,lx2*ly2)
27 ENDIF
28C
29 RETURN
30 END
31C-----------------------------------------------------------------------
32 SUBROUTINE STNRINV
33C
34C Calculate 2nd and 3rd strain-rate invariants
35C
36 INCLUDE 'SIZE'
37 INCLUDE 'SOLN'
38 INCLUDE 'TSTEP'
39 common /screv/ ei2(lx1,ly1,lz1,lelt)
40 $ , ei3(lx1,ly1,lz1,lelt)
41 common /ctmp1/ exx(lx1,ly1,lz1,lelt)
42 $ , exy(lx1,ly1,lz1,lelt)
43 $ , eyy(lx1,ly1,lz1,lelt)
44 $ , ezz(lx1,ly1,lz1,lelt)
45 common /ctmp0/ exz(lx1,ly1,lz1,lelt)
46 $ , eyz(lx1,ly1,lz1,lelt)
47c
48 NTOT1 = lx1*ly1*lz1*NELV
49 CALL RZERO (EI2,NTOT1)
50 CALL RZERO (EI3,NTOT1)
51 IF (ISTEP.EQ.0) RETURN
52C
53 MATMOD = 0
54 CALL STNRATE (VX,VY,VZ,NELV,MATMOD)
55C
56 IF (ldim.EQ.2) THEN
57 CALL COL3 (EI2,EXX,EYY,NTOT1)
58 CALL SUBCOL3 (EI2,EXY,EXY,NTOT1)
59 CALL RZERO (EI3,NTOT1)
60 ELSE
61 CONST = 2.0
62 CALL COL4 (EI3,EXX,EYY,EZZ,NTOT1)
63 CALL COL4 (EI2,EXY,EXZ,EYZ,NTOT1)
64 CALL ADD2S2 (EI3,EI2,CONST,NTOT1)
65 CALL SUBCOL4 (EI3,EXX,EYZ,EYZ,NTOT1)
66 CALL SUBCOL4 (EI3,EYY,EXZ,EXZ,NTOT1)
67 CALL SUBCOL4 (EI3,EZZ,EXY,EXY,NTOT1)
68 CALL COL3 (EI2,EXX,EYY,NTOT1)
69 CALL ADDCOL3 (EI2,EXX,EZZ,NTOT1)
70 CALL ADDCOL3 (EI2,EYY,EZZ,NTOT1)
71 CALL SUBCOL3 (EI2,EXY,EXY,NTOT1)
72 CALL SUBCOL3 (EI2,EXZ,EXZ,NTOT1)
73 CALL SUBCOL3 (EI2,EYZ,EYZ,NTOT1)
74 ENDIF
75C
76 RETURN
77 END
78C-----------------------------------------------------------------------
79 SUBROUTINE OPDOT (DP,A1,A2,A3,B1,B2,B3,N)
80C
81 INCLUDE 'SIZE'
82C
83 DIMENSION DP(LX1,LY1,LZ1,1)
84 $ , A1(LX1,LY1,LZ1,1)
85 $ , A2(LX1,LY1,LZ1,1)
86 $ , A3(LX1,LY1,LZ1,1)
87 $ , B1(LX1,LY1,LZ1,1)
88 $ , B2(LX1,LY1,LZ1,1)
89 $ , B3(LX1,LY1,LZ1,1)
90C
91 IF (ldim.EQ.2) THEN
92 CALL VDOT2 (DP,A1,A2,B1,B2,N)
93 ELSE
94 CALL VDOT3 (DP,A1,A2,A3,B1,B2,B3,N)
95 ENDIF
96C
97 RETURN
98 END
99C-----------------------------------------------------------------------
100 SUBROUTINE OPADDS (A1,A2,A3,B1,B2,B3,CONST,N,ISC)
101C
102 INCLUDE 'SIZE'
103C
104 DIMENSION A1(LX1,LY1,LZ1,1)
105 $ , A2(LX1,LY1,LZ1,1)
106 $ , A3(LX1,LY1,LZ1,1)
107 $ , B1(LX1,LY1,LZ1,1)
108 $ , B2(LX1,LY1,LZ1,1)
109 $ , B3(LX1,LY1,LZ1,1)
110C
111 IF (ISC.EQ.1) THEN
112 CALL ADD2S1 (A1,B1,CONST,N)
113 CALL ADD2S1 (A2,B2,CONST,N)
114 IF (ldim.EQ.3) CALL ADD2S1 (A3,B3,CONST,N)
115 ELSEIF (ISC.EQ.2) THEN
116 CALL ADD2S2 (A1,B1,CONST,N)
117 CALL ADD2S2 (A2,B2,CONST,N)
118 IF (ldim.EQ.3) CALL ADD2S2 (A3,B3,CONST,N)
119 ENDIF
120C
121 RETURN
122 END
123C-----------------------------------------------------------------------
124 SUBROUTINE FACEXS (A,B,IFACE1,IOP)
125C
126C IOP = 0
127C Extract scalar A from B on face IFACE1.
128C
129C IOP = 1
130C Extract scalar B from A on face IFACE1.
131C
132C A has the (NX,NY,NFACE) data structure
133C B has the (NX,NY,NZ) data structure
134C IFACE1 is in the preprocessor notation
135C IFACE is the dssum notation.
136C
137 INCLUDE 'SIZE'
138 INCLUDE 'TOPOL'
139C
140 DIMENSION A(LX1,LY1),B(LX1,LY1,LZ1)
141C
142 CALL DSSET(lx1,ly1,lz1)
143 IFACE = EFACE1(IFACE1)
144 JS1 = SKPDAT(1,IFACE)
145 JF1 = SKPDAT(2,IFACE)
146 JSKIP1 = SKPDAT(3,IFACE)
147 JS2 = SKPDAT(4,IFACE)
148 JF2 = SKPDAT(5,IFACE)
149 JSKIP2 = SKPDAT(6,IFACE)
150C
151 I = 0
152 IF (IOP.EQ.0) THEN
153 DO 100 J2=JS2,JF2,JSKIP2
154 DO 100 J1=JS1,JF1,JSKIP1
155 I = I+1
156 A(I,1) = B(J1,J2,1)
157 100 CONTINUE
158 ELSE
159 DO 150 J2=JS2,JF2,JSKIP2
160 DO 150 J1=JS1,JF1,JSKIP1
161 I = I+1
162 B(J1,J2,1) = A(I,1)
163 150 CONTINUE
164 ENDIF
165C
166 RETURN
167 END
168C-----------------------------------------------------------------------
169 SUBROUTINE FACEXV (A1,A2,A3,B1,B2,B3,IFACE1,IOP)
170C
171C IOP = 0
172C Extract vector (A1,A2,A3) from (B1,B2,B3) on face IFACE1.
173C
174C IOP = 1
175C Extract vector (B1,B2,B3) from (A1,A2,A3) on face IFACE1.
176C
177C A1, A2, A3 have the (NX,NY,NFACE) data structure
178C B1, B2, B3 have the (NX,NY,NZ) data structure
179C IFACE1 is in the preprocessor notation
180C IFACE is the dssum notation.
181C
182 INCLUDE 'SIZE'
183 INCLUDE 'TOPOL'
184C
185 DIMENSION A1(LX1,LY1),A2(LX1,LY1),A3(LX1,LY1),
186 $ B1(LX1,LY1,LZ1),B2(LX1,LY1,LZ1),B3(LX1,LY1,LZ1)
187C
188 CALL DSSET(lx1,ly1,lz1)
189 IFACE = EFACE1(IFACE1)
190 JS1 = SKPDAT(1,IFACE)
191 JF1 = SKPDAT(2,IFACE)
192 JSKIP1 = SKPDAT(3,IFACE)
193 JS2 = SKPDAT(4,IFACE)
194 JF2 = SKPDAT(5,IFACE)
195 JSKIP2 = SKPDAT(6,IFACE)
196 I = 0
197C
198 IF (IOP.EQ.0) THEN
199 DO 100 J2=JS2,JF2,JSKIP2
200 DO 100 J1=JS1,JF1,JSKIP1
201 I = I+1
202 A1(I,1) = B1(J1,J2,1)
203 A2(I,1) = B2(J1,J2,1)
204 A3(I,1) = B3(J1,J2,1)
205 100 CONTINUE
206 ELSE
207 DO 150 J2=JS2,JF2,JSKIP2
208 DO 150 J1=JS1,JF1,JSKIP1
209 I = I+1
210 B1(J1,J2,1) = A1(I,1)
211 B2(J1,J2,1) = A2(I,1)
212 B3(J1,J2,1) = A3(I,1)
213 150 CONTINUE
214 ENDIF
215C
216 RETURN
217 END
218C-----------------------------------------------------------------------
219 SUBROUTINE FACSUB2 (A1,A2,A3,B1,B2,B3,IFACE1)
220C
221C Subtract B1,B2,B3 from A1,A2,A3 on surface IFACE1 of element IE.
222C
223C A1, A2, A3 have the (NX,NY,NZ) data structure
224C B1, B2, B3 have the (NX,NY,NFACE) data structure
225C IFACE1 is in the preprocessor notation
226C IFACE is the dssum notation.
227C
228 INCLUDE 'SIZE'
229 INCLUDE 'TOPOL'
230C
231 DIMENSION A1(LX1,LY1,LZ1),A2(LX1,LY1,LZ1),A3(LX1,LY1,LZ1),
232 $ B1(LX1,LY1),B2(LX1,LY1),B3(LX1,LY1)
233C
234 CALL DSSET(lx1,ly1,lz1)
235 IFACE = EFACE1(IFACE1)
236 JS1 = SKPDAT(1,IFACE)
237 JF1 = SKPDAT(2,IFACE)
238 JSKIP1 = SKPDAT(3,IFACE)
239 JS2 = SKPDAT(4,IFACE)
240 JF2 = SKPDAT(5,IFACE)
241 JSKIP2 = SKPDAT(6,IFACE)
242C
243 I = 0
244 DO 100 J2=JS2,JF2,JSKIP2
245 DO 100 J1=JS1,JF1,JSKIP1
246 I = I+1
247 A1(J1,J2,1) = A1(J1,J2,1) - B1(I,1)
248 A2(J1,J2,1) = A2(J1,J2,1) - B2(I,1)
249 A3(J1,J2,1) = A3(J1,J2,1) - B3(I,1)
250 100 CONTINUE
251C
252 RETURN
253 END
254C-----------------------------------------------------------------------
255 SUBROUTINE GAMMASF (H1,H2)
256C-----------------------------------------------------------------------
257C
258C Compute lagrest eigenvalue of coupled Helmholtz operator
259C
260C-----------------------------------------------------------------------
261 INCLUDE 'SIZE'
262 INCLUDE 'EIGEN'
263 INCLUDE 'INPUT'
264 INCLUDE 'MASS'
265 INCLUDE 'MVGEOM'
266 INCLUDE 'SOLN'
267 INCLUDE 'TSTEP'
268 INCLUDE 'WZ'
269 common /scrmg/ ae1(lx1,ly1,lz1,lelv)
270 $ , ae2(lx1,ly1,lz1,lelv)
271 $ , ae3(lx1,ly1,lz1,lelv)
272 common /scruz/ e1(lx1,ly1,lz1,lelv)
273 $ , e2(lx1,ly1,lz1,lelv)
274 $ , e3(lx1,ly1,lz1,lelv)
275C
276 DIMENSION H1(LX1,LY1,LZ1,1),H2(LX1,LY1,LZ1,1)
277C
278 NTOT1 = lx1*ly1*lz1*NELV
279 IMESH = 1
280 MATMOD = 0
281C
282 IF (ISTEP.EQ.0) THEN
283 EIGGA = 0.0
284 CALL STX1SF
285 ELSE
286 CALL COPY (E1,EV1,NTOT1)
287 CALL COPY (E2,EV2,NTOT1)
288 IF (ldim.EQ.3) CALL COPY (E3,EV3,NTOT1)
289 ENDIF
290C
291 EVNEW = EIGGA
292C
293 DO 1000 ITER=1,NMXE
294C
295 CALL AXHMSF (AE1,AE2,AE3,E1,E2,E3,H1,H2,MATMOD)
296 CALL RMASK (AE1,AE2,AE3,NELV)
297 CALL OPDSSUM (AE1,AE2,AE3)
298C
299 EVOLD = EVNEW
300 EVNEW = GLSC3(E1,AE1,VMULT,NTOT1) + GLSC3(E2,AE2,VMULT,NTOT1)
301 IF (ldim.EQ.3) EVNEW = EVNEW + GLSC3(E3,AE3,VMULT,NTOT1)
302 CRIT = ABS( (EVNEW - EVOLD)/EVNEW )
303 IF ( CRIT .LT. TOLEV ) GOTO 2000
304C
305 CALL COL3 (E1,BINVM1,AE1,NTOT1)
306 CALL COL3 (E2,BINVM1,AE2,NTOT1)
307 IF (ldim.EQ.3) CALL COL3 (E3,BINVM1,AE3,NTOT1)
308 XX = GLSC3(E1,AE1,VMULT,NTOT1) + GLSC3(E2,AE2,VMULT,NTOT1)
309 IF (ldim.EQ.3) XX = XX + GLSC3(E3,AE3,VMULT,NTOT1)
310 IF (XX .LT. 0.0) GO TO 9000
311C
312 XNORM=1./SQRT( XX )
313 CALL CMULT (E1,XNORM,NTOT1)
314 CALL CMULT (E2,XNORM,NTOT1)
315 IF (ldim.EQ.3) CALL CMULT (E3,XNORM,NTOT1)
316C
317 1000 CONTINUE
318C
319C Save eigenvalue for all cases.
320C Save eigenvectors for deforming geometries.
321C
322 2000 EIGGA = EVNEW
323 IF (IFMVBD) THEN
324 CALL COPY (EV1,E1,NTOT1)
325 CALL COPY (EV2,E2,NTOT1)
326 IF (ldim.EQ.3) CALL COPY (EV3,E3,NTOT1)
327 ENDIF
328C
329 RETURN
330C
331 9000 CONTINUE
332 IF (NID.EQ.0)
333 $WRITE ( 6,*) ' Non +ve def. A-operator detected during eigenvalue
334 $computation : tran(x)Ax =',XX
335 CALL EMERXIT
336 call exitt
337C
338 END
339C-----------------------------------------------------------------------
340 SUBROUTINE CMULT2 (A,B,CONST,N)
341 DIMENSION A(1),B(1)
342 DO 100 I=1,N
343 A(I)=B(I)*CONST
344 100 CONTINUE
345 RETURN
346 END
347C-----------------------------------------------------------------------
348 SUBROUTINE ADD3S (A,B,C,CONST,N)
349 DIMENSION A(1),B(1),C(1)
350 DO 100 I=1,N
351 A(I)=B(I)+CONST*C(I)
352 100 CONTINUE
353 RETURN
354 END
355C-----------------------------------------------------------------------
356 SUBROUTINE EMERXIT
357C
358 INCLUDE 'SIZE'
359 INCLUDE 'INPUT'
360 INCLUDE 'TSTEP'
361 INCLUDE 'PARALLEL'
362C
363C Try to hang in there on the first few time steps (pff 8/92)
364 IF (IFTRAN.AND.ISTEP.LT.9) RETURN
365C
366 LASTEP = 1
367 CALL PREPOST(.true.,' ')
368C
369 IF (NP.EQ.1) THEN
370 WRITE (6,*) ' '
371 WRITE (6,*)
372 $ ' Emergency exit:',ISTEP,' time =',TIME
373 WRITE (6,*)
374 $ ' Latest solution and data are dumped for post-processing.'
375 WRITE (6,*) ' *** STOP ***'
376 ELSE
377 WRITE (6,*) ' '
378 WRITE (6,*) NID,
379 $ ' Emergency exit:',ISTEP,' time =',TIME
380 WRITE (6,*)
381 $ ' Latest solution and data are dumped for post-processing.'
382 WRITE (6,*) ' *** STOP ***'
383 ENDIF
384C
385 call runstat
386c
387 call exitt
388 RETURN
389 END
390C-----------------------------------------------------------------------
391 SUBROUTINE FACCVS (A1,A2,A3,B,IFACE1)
392C
393C Collocate scalar B with vector A, components A1,A2,A3,
394C on the surface IFACE1 of an element.
395C
396C A1,A2,A3 have the (NX,NY,NZ) data structure
397C B has the (NX,NY,IFACE) data structure
398C IFACE1 is in the preprocessor notation
399C IFACE is the dssum notation.
400C
401 INCLUDE 'SIZE'
402 INCLUDE 'TOPOL'
403 DIMENSION A1(LX1,LY1,LZ1),A2(LX1,LY1,LZ1),A3(LX1,LY1,LZ1),
404 $ B(LX1,LY1)
405C
406C Set up counters
407C
408 CALL DSSET(lx1,ly1,lz1)
409 IFACE = EFACE1(IFACE1)
410 JS1 = SKPDAT(1,IFACE)
411 JF1 = SKPDAT(2,IFACE)
412 JSKIP1 = SKPDAT(3,IFACE)
413 JS2 = SKPDAT(4,IFACE)
414 JF2 = SKPDAT(5,IFACE)
415 JSKIP2 = SKPDAT(6,IFACE)
416 I = 0
417C
418 IF (ldim.EQ.2) THEN
419 DO 100 J2=JS2,JF2,JSKIP2
420 DO 100 J1=JS1,JF1,JSKIP1
421 I = I+1
422 A1(J1,J2,1) = A1(J1,J2,1)*B(I,1)
423 A2(J1,J2,1) = A2(J1,J2,1)*B(I,1)
424 100 CONTINUE
425 ELSE
426 DO 200 J2=JS2,JF2,JSKIP2
427 DO 200 J1=JS1,JF1,JSKIP1
428 I = I+1
429 A1(J1,J2,1) = A1(J1,J2,1)*B(I,1)
430 A2(J1,J2,1) = A2(J1,J2,1)*B(I,1)
431 A3(J1,J2,1) = A3(J1,J2,1)*B(I,1)
432 200 CONTINUE
433 ENDIF
434C
435 RETURN
436 END
437C-----------------------------------------------------------------------
438 SUBROUTINE STX1SF
439C------------------------------------------------------------------
440C
441C Compute startvector for finding an eigenvalue on mesh 1.
442C Normalization: XT*B*X = 1
443C
444C------------------------------------------------------------------
445 INCLUDE 'SIZE'
446 INCLUDE 'MASS'
447 INCLUDE 'SOLN'
448 common /scrmg/ ae1(lx1,ly1,lz1,lelv)
449 $ , ae2(lx1,ly1,lz1,lelv)
450 $ , ae3(lx1,ly1,lz1,lelv)
451 common /scruz/ e1(lx1,ly1,lz1,lelv)
452 $ , e2(lx1,ly1,lz1,lelv)
453 $ , e3(lx1,ly1,lz1,lelv)
454C
455 NTOT1 = lx1*ly1*lz1*NELV
456 CALL RZERO3 (E1 ,E2 ,E3 ,NTOT1)
457 CALL RZERO3 (AE1,AE2,AE3,NTOT1)
458C
459 CALL COPY (E1,BM1,NTOT1)
460 CALL COPY (E2,BM1,NTOT1)
461 IF (ldim.EQ.3) CALL COPY (E3,BM1,NTOT1)
462C
463 CALL RMASK (E1,E2,E3,NELV)
464 CALL COL3 (AE1,BM1,E1,NTOT1)
465 CALL COL3 (AE2,BM1,E2,NTOT1)
466 IF (ldim.EQ.3) CALL COL3 (AE3,BM1,E3,NTOT1)
467C
468 CALL OPDSSUM (AE1,AE2,AE3)
469C
470 XX = GLSC3 (E1,AE1,VMULT,NTOT1) + GLSC3 (E2,AE2,VMULT,NTOT1)
471 IF (ldim.EQ.3) XX = XX + GLSC3 (E3,AE3,VMULT,NTOT1)
472 XNORM = 1./SQRT(XX)
473 CALL CMULT (E1,XNORM,NTOT1)
474 CALL CMULT (E2,XNORM,NTOT1)
475 IF (ldim.EQ.3) CALL CMULT (E3,XNORM,NTOT1)
476
477c call exitti ('quit in stx1sf$,',nel)
478
479C
480 RETURN
481 END
482C-----------------------------------------------------------------------
483 SUBROUTINE SOLVEL
484C
485 INCLUDE 'SIZE'
486 INCLUDE 'GEOM'
487 INCLUDE 'SOLN'
488 INCLUDE 'TSTEP'
489C
490 DO 100 IEL=1,NELV
491 DO 100 K=1,lz1
492 DO 100 J=1,ly1
493 DO 100 I=1,lx1
494 CALL VSOLN (VX (I,J,K,IEL),VY (I,J,K,IEL),VZ (I,J,K,IEL),
495 $ XM1(I,J,K,IEL),YM1(I,J,K,IEL),ZM1(I,J,K,IEL),PI)
496 100 CONTINUE
497C
498 RETURN
499 END
500C-----------------------------------------------------------------------
501 SUBROUTINE VSOLN (UX,UY,UZ,X,Y,Z,PI)
502C
503C URR=(1.-0.75/SQRT(X**2+Y**2)+0.0625/(SQRT(X**2+Y**2)**3))*(X/
504C $ SQRT(X**2+Y**2))
505C UTETA=-(1.-0.375/SQRT(X**2+Y**2)-0.03125/(SQRT(X**2+Y**2)**3))*
506C $ (Y/SQRT(X**2+Y**2))
507C UX=URR*(X/SQRT(X**2+Y**2))-UTETA*(Y/SQRT(X**2+Y**2))
508C UY=URR*(Y/SQRT(X**2+Y**2))+UTETA*(X/SQRT(X**2+Y**2))
509C
510C UX = 2.*COS( PI*X )
511C UY = PI*Y*SIN( PI*X )
512C
513 UX = 0.0
514 UY = 0.0
515C
516 RETURN
517 END
518C-----------------------------------------------------------------------
519 SUBROUTINE SOLPRES
520C
521 INCLUDE 'SIZE'
522 INCLUDE 'GEOM'
523 INCLUDE 'SOLN'
524 INCLUDE 'TSTEP'
525C
526 DO 100 IEL=1,NELV
527 DO 100 K=1,lz2
528 DO 100 J=1,ly2
529 DO 100 I=1,lx2
530 CALL PRSOLN (PR (I,J,K,IEL),XM2(I,J,K,IEL),YM2(I,J,K,IEL),
531 * ZM2(I,J,K,IEL),PI)
532 100 CONTINUE
533C
534 RETURN
535 END
536C-----------------------------------------------------------------------
537 SUBROUTINE PRSOLN (P,X,Y,Z,PI)
538C
539C R = SQRT( X**2 + Y**2 )
540C CS = X/R
541C P = -0.75 * CS / R**2
542C
543C P = -SIN( PI*X )*COS( PI*Y )
544C
545 P = 0.0
546C
547 RETURN
548 END
549C-----------------------------------------------------------------------
550 SUBROUTINE PRINTEL (TA,A,IEL)
551C
552 INCLUDE 'SIZE'
553 DIMENSION TA(LX1,LY1,LZ1,LELT)
554 CHARACTER A*10
555C
556 lz1I = 1
557 lz1J = lz1
558 lz1INC = 1
559C
560 WRITE (21,*) 'ELEMENT NUMBER ',IEL
561 DO 101 IPL=lz1I,lz1J,lz1INC
562 CALL OUTM1 (TA,A,lz1,IEL,IPL)
563 101 CONTINUE
564C
565 RETURN
566 END
567C-----------------------------------------------------------------------
568 SUBROUTINE PRINTV (TA,A,NEL)
569C-----------------------------------------------------------------------
570C
571C Store the results
572C
573C-----------------------------------------------------------------------
574 INCLUDE 'SIZE'
575 DIMENSION TA(LX1,LY1,LZ1,LELT)
576 CHARACTER A*10
577C
578 lz1I = 1
579 lz1J = lz1
580 lz1INC = 1
581C
582 DO 9001 IEL = 1,NEL
583 WRITE (21,*) 'ELEMENT NUMBER ',IEL
584 DO 101 IPL=lz1I,lz1J,lz1INC
585 CALL OUTM1 (TA,A,lz1,IEL,IPL)
586 101 CONTINUE
587 9001 CONTINUE
588C
589 RETURN
590 END
591C-----------------------------------------------------------------------
592 SUBROUTINE OUTF1 (X,TXT,IEL,IFC)
593 INCLUDE 'SIZE'
594 DIMENSION X(LX1,LZ1,6,LELT)
595 CHARACTER*10 TXT
596C
597 NFACE = 2*ldim
598 NZI = lz1
599 NZJ = 1
600 NZINC = -1
601 NXI = 1
602 NXJ = lx1
603 NXINC = 1
604C
605 WRITE(21,106) TXT,IFC,NFACE
606 DO 100 J=NZI,NZJ,NZINC
607 WRITE(21,105) (X(I,J,IFC,IEL),I=NXI,NXJ,NXINC)
608 100 CONTINUE
609C
610 105 FORMAT(5E15.6)
611 106 FORMAT(///,5X,' ^ ',/,
612 $ 5X,' S | ',/,
613 $ 5X,' | ',A10,/,
614 $ 5X,' +----> ','Plane = ',I2,'/',I2,/,
615 $ 5X,' R ',/)
616C
617 RETURN
618 END
619C-----------------------------------------------------------------------
620 SUBROUTINE OUTM1 (X,TXT,NP,IEL,IP)
621 INCLUDE 'SIZE'
622 DIMENSION X(LX1,LY1,LZ1,LELT)
623 CHARACTER*10 TXT
624C
625 NYI = ly1
626 NYJ = 1
627 NYINC = -1
628 NXI = 1
629 NXJ = lx1
630 NXINC = 1
631C
632 WRITE(6,106) TXT,IP,NP
633 DO 100 J=NYI,NYJ,NYINC
634 WRITE(6,105) (X(I,J,IP,IEL),I=NXI,NXJ,NXINC)
635 100 CONTINUE
636C
637c 105 FORMAT(1p8e10.3)
638 105 FORMAT(8f10.3)
639 106 FORMAT(///,5X,' ^ ',/,
640 $ 5X,' Y | ',/,
641 $ 5X,' | ',A10,/,
642 $ 5X,' +----> ','Plane = ',I2,'/',I2,/,
643 $ 5X,' X ',/)
644C
645 RETURN
646 END
647C-----------------------------------------------------------------------
648 SUBROUTINE OUTM2 (X,TXT,NP,IEL,IP)
649 INCLUDE 'SIZE'
650 DIMENSION X(LX2,LY2,LZ2,LELV)
651 CHARACTER*10 TXT
652C
653 NYI = ly2
654 NYJ = 1
655 NYINC = -1
656 NXI = 1
657 NXJ = lx2
658 NXINC = 1
659C
660 WRITE(21,106) TXT,IP,NP
661 DO 100 J=NYI,NYJ,NYINC
662 WRITE(21,105) (X(I,J,IP,IEL),I=NXI,NXJ,NXINC)
663 100 CONTINUE
664C
665 105 FORMAT(5E15.6)
666 106 FORMAT(///,5X,' ^ ',/,
667 $ 5X,' Y | ',/,
668 $ 5X,' | ',A10,/,
669 $ 5X,' +----> ','Plane = ',I2,'/',I2,/,
670 $ 5X,' X ',/)
671C
672 RETURN
673 END
674C-----------------------------------------------------------------------
675 SUBROUTINE STSMASK (C1MASK,C2MASK,C3MASK)
676C
677 INCLUDE 'SIZE'
678 INCLUDE 'GEOM'
679 INCLUDE 'TSTEP'
680 include 'INPUT'
681 common /screv/ hfmask(lx1,lz1,6,lelt)
682 $ , hvmask(lx1,ly1,lz1,lelt)
683C
684 DIMENSION C1MASK(LX1,LY1,LZ1,1)
685 $ , C2MASK(LX1,LY1,LZ1,1)
686 $ , C3MASK(LX1,LY1,LZ1,1)
687 INTEGER IMDATA
688 SAVE IMDATA
689 DATA IMDATA /0/
690C
691 IFLD = IFIELD
692 NEL = NELFLD(IFIELD)
693C
694 IF (IMDATA.EQ.0) THEN
695 CALL SETCDAT
696 IMDATA=1
697 ENDIF
698C
699 IF (IFLD.EQ.1) CALL SKIPCNR (NEL)
700 CALL SETHMSK (HVMASK,HFMASK,IFLD,NEL)
701 CALL SETMLOG (HVMASK,HFMASK,IFLD,NEL)
702 CALL SETMASK (C1MASK,C2MASK,C3MASK,HVMASK,NEL)
703 IF (IFLMSF(IFLD)) CALL SETCSYS (HVMASK,HFMASK,NEL)
704 IF (IFLD.EQ.0) CALL FIXWMSK (C2MASK,C3MASK,HVMASK,HFMASK,NEL)
705
706 if (ifaxis.and.ifld.eq.1) call fixmska (c1mask,c2mask,c3mask)
707C
708 RETURN
709 END
710C-----------------------------------------------------------------------
711 SUBROUTINE UPDMSYS (IFLD)
712C
713 INCLUDE 'SIZE'
714 INCLUDE 'GEOM'
715 INCLUDE 'TSTEP'
716 common /screv/ hfmask(lx1,lz1,6,lelt)
717 $ , hvmask(lx1,ly1,lz1,lelt)
718C
719 IF (.NOT.IFLMSF(IFLD)) RETURN
720C
721 NEL = NELFLD(IFLD)
722 CALL SETHMSK (HVMASK,HFMASK,IFLD,NEL)
723 CALL SETCSYS (HVMASK,HFMASK,NEL)
724C
725 RETURN
726 END
727C-----------------------------------------------------------------------
728 SUBROUTINE SETHMSK (HVMASK,HFMASK,IFLD,NEL)
729C
730 INCLUDE 'SIZE'
731 INCLUDE 'INPUT'
732 INCLUDE 'TSTEP'
733C
734 DIMENSION HVMASK(LX1,LY1,LZ1,1)
735 $ , HFMASK(LX1,LZ1,6,1)
736 CHARACTER CB*3
737C
738 NTOT1 = lx1*ly1*lz1*NEL
739 NXZ1 = lx1*lz1
740 NTOTF = NXZ1*6*NEL
741 NFACE = 2*ldim
742 CONST = 5.0
743 CALL CFILL (HVMASK,CONST,NTOT1)
744 CALL CFILL (HFMASK,CONST,NTOTF)
745C
746 IF (IFLD.EQ.1) THEN
747C
748 DO 110 IEL=1,NEL
749 DO 110 IFC=1,NFACE
750 CB=CBC(IFC,IEL,IFLD)
751 IF (CB.EQ.'ON ' .OR. CB.EQ.'on ' .or.
752 $ CB.EQ.'MM ' .OR. CB.EQ.'mm ' ) THEN
753 CALL FACEV (HVMASK,IEL,IFC,3.0,lx1,ly1,lz1)
754 CALL CFILL (HFMASK(1,1,IFC,IEL),3.0,NXZ1)
755 ENDIF
756 110 CONTINUE
757C
758 DO 120 IEL=1,NEL
759 DO 120 IFC=1,NFACE
760 CB=CBC(IFC,IEL,IFLD)
761 IF (CB.EQ.'SYM' .OR. CB.EQ.'A ' .OR. CB.EQ.'WS ' .OR.
762 $ CB.EQ.'ws ' .OR. CB.EQ.'WSL' .OR. CB.EQ.'wsl' .OR.
763 $ CB.EQ.'SH ' .OR. CB.EQ.'sh ' .OR. CB.EQ.'SHL' .OR.
764 $ CB.EQ.'shl') THEN
765 CALL FACEV (HVMASK,IEL,IFC,2.0,lx1,ly1,lz1)
766 CALL CFILL (HFMASK(1,1,IFC,IEL),2.0,NXZ1)
767 ENDIF
768 120 CONTINUE
769C
770 DO 130 IEL=1,NEL
771 DO 130 IFC=1,NFACE
772 CB=CBC(IFC,IEL,IFLD)
773 IF (CB.EQ.'MF ' .OR. CB.EQ.'V ' .OR. CB.EQ.'v ' .OR.
774 $ CB.EQ.'VL ' .OR. CB.EQ.'vl ' .OR. CB(1:2).EQ.'mv') THEN
775 CALL FACEV (HVMASK,IEL,IFC,1.0,lx1,ly1,lz1)
776 CALL CFILL (HFMASK(1,1,IFC,IEL),1.0,NXZ1)
777 ENDIF
778 130 CONTINUE
779C
780 DO 140 IEL=1,NEL
781 DO 140 IFC=1,NFACE
782 CB=CBC(IFC,IEL,IFLD)
783 IF (CB.EQ.'W ') THEN
784 CALL FACEV (HVMASK,IEL,IFC,0.0,lx1,ly1,lz1)
785 CALL CFILL (HFMASK(1,1,IFC,IEL),0.0,NXZ1)
786 ENDIF
787 140 CONTINUE
788C
789 ELSE
790C
791 DO 210 IEL=1,NEL
792 DO 210 IFC=1,NFACE
793 CB=CBC(IFC,IEL,IFLD)
794 IF (CB.EQ.'SYM') THEN
795 CALL FACEV (HVMASK,IEL,IFC,2.0,lx1,ly1,lz1)
796 CALL CFILL (HFMASK(1,1,IFC,IEL),2.0,NXZ1)
797 ENDIF
798 210 CONTINUE
799C
800c write(6,*) 'MASK this is ifield:',ifield
801 DO 220 IEL=1,NEL
802 DO 220 IFC=1,NFACE
803 CB=CBC(IFC,IEL,IFLD)
804 IF (CB(1:1).EQ.'M' .OR. CB(1:1).EQ.'m') THEN
805c CALL FACEV (HVMASK,IEL,IFC,1.0,lx1,ly1,lz1)
806c CALL CFILL (HFMASK(1,1,IFC,IEL),1.0,NXZ1)
807 CALL FACEV (HVMASK,IEL,IFC,2.0,lx1,ly1,lz1)
808 CALL CFILL (HFMASK(1,1,IFC,IEL),2.0,NXZ1)
809 ENDIF
810 220 CONTINUE
811C
812 DO 230 IEL=1,NEL
813 DO 230 IFC=1,NFACE
814 CB=CBC(IFC,IEL,IFLD)
815 IF (CB.EQ.'FIX') THEN
816 CALL FACEV (HVMASK,IEL,IFC,0.0,lx1,ly1,lz1)
817 CALL CFILL (HFMASK(1,1,IFC,IEL),0.0,NXZ1)
818 ENDIF
819 230 CONTINUE
820C
821 ENDIF
822C
823 CALL DSOP (HVMASK,'MNA',lx1,ly1,lz1)
824C
825 RETURN
826 END
827C-----------------------------------------------------------------------
828 SUBROUTINE SKIPCNR (NEL)
829C
830 INCLUDE 'SIZE'
831 INCLUDE 'GEOM'
832 INCLUDE 'INPUT'
833 common /indxfc/ mcrfc(4,6)
834 $ , MFCCR(3,8)
835 $ , MEGCR(3,8)
836 $ , MFCEG(2,12)
837 $ , MCREG(2,12)
838 $ , MCRRST(3,8)
839 $ , MIDRST(3,12)
840 $ , MCRIND(8)
841 $ , MEDIND(2,4)
842 $ , NTEFC(2,12)
843 $ , NTCRF(2,3)
844C
845 NFACE = 2*ldim
846 NCRFC = NFACE - 2
847 NMXCR = 8*NEL
848 CALL LFALSE (IFNSKP,NMXCR)
849C
850 DO 100 IEL=1,NEL
851 DO 100 IFC=1,NFACE
852 IF (CDOF(IFC,IEL).EQ.'1') THEN
853 ICR=MCRFC(1,IFC)
854 IFNSKP(ICR,IEL)=.TRUE.
855 ELSEIF (CDOF(IFC,IEL).EQ.'2') THEN
856 ICR=MCRFC(2,IFC)
857 IFNSKP(ICR,IEL)=.TRUE.
858 ELSEIF (CDOF(IFC,IEL).EQ.'3') THEN
859 ICR=MCRFC(3,IFC)
860 IFNSKP(ICR,IEL)=.TRUE.
861 ELSEIF (CDOF(IFC,IEL).EQ.'4') THEN
862 ICR=MCRFC(4,IFC)
863 IFNSKP(ICR,IEL)=.TRUE.
864 ENDIF
865 IF (CDOF(IFC,IEL).EQ.'*') THEN
866 DO 160 ICRFC=1,NCRFC
867 ICR=MCRFC(ICRFC,IFC)
868 IFNSKP(ICR,IEL)=.TRUE.
869 160 CONTINUE
870 ENDIF
871 100 CONTINUE
872C
873 RETURN
874 END
875C-----------------------------------------------------------------------
876 SUBROUTINE SETMASK (C1MASK,C2MASK,C3MASK,HVMASK,NEL)
877C
878 INCLUDE 'SIZE'
879 INCLUDE 'TOTAL' !XXXXX
880C
881 DIMENSION HVMASK (LX1,LY1,LZ1,1)
882 $ , C1MASK(LX1,LY1,LZ1,1)
883 $ , C2MASK(LX1,LY1,LZ1,1)
884 $ , C3MASK(LX1,LY1,LZ1,1)
885C
886 NTOT1 = lx1*ly1*lz1*NEL
887 CALL RZERO3 (C1MASK,C2MASK,C3MASK,NTOT1)
888C
889 DO 100 IEL=1,NEL
890 DO 100 IZ=1,lz1
891 DO 100 IY=1,ly1
892 DO 100 IX=1,lx1
893 HMV=ABS( HVMASK(IX,IY,IZ,IEL) )
894 IF (HMV .GT. 2.9) THEN
895 C1MASK(IX,IY,IZ,IEL) = 1.0
896 ENDIF
897 IF ((HMV.GT.1.9 .AND. HMV.LT.2.1) .OR. HMV.GT.4.9) THEN
898 C2MASK(IX,IY,IZ,IEL) = 1.0
899 ENDIF
900 100 CONTINUE
901C
902 IF (ldim.EQ.3) CALL COPY (C3MASK,C2MASK,NTOT1)
903C
904 RETURN
905 END
906C-----------------------------------------------------------------------
907 SUBROUTINE SETMLOG (HVMASK,HFMASK,IFLD,NEL)
908C
909 INCLUDE 'SIZE'
910 INCLUDE 'GEOM'
911 common /indxfc/ mcrfc(4,6)
912 $ , MFCCR(3,8)
913 $ , MEGCR(3,8)
914 $ , MFCEG(2,12)
915 $ , MCREG(2,12)
916 $ , MCRRST(3,8)
917 $ , MIDRST(3,12)
918 $ , MCRIND(8)
919 $ , MEDIND(2,4)
920 $ , NTEFC(2,12)
921 $ , NTCRF(2,3)
922C
923 DIMENSION HVMASK(LX1,LY1,LZ1,1)
924 $ , HFMASK(LX1,LZ1,6,1)
925C
926 NFACE = 2*ldim
927 NEDGE = 12
928 NCRNR = 2**ldim
929 NTOTF = NFACE*NEL
930 NTOTS = NEDGE*NEL
931 NTOTC = NCRNR*NEL
932 EPSA = 1.E-6
933C
934 IFLMSF(IFLD) = .FALSE.
935 IFLMSE(IFLD) = .FALSE.
936 IFLMSC(IFLD) = .FALSE.
937 CALL LFALSE (IFMSFC(1,1,IFLD),NTOTF)
938 CALL LFALSE (IFMSEG(1,1,IFLD),NTOTS)
939 CALL LFALSE (IFMSCR(1,1,IFLD),NTOTC)
940C
941 DO 100 IEL=1,NEL
942 DO 100 IFC=1,NFACE
943 HMF = ABS( HFMASK(1,1,IFC,IEL) )
944 IF (HMF .GT. 1.9 .AND. HMF .LT. 3.1 ) THEN
945 IFLMSF(IFLD) = .TRUE.
946 IFMSFC(IFC,IEL,IFLD) = .TRUE.
947 ENDIF
948 100 CONTINUE
949 CALL GLLOG(IFLMSF(IFLD),.TRUE.)
950C
951 IF (ldim.EQ.3) THEN
952 DO 200 IEL=1,NEL
953 DO 200 ISD=1,NEDGE
954 IX = MIDRST(1,ISD)
955 IY = MIDRST(2,ISD)
956 IZ = MIDRST(3,ISD)
957 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
958 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 200
959 IDIFF = 0
960 DO 220 II=1,2
961 IFC = MFCEG(II,ISD)
962 HMF = ABS( HFMASK(1,1,IFC,IEL) )
963 IF (ABS(HMV - HMF) .GT. EPSA) IDIFF=IDIFF + 1
964 220 CONTINUE
965 IF (IDIFF.EQ.2) THEN
966 IFLMSE(IFLD) = .TRUE.
967 IFMSEG(ISD,IEL,IFLD) = .TRUE.
968 ENDIF
969 200 CONTINUE
970 CALL GLLOG(IFLMSE(IFLD),.TRUE.)
971 ENDIF
972C
973 DO 300 IEL=1,NEL
974 DO 300 ICR=1,NCRNR
975 IX = MCRRST(1,ICR)
976 IY = MCRRST(2,ICR)
977 IZ = MCRRST(3,ICR)
978 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
979 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 300
980 IDIFF = 0
981 DO 330 II=1,ldim
982 IFC = MFCCR(II,ICR)
983 HMF = ABS( HFMASK(1,1,IFC,IEL) )
984 IF (ABS(HMV - HMF) .GT. EPSA) IDIFF=IDIFF + 1
985 330 CONTINUE
986 IF (ldim.EQ.3) THEN
987 DO 360 II=1,ldim
988 ISD = MEGCR(II,ICR)
989 IXS = MIDRST(1,ISD)
990 IYS = MIDRST(2,ISD)
991 IZS = MIDRST(3,ISD)
992 HMS = ABS( HVMASK(IXS,IYS,IZS,IEL) )
993 IF (ABS(HMV - HMS) .GT. EPSA) IDIFF=IDIFF + 1
994 360 CONTINUE
995 ENDIF
996 IF ( (ldim.EQ.2 .AND. IDIFF.EQ.2) .OR.
997 $ (ldim.EQ.3 .AND. IDIFF.EQ.6) ) THEN
998 IFLMSC(IFLD) = .TRUE.
999 IFMSCR(ICR,IEL,IFLD) = .TRUE.
1000 ENDIF
1001 300 CONTINUE
1002 CALL GLLOG(IFLMSC(IFLD),.TRUE.)
1003C
1004 RETURN
1005 END
1006C-----------------------------------------------------------------------
1007 SUBROUTINE SETCSYS (HVMASK,HFMASK,NEL)
1008C
1009 INCLUDE 'SIZE'
1010 INCLUDE 'GEOM'
1011 INCLUDE 'TSTEP'
1012 DIMENSION HVMASK(LX1,LY1,LZ1,1)
1013 $ , HFMASK(LX1,LZ1,6,1)
1014C
1015 NFACE = 2*ldim
1016 NTOT1 = lx1*ly1*lz1*NEL
1017C
1018 CALL RZERO3 (VNX,VNY,VNZ,NTOT1)
1019 CALL RZERO3 (V1X,V1Y,V1Z,NTOT1)
1020 CALL RZERO3 (V2X,V2Y,V2Z,NTOT1)
1021C
1022 DO 10 IEL=1,NEL
1023 DO 10 IFC=1,NFACE
1024 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1025 IF (HMF .GT. 1.9 .AND. HMF .LT. 3.1)
1026 $ CALL FACEXV(UNX(1,1,IFC,IEL),UNY(1,1,IFC,IEL),UNZ(1,1,IFC,IEL),
1027 $ VNX(1,1,1,IEL),VNY(1,1,1,IEL),VNZ(1,1,1,IEL),IFC,1)
1028 10 CONTINUE
1029C
1030 IF (ldim.EQ.2) THEN
1031 CALL COMAVN2 (HVMASK,HFMASK,NEL)
1032 ELSE
1033 CALL COMAVN3 (HVMASK,HFMASK,NEL)
1034 ENDIF
1035C
1036 RETURN
1037 END
1038C-----------------------------------------------------------------------
1039 SUBROUTINE COMAVN2 (HVMASK,HFMASK,NEL)
1040C
1041 INCLUDE 'SIZE'
1042 INCLUDE 'GEOM'
1043 common /indxfc/ mcrfc(4,6)
1044 $ , MFCCR(3,8)
1045 $ , MEGCR(3,8)
1046 $ , MFCEG(2,12)
1047 $ , MCREG(2,12)
1048 $ , MCRRST(3,8)
1049 $ , MIDRST(3,12)
1050 $ , MCRIND(8)
1051 $ , MEDIND(2,4)
1052 $ , NTEFC(2,12)
1053 $ , NTCRF(2,3)
1054C
1055 DIMENSION HVMASK(LX1,LY1,LZ1,1)
1056 $ , HFMASK(LX1,LZ1,6,1)
1057C
1058 NTOT1 = lx1*ly1*lz1*NEL
1059 NFACE = 2*ldim
1060 NCRNR = 2**ldim
1061 EPSA = 1.0E-06
1062 CALL RZERO (VNZ,NTOT1)
1063C
1064 IZ = 1
1065 DO 100 IEL=1,NEL
1066 DO 100 ICR=1,NCRNR
1067 IX = MCRRST(1,ICR)
1068 IY = MCRRST(2,ICR)
1069 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1070 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 100
1071 VNX(IX,IY,IZ,IEL) = 0.0
1072 VNY(IX,IY,IZ,IEL) = 0.0
1073 100 CONTINUE
1074C
1075 DO 200 IEL=1,NEL
1076 DO 200 ICR=1,NCRNR
1077 IF (IFNSKP(ICR,IEL)) GOTO 200
1078 IX = MCRRST(1,ICR)
1079 IY = MCRRST(2,ICR)
1080 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1081 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 200
1082 DO 220 II=1,2
1083 IFC = MFCCR(II,ICR)
1084 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1085 IF (ABS(HMV - HMF) .LT. EPSA) THEN
1086 IR = MCRRST(II,ICR)
1087 VNX(IX,IY,IZ,IEL)=VNX(IX,IY,IZ,IEL) + UNX(IR,IZ,IFC,IEL)
1088 VNY(IX,IY,IZ,IEL)=VNY(IX,IY,IZ,IEL) + UNY(IR,IZ,IFC,IEL)
1089 ENDIF
1090 220 CONTINUE
1091 200 CONTINUE
1092C
1093 CALL DSSUM (VNX,lx1,ly1,lz1)
1094 CALL DSSUM (VNY,lx1,ly1,lz1)
1095 CALL UNITVEC (VNX,VNY,VNZ,NTOT1)
1096C
1097 CALL COPY (V1Y,VNX,NTOT1)
1098 CALL COPY (V1X,VNY,NTOT1)
1099 CALL CHSIGN (V1X,NTOT1)
1100C
1101 RETURN
1102 END
1103C-----------------------------------------------------------------------
1104 SUBROUTINE COMAVN3 (HVMASK,HFMASK,NEL)
1105C
1106 INCLUDE 'SIZE'
1107 INCLUDE 'GEOM'
1108 common /scrcg/ vnmag(lx1,ly1,lz1,lelt)
1109 common /indxfc/ mcrfc(4,6)
1110 $ , MFCCR(3,8)
1111 $ , MEGCR(3,8)
1112 $ , MFCEG(2,12)
1113 $ , MCREG(2,12)
1114 $ , MCRRST(3,8)
1115 $ , MIDRST(3,12)
1116 $ , MCRIND(8)
1117 $ , MEDIND(2,4)
1118 $ , NTEFC(2,12)
1119 $ , NTCRF(2,3)
1120C
1121 DIMENSION HVMASK(LX1,LY1,LZ1,1)
1122 $ , HFMASK(LX1,LZ1,6,1)
1123C
1124 NTOT1 = lx1*ly1*lz1*NEL
1125 NFACE = 2*ldim
1126 NCRNR = 2**ldim
1127 NEDGE = 12
1128 NMID =(lx1 + 1)/2
1129 NXM1 = lx1 - 1
1130 EPSA = 1.0E-06
1131 EPSN = 1.0E-03
1132C
1133 DO 100 IEL=1,NEL
1134 DO 100 ISD=1,NEDGE
1135 IX = MIDRST(1,ISD)
1136 IY = MIDRST(2,ISD)
1137 IZ = MIDRST(3,ISD)
1138 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1139 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 100
1140 CALL EDGINDV (LV1,LV2,LVSKIP,ISD)
1141 DO 120 I=2,NXM1
1142 LV = LV1 + (I-1)*LVSKIP
1143 VNX(LV,1,1,IEL) = 0.0
1144 VNY(LV,1,1,IEL) = 0.0
1145 VNZ(LV,1,1,IEL) = 0.0
1146 120 CONTINUE
1147 100 CONTINUE
1148C
1149 DO 150 IEL=1,NEL
1150 DO 150 ICR=1,NCRNR
1151 IX = MCRRST(1,ICR)
1152 IY = MCRRST(2,ICR)
1153 IZ = MCRRST(3,ICR)
1154 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1155 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 150
1156 VNX(IX,IY,IZ,IEL) = 0.0
1157 VNY(IX,IY,IZ,IEL) = 0.0
1158 VNZ(IX,IY,IZ,IEL) = 0.0
1159 150 CONTINUE
1160C
1161C (1) All Edges
1162C
1163 DO 200 IEL=1,NEL
1164 DO 200 ISD=1,NEDGE
1165 IX = MIDRST(1,ISD)
1166 IY = MIDRST(2,ISD)
1167 IZ = MIDRST(3,ISD)
1168 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1169 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 200
1170 DO 220 II=1,2
1171 IFC = MFCEG(II,ISD)
1172 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1173 IF (ABS(HMV - HMF) .LT. EPSA) THEN
1174 CALL EDGINDV (LV1,LV2,LVSKIP,ISD)
1175 CALL EDGINDF (LF1,LF2,LFSKIP,ISD,II)
1176 DO 240 I=2,NXM1
1177 LV = LV1 + (I-1)*LVSKIP
1178 LF = LF1 + (I-1)*LFSKIP
1179 VNX(LV,1,1,IEL)=VNX(LV,1,1,IEL)+UNX(LF,1,IFC,IEL)
1180 VNY(LV,1,1,IEL)=VNY(LV,1,1,IEL)+UNY(LF,1,IFC,IEL)
1181 VNZ(LV,1,1,IEL)=VNZ(LV,1,1,IEL)+UNZ(LF,1,IFC,IEL)
1182 240 CONTINUE
1183 ENDIF
1184 220 CONTINUE
1185 200 CONTINUE
1186C
1187C (2) All Corners
1188C
1189 DO 300 IEL=1,NEL
1190 DO 300 ICR=1,NCRNR
1191 IX = MCRRST(1,ICR)
1192 IY = MCRRST(2,ICR)
1193 IZ = MCRRST(3,ICR)
1194 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1195 IF (HMV .LT. 1.9 .OR. HMV .GT. 3.1) GOTO 300
1196 DO 320 II=1,3
1197 IFC = MFCCR(II,ICR)
1198 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1199 IF (ABS(HMV - HMF) .LT. EPSA) THEN
1200 IRA = NTCRF(1,II)
1201 ISA = NTCRF(2,II)
1202 IR = MCRRST(IRA,ICR)
1203 IS = MCRRST(ISA,ICR)
1204 VNX(IX,IY,IZ,IEL)=VNX(IX,IY,IZ,IEL)+UNX(IR,IS,IFC,IEL)
1205 VNY(IX,IY,IZ,IEL)=VNY(IX,IY,IZ,IEL)+UNY(IR,IS,IFC,IEL)
1206 VNZ(IX,IY,IZ,IEL)=VNZ(IX,IY,IZ,IEL)+UNZ(IR,IS,IFC,IEL)
1207 ENDIF
1208 320 CONTINUE
1209 300 CONTINUE
1210C
1211 CALL DSSUM (VNX,lx1,ly1,lz1)
1212 CALL DSSUM (VNY,lx1,ly1,lz1)
1213 CALL DSSUM (VNZ,lx1,ly1,lz1)
1214 CALL UNITVEC (VNX,VNY,VNZ,NTOT1)
1215 CALL VDOT3 (VNMAG,VNX,VNY,VNZ,VNX,VNY,VNZ,NTOT1)
1216C
1217 DO 500 IEL=1,NEL
1218 DO 500 IFC=1,NFACE
1219 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1220 IF (HMF .LT. 1.9 .OR. HMF .GT. 3.1) GOTO 500
1221 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
1222 DO 520 J2=JS2,JF2,JSKIP2
1223 DO 520 J1=JS1,JF1,JSKIP1
1224 IF (VNMAG(J1,J2,1,IEL) .LT. EPSA) GOTO 520
1225 VlzdIF = ABS(VNZ(J1,J2,1,IEL)) - 1.0
1226 IF (ABS(VlzdIF) .LT. EPSN) THEN
1227 V1X(J1,J2,1,IEL) = 1.0
1228 V1Y(J1,J2,1,IEL) = 0.0
1229 V1Z(J1,J2,1,IEL) = 0.0
1230 ELSE
1231 SSN = SQRT(VNX(J1,J2,1,IEL)**2 + VNY(J1,J2,1,IEL)**2)
1232 V1X(J1,J2,1,IEL) = -VNY(J1,J2,1,IEL) / SSN
1233 V1Y(J1,J2,1,IEL) = VNX(J1,J2,1,IEL) / SSN
1234 V1Z(J1,J2,1,IEL) = 0.0
1235 ENDIF
1236 520 CONTINUE
1237 500 CONTINUE
1238C
1239 CALL DSSUM (V1X,lx1,ly1,lz1)
1240 CALL DSSUM (V1Y,lx1,ly1,lz1)
1241 CALL DSSUM (V1Z,lx1,ly1,lz1)
1242 CALL UNITVEC (V1X,V1Y,V1Z,NTOT1)
1243C
1244 CALL VCROSS (V2X,V2Y,V2Z,VNX,VNY,VNZ,V1X,V1Y,V1Z,NTOT1)
1245C
1246 RETURN
1247 END
1248C-----------------------------------------------------------------------
1249 SUBROUTINE FIXWMSK (W2MASK,W3MASK,HVMASK,HFMASK,NEL)
1250C
1251 INCLUDE 'SIZE'
1252 DIMENSION HVMASK(LX1,LY1,LZ1,1)
1253 $ , HFMASK(LX1,LZ1,6,1)
1254C
1255 DIMENSION W2MASK(LX1,LY1,LZ1,1)
1256 $ , W3MASK(LX1,LY1,LZ1,1)
1257C
1258 IF (ldim.EQ.2) THEN
1259 CALL FXWMS2 (W2MASK,HVMASK,HFMASK,NEL)
1260 ELSE
1261 CALL FXWMS3 (W2MASK,W3MASK,HVMASK,HFMASK,NEL)
1262 ENDIF
1263C
1264 CALL DSOP(W2MASK,'MUL',lx1,ly1,lz1)
1265 IF (ldim.EQ.3) CALL DSOP(W3MASK,'MUL',lx1,ly1,lz1)
1266C
1267 RETURN
1268 END
1269C-----------------------------------------------------------------------
1270 SUBROUTINE FXWMS2 (W2MASK,HVMASK,HFMASK,NEL)
1271C
1272 INCLUDE 'SIZE'
1273 INCLUDE 'GEOM'
1274 common /indxfc/ mcrfc(4,6)
1275 $ , MFCCR(3,8)
1276 $ , MEGCR(3,8)
1277 $ , MFCEG(2,12)
1278 $ , MCREG(2,12)
1279 $ , MCRRST(3,8)
1280 $ , MIDRST(3,12)
1281 $ , MCRIND(8)
1282 $ , MEDIND(2,4)
1283 $ , NTEFC(2,12)
1284 $ , NTCRF(2,3)
1285C
1286 DIMENSION W2MASK(LX1,LY1,LZ1,1)
1287 $ , HVMASK(LX1,LY1,LZ1,1)
1288 $ , HFMASK(LX1,LZ1,6,1)
1289C
1290 NCRNR = 2**ldim
1291 EPSA = 1.0E-06
1292C
1293 IZ = 1
1294 DO 100 IEL=1,NEL
1295 DO 100 ICR=1,NCRNR
1296 IX = MCRRST(1,ICR)
1297 IY = MCRRST(2,ICR)
1298 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1299 IF (HMV .LT. 1.9 .OR. HMV .GT. 2.1) GOTO 100
1300 DO 120 II=1,2
1301 IFC = MFCCR(II,ICR)
1302 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1303 IF (ABS(HMV - HMF) .LT. EPSA) THEN
1304 IR = MCRRST(II,ICR)
1305 DOT = VNX(IX,IY,IZ,IEL)*UNX(IR,IZ,IFC,IEL) +
1306 $ VNY(IX,IY,IZ,IEL)*UNY(IR,IZ,IFC,IEL)
1307 IF (DOT .LT. 0.99) THEN
1308 W2MASK(IX,IY,IZ,IEL) = 0.0
1309 GOTO 100
1310 ENDIF
1311 ENDIF
1312 120 CONTINUE
1313 100 CONTINUE
1314C
1315 RETURN
1316 END
1317C-----------------------------------------------------------------------
1318 SUBROUTINE FXWMS3 (W2MASK,W3MASK,HVMASK,HFMASK,NEL)
1319C
1320 INCLUDE 'SIZE'
1321 INCLUDE 'GEOM'
1322 common /indxfc/ mcrfc(4,6)
1323 $ , MFCCR(3,8)
1324 $ , MEGCR(3,8)
1325 $ , MFCEG(2,12)
1326 $ , MCREG(2,12)
1327 $ , MCRRST(3,8)
1328 $ , MIDRST(3,12)
1329 $ , MCRIND(8)
1330 $ , MEDIND(2,4)
1331 $ , NTEFC(2,12)
1332 $ , NTCRF(2,3)
1333C
1334 DIMENSION W2MASK(LX1,LY1,LZ1,1)
1335 $ , W3MASK(LX1,LY1,LZ1,1)
1336 $ , HVMASK(LX1,LY1,LZ1,1)
1337 $ , HFMASK(LX1,LZ1,6,1)
1338C
1339 NCRNR = 2**ldim
1340 NEDGE = 12
1341 NMID = (lx1 + 1)/2
1342 EPSA = 1.0E-06
1343C
1344 DO 100 IEL=1,NEL
1345 DO 100 ISD=1,12
1346 IX = MIDRST(1,ISD)
1347 IY = MIDRST(2,ISD)
1348 IZ = MIDRST(3,ISD)
1349 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1350 IF (HMV .LT. 1.9 .OR. HMV .GT. 2.1) GOTO 100
1351 DO 120 II=1,2
1352 IFC = MFCEG(II,ISD)
1353 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1354 IF (ABS(HMV - HMF) .LT. EPSA) THEN
1355 CALL EDGINDF (LF1,LF2,LFSKIP,ISD,II)
1356 LF = LF1 + (NMID-1)*LFSKIP
1357 DOT = VNX(IX,IY,IZ,IEL)*UNX(LF,1,IFC,IEL) +
1358 $ VNY(IX,IY,IZ,IEL)*UNY(LF,1,IFC,IEL) +
1359 $ VNZ(IX,IY,IZ,IEL)*UNZ(LF,1,IFC,IEL)
1360 IF (DOT .LT. 0.99) THEN
1361 CALL EDGINDV (LV1,LV2,LVSKIP,ISD)
1362 DO 140 LV=LV1,LV2,LVSKIP
1363 W3MASK(LV,1,1,IEL) = 0.0
1364 140 CONTINUE
1365 ENDIF
1366 ENDIF
1367 120 CONTINUE
1368 100 CONTINUE
1369C
1370C All Corners
1371C
1372 DO 300 IEL=1,NEL
1373 DO 300 ICR=1,NCRNR
1374 IX = MCRRST(1,ICR)
1375 IY = MCRRST(2,ICR)
1376 IZ = MCRRST(3,ICR)
1377 HMV = ABS( HVMASK(IX,IY,IZ,IEL) )
1378 IF (HMV .LT. 1.9 .OR. HMV .GT. 2.1) GOTO 300
1379 DO 320 II=1,3
1380 IFC = MFCCR(II,ICR)
1381 HMF = ABS( HFMASK(1,1,IFC,IEL) )
1382 IF (ABS(HMV - HMF) .LT. EPSA) THEN
1383 IRA = NTCRF(1,II)
1384 ISA = NTCRF(2,II)
1385 IR = MCRRST(IRA,ICR)
1386 IS = MCRRST(ISA,ICR)
1387 DOT = VNX(IX,IY,IZ,IEL)*UNX(IR,IS,IFC,IEL) +
1388 $ VNY(IX,IY,IZ,IEL)*UNY(IR,IS,IFC,IEL) +
1389 $ VNZ(IX,IY,IZ,IEL)*UNZ(IR,IS,IFC,IEL)
1390 ENDIF
1391 IF (DOT .LT. 0.99) THEN
1392 W2MASK(IX,IY,IZ,IEL) = 0.0
1393 GOTO 300
1394 ENDIF
1395 320 CONTINUE
1396 300 CONTINUE
1397C
1398 RETURN
1399 END
1400C-----------------------------------------------------------------------
1401 SUBROUTINE SETCDAT
1402C
1403 INCLUDE 'SIZE'
1404 common /indxfc/ mcrfc(4,6)
1405 $ , MFCCR(3,8)
1406 $ , MEGCR(3,8)
1407 $ , MFCEG(2,12)
1408 $ , MCREG(2,12)
1409 $ , MCRRST(3,8)
1410 $ , MIDRST(3,12)
1411 $ , MCRIND(8)
1412 $ , MEDIND(2,4)
1413 $ , NTEFC(2,12)
1414 $ , NTCRF(2,3)
1415C
1416 NMID = (lx1 +1)/2
1417C
1418C Corners on faces
1419C
1420 MCRFC(1,1) = 1
1421 MCRFC(2,1) = 2
1422 MCRFC(3,1) = 6
1423 MCRFC(4,1) = 5
1424 MCRFC(1,2) = 2
1425 MCRFC(2,2) = 3
1426 MCRFC(3,2) = 7
1427 MCRFC(4,2) = 6
1428 MCRFC(1,3) = 3
1429 MCRFC(2,3) = 4
1430 MCRFC(3,3) = 8
1431 MCRFC(4,3) = 7
1432 MCRFC(1,4) = 4
1433 MCRFC(2,4) = 1
1434 MCRFC(3,4) = 5
1435 MCRFC(4,4) = 8
1436 MCRFC(1,5) = 4
1437 MCRFC(2,5) = 3
1438 MCRFC(3,5) = 2
1439 MCRFC(4,5) = 1
1440 MCRFC(1,6) = 5
1441 MCRFC(2,6) = 6
1442 MCRFC(3,6) = 7
1443 MCRFC(4,6) = 8
1444C
1445C Faces at corners
1446C
1447 MFCCR(1,1) = 4
1448 MFCCR(2,1) = 1
1449 MFCCR(3,1) = 5
1450 MFCCR(1,2) = 1
1451 MFCCR(2,2) = 2
1452 MFCCR(3,2) = 5
1453 MFCCR(1,3) = 2
1454 MFCCR(2,3) = 3
1455 MFCCR(3,3) = 5
1456 MFCCR(1,4) = 3
1457 MFCCR(2,4) = 4
1458 MFCCR(3,4) = 5
1459 MFCCR(1,5) = 4
1460 MFCCR(2,5) = 1
1461 MFCCR(3,5) = 6
1462 MFCCR(1,6) = 1
1463 MFCCR(2,6) = 2
1464 MFCCR(3,6) = 6
1465 MFCCR(1,7) = 2
1466 MFCCR(2,7) = 3
1467 MFCCR(3,7) = 6
1468 MFCCR(1,8) = 3
1469 MFCCR(2,8) = 4
1470 MFCCR(3,8) = 6
1471C
1472C Edges at corners
1473C
1474 MEGCR(1,1) = 4
1475 MEGCR(2,1) = 1
1476 MEGCR(3,1) = 9
1477 MEGCR(1,2) = 1
1478 MEGCR(2,2) = 2
1479 MEGCR(3,2) = 10
1480 MEGCR(1,3) = 2
1481 MEGCR(2,3) = 3
1482 MEGCR(3,3) = 11
1483 MEGCR(1,4) = 3
1484 MEGCR(2,4) = 4
1485 MEGCR(3,4) = 12
1486 MEGCR(1,5) = 8
1487 MEGCR(2,5) = 5
1488 MEGCR(3,5) = 9
1489 MEGCR(1,6) = 5
1490 MEGCR(2,6) = 6
1491 MEGCR(3,6) = 10
1492 MEGCR(1,7) = 6
1493 MEGCR(2,7) = 7
1494 MEGCR(3,7) = 11
1495 MEGCR(1,8) = 7
1496 MEGCR(2,8) = 8
1497 MEGCR(3,8) = 12
1498C
1499C Faces on edges
1500C
1501 MFCEG(1,1) = 1
1502 MFCEG(2,1) = 5
1503 MFCEG(1,2) = 2
1504 MFCEG(2,2) = 5
1505 MFCEG(1,3) = 3
1506 MFCEG(2,3) = 5
1507 MFCEG(1,4) = 4
1508 MFCEG(2,4) = 5
1509 MFCEG(1,5) = 1
1510 MFCEG(2,5) = 6
1511 MFCEG(1,6) = 2
1512 MFCEG(2,6) = 6
1513 MFCEG(1,7) = 3
1514 MFCEG(2,7) = 6
1515 MFCEG(1,8) = 4
1516 MFCEG(2,8) = 6
1517 MFCEG(1,9) = 4
1518 MFCEG(2,9) = 1
1519 MFCEG(1,10) = 1
1520 MFCEG(2,10) = 2
1521 MFCEG(1,11) = 2
1522 MFCEG(2,11) = 3
1523 MFCEG(1,12) = 3
1524 MFCEG(2,12) = 4
1525C
1526C Corners at edges
1527C
1528 MCREG(1,1) = 1
1529 MCREG(2,1) = 2
1530 MCREG(1,2) = 2
1531 MCREG(2,2) = 3
1532 MCREG(1,3) = 4
1533 MCREG(2,3) = 3
1534 MCREG(1,4) = 1
1535 MCREG(2,4) = 4
1536 MCREG(1,5) = 5
1537 MCREG(2,5) = 6
1538 MCREG(1,6) = 6
1539 MCREG(2,6) = 7
1540 MCREG(1,7) = 8
1541 MCREG(2,7) = 7
1542 MCREG(1,8) = 5
1543 MCREG(2,8) = 8
1544 MCREG(1,9) = 1
1545 MCREG(2,9) = 5
1546 MCREG(1,10) = 2
1547 MCREG(2,10) = 6
1548 MCREG(1,11) = 3
1549 MCREG(2,11) = 7
1550 MCREG(1,12) = 4
1551 MCREG(2,12) = 8
1552C
1553C Corner indices (Vol array)
1554C
1555 MCRRST(1,1) = 1
1556 MCRRST(2,1) = 1
1557 MCRRST(3,1) = 1
1558 MCRRST(1,2) = lx1
1559 MCRRST(2,2) = 1
1560 MCRRST(3,2) = 1
1561 MCRRST(1,3) = lx1
1562 MCRRST(2,3) = lx1
1563 MCRRST(3,3) = 1
1564 MCRRST(1,4) = 1
1565 MCRRST(2,4) = lx1
1566 MCRRST(3,4) = 1
1567 MCRRST(1,5) = 1
1568 MCRRST(2,5) = 1
1569 MCRRST(3,5) = lx1
1570 MCRRST(1,6) = lx1
1571 MCRRST(2,6) = 1
1572 MCRRST(3,6) = lx1
1573 MCRRST(1,7) = lx1
1574 MCRRST(2,7) = lx1
1575 MCRRST(3,7) = lx1
1576 MCRRST(1,8) = 1
1577 MCRRST(2,8) = lx1
1578 MCRRST(3,8) = lx1
1579C
1580C Mid-edge indcies (Vol array)
1581C
1582 MIDRST(1,1) = NMID
1583 MIDRST(1,2) = lx1
1584 MIDRST(1,3) = NMID
1585 MIDRST(1,4) = 1
1586 MIDRST(1,5) = NMID
1587 MIDRST(1,6) = lx1
1588 MIDRST(1,7) = NMID
1589 MIDRST(1,8) = 1
1590 MIDRST(1,9) = 1
1591 MIDRST(1,10) = lx1
1592 MIDRST(1,11) = lx1
1593 MIDRST(1,12) = 1
1594 MIDRST(2,1) = 1
1595 MIDRST(2,2) = NMID
1596 MIDRST(2,3) = lx1
1597 MIDRST(2,4) = NMID
1598 MIDRST(2,5) = 1
1599 MIDRST(2,6) = NMID
1600 MIDRST(2,7) = lx1
1601 MIDRST(2,8) = NMID
1602 MIDRST(2,9) = 1
1603 MIDRST(2,10) = 1
1604 MIDRST(2,11) = lx1
1605 MIDRST(2,12) = lx1
1606 MIDRST(3,1) = 1
1607 MIDRST(3,2) = 1
1608 MIDRST(3,3) = 1
1609 MIDRST(3,4) = 1
1610 MIDRST(3,5) = lx1
1611 MIDRST(3,6) = lx1
1612 MIDRST(3,7) = lx1
1613 MIDRST(3,8) = lx1
1614 MIDRST(3,9) = NMID
1615 MIDRST(3,10) = NMID
1616 MIDRST(3,11) = NMID
1617 MIDRST(3,12) = NMID
1618C
1619C 1-D corners indices (Vol array)
1620C
1621 MCRIND(1) = 1
1622 MCRIND(2) = lx1
1623 MCRIND(3) = lx1**2
1624 MCRIND(7) = lx1**3
1625 MCRIND(4) = MCRIND(3) - lx1 + 1
1626 MCRIND(5) = MCRIND(7) - MCRIND(3) + 1
1627 MCRIND(6) = MCRIND(5) + lx1 - 1
1628 MCRIND(8) = MCRIND(7) - lx1 + 1
1629C
1630C 1-D edge indices (Face array)
1631C
1632 MEDIND(1,1) = 1
1633 MEDIND(2,1) = lx1
1634 MEDIND(1,2) = lx1**2 - lx1 + 1
1635 MEDIND(2,2) = lx1**2
1636 MEDIND(1,3) = 1
1637 MEDIND(2,3) = MEDIND(1,2)
1638 MEDIND(1,4) = lx1
1639 MEDIND(2,4) = lx1**2
1640C
1641C 1-D edge index type (Face array)
1642C
1643 NTEFC(1,1) = 1
1644 NTEFC(2,1) = 1
1645 NTEFC(1,2) = 1
1646 NTEFC(2,2) = 4
1647 NTEFC(1,3) = 1
1648 NTEFC(2,3) = 2
1649 NTEFC(1,4) = 1
1650 NTEFC(2,4) = 3
1651 NTEFC(1,5) = 2
1652 NTEFC(2,5) = 1
1653 NTEFC(1,6) = 2
1654 NTEFC(2,6) = 4
1655 NTEFC(1,7) = 2
1656 NTEFC(2,7) = 2
1657 NTEFC(1,8) = 2
1658 NTEFC(2,8) = 3
1659 NTEFC(1,9) = 3
1660 NTEFC(2,9) = 3
1661 NTEFC(1,10) = 4
1662 NTEFC(2,10) = 3
1663 NTEFC(1,11) = 4
1664 NTEFC(2,11) = 4
1665 NTEFC(1,12) = 3
1666 NTEFC(2,12) = 4
1667C
1668C Corner index address on face in MCRRST
1669C
1670 NTCRF(1,1) = 1
1671 NTCRF(2,1) = 3
1672 NTCRF(1,2) = 2
1673 NTCRF(2,2) = 3
1674 NTCRF(1,3) = 1
1675 NTCRF(2,3) = 2
1676C
1677 RETURN
1678 END
1679C-----------------------------------------------------------------------
1680 SUBROUTINE EDGINDF (LF1,LF2,LFSKIP,ISD,IFCN)
1681C
1682 INCLUDE 'SIZE'
1683 common /indxfc/ mcrfc(4,6)
1684 $ , MFCCR(3,8)
1685 $ , MEGCR(3,8)
1686 $ , MFCEG(2,12)
1687 $ , MCREG(2,12)
1688 $ , MCRRST(3,8)
1689 $ , MIDRST(3,12)
1690 $ , MCRIND(8)
1691 $ , MEDIND(2,4)
1692 $ , NTEFC(2,12)
1693 $ , NTCRF(2,3)
1694C
1695 ITYP = NTEFC(IFCN,ISD)
1696C
1697 LF1 = MEDIND(1,ITYP)
1698 LF2 = MEDIND(2,ITYP)
1699C
1700 LFSKIP = 1
1701 IF (ITYP .GE. 3) LFSKIP = lx1
1702C
1703 RETURN
1704 END
1705C-----------------------------------------------------------------------
1706 SUBROUTINE EDGINDV (LV1,LV2,LVSKIP,ISD)
1707C
1708 INCLUDE 'SIZE'
1709 common /indxfc/ mcrfc(4,6)
1710 $ , mfccr(3,8)
1711 $ , MEGCR(3,8)
1712 $ , MFCEG(2,12)
1713 $ , MCREG(2,12)
1714 $ , MCRRST(3,8)
1715 $ , MIDRST(3,12)
1716 $ , MCRIND(8)
1717 $ , MEDIND(2,4)
1718 $ , NTEFC(2,12)
1719 $ , NTCRF(2,3)
1720C
1721 IODD = ISD - ISD/2*2
1722 ICR1 = MCREG(1,ISD)
1723 ICR2 = MCREG(2,ISD)
1724C
1725 LV1 = MCRIND(ICR1)
1726 LV2 = MCRIND(ICR2)
1727C
1728 IF (ISD .GE. 9) THEN
1729 LVSKIP = lx1**2
1730 ELSE
1731 IF (IODD.EQ.0) THEN
1732 LVSKIP = lx1
1733 ELSE
1734 LVSKIP = 1
1735 ENDIF
1736 ENDIF
1737C
1738 RETURN
1739 END
1740C-----------------------------------------------------------------------
1741 SUBROUTINE SETCDOF
1742C
1743 INCLUDE 'SIZE'
1744 INCLUDE 'INPUT'
1745C
1746 NFACE = 2*ldim
1747C
1748 DO 100 IEL=1,NELT
1749 DO 100 IFC=1,NFACE
1750 CDOF(IFC,IEL)=CBC(IFC,IEL,0)(1:1)
1751 100 CONTINUE
1752C
1753 RETURN
1754 END
1755C-----------------------------------------------------------------------
1756 SUBROUTINE AMASK (VB1,VB2,VB3,V1,V2,V3,NEL)
1757C
1758 INCLUDE 'SIZE'
1759 INCLUDE 'GEOM'
1760 INCLUDE 'INPUT'
1761 INCLUDE 'SOLN'
1762 INCLUDE 'TSTEP'
1763 common /scrsf/ a1mask(lx1,ly1,lz1,lelt)
1764 $ , a2mask(lx1,ly1,lz1,lelt)
1765 $ , a3mask(lx1,ly1,lz1,lelt)
1766 common /ctmp0/ wa(lx1,ly1,lz1,lelt)
1767C
1768 DIMENSION VB1(LX1,LY1,LZ1,1)
1769 $ , VB2(LX1,LY1,LZ1,1)
1770 $ , VB3(LX1,LY1,LZ1,1)
1771 $ , V1(LX1,LY1,LZ1,1)
1772 $ , V2(LX1,LY1,LZ1,1)
1773 $ , V3(LX1,LY1,LZ1,1)
1774C
1775 NTOT1 = lx1*ly1*lz1*NEL
1776 CALL RONE (WA,NTOT1)
1777 CALL COPY (VB1,V1,NTOT1)
1778 CALL COPY (VB2,V2,NTOT1)
1779 IF (ldim.EQ.3) CALL COPY (VB3,V3,NTOT1)
1780C
1781 IF (IFIELD.EQ.1) THEN
1782 CALL SUB3 (A1MASK,WA,V1MASK,NTOT1)
1783 CALL SUB3 (A2MASK,WA,V2MASK,NTOT1)
1784 IF (ldim.EQ.3) CALL SUB3 (A3MASK,WA,V3MASK,NTOT1)
1785 ELSEIF (IFIELD.EQ.ifldmhd) THEN
1786 CALL SUB3 (A1MASK,WA,B1MASK,NTOT1)
1787 CALL SUB3 (A2MASK,WA,B2MASK,NTOT1)
1788 IF (ldim.EQ.3) CALL SUB3 (A3MASK,WA,B3MASK,NTOT1)
1789 ELSE
1790 CALL SUB3 (A1MASK,WA,W1MASK,NTOT1)
1791 CALL SUB3 (A2MASK,WA,W2MASK,NTOT1)
1792 IF (ldim.EQ.3) CALL SUB3 (A3MASK,WA,W3MASK,NTOT1)
1793 ENDIF
1794C
1795 CALL QMASK (VB1,VB2,VB3,A1MASK,A2MASK,A3MASK,NEL)
1796C
1797 RETURN
1798 END
1799C-----------------------------------------------------------------------
1800 SUBROUTINE RMASK (R1,R2,R3,NEL)
1801C
1802 INCLUDE 'SIZE'
1803 INCLUDE 'INPUT'
1804 INCLUDE 'SOLN'
1805 INCLUDE 'TSTEP'
1806 INCLUDE 'MVGEOM'
1807C
1808 DIMENSION R1 (LX1,LY1,LZ1,1)
1809 $ , R2 (LX1,LY1,LZ1,1)
1810 $ , R3 (LX1,LY1,LZ1,1)
1811C
1812c if (ifsplit) then
1813c call opmask(r1,r2,r3)
1814c return
1815c endif
1816
1817c call outfldro (v1mask,'v1mask rmk',0)
1818c call outfldro (v2mask,'v2mask rmk',1)
1819
1820 IF (IFIELD.EQ.1) THEN
1821 CALL QMASK (R1,R2,R3,V1MASK,V2MASK,V3MASK,NEL)
1822 ELSEIF (ifield.eq.ifldmhd) then
1823 CALL QMASK (R1,R2,R3,B1MASK,B2MASK,B3MASK,NEL)
1824 ELSE
1825 CALL QMASK (R1,R2,R3,W1MASK,W2MASK,W3MASK,NEL)
1826 ENDIF
1827
1828c call outfldro (r1,'r1 rmask',0)
1829c call outfldro (r2,'r2 rmask',1)
1830c call exitti ('quit in rmask$,',nel)
1831
1832 RETURN
1833 END
1834C-----------------------------------------------------------------------
1835 SUBROUTINE QMASK (R1,R2,R3,R1MASK,R2MASK,R3MASK,NEL)
1836C
1837 INCLUDE 'SIZE'
1838 INCLUDE 'GEOM'
1839 INCLUDE 'TSTEP'
1840 common /ctmp1/ s1(lx1,ly1,lz1,lelt)
1841 $ , s2(lx1,ly1,lz1,lelt)
1842 $ , s3(lx1,ly1,lz1,lelt)
1843C
1844 DIMENSION R1(LX1,LY1,LZ1,1)
1845 $ , R2(LX1,LY1,LZ1,1)
1846 $ , R3(LX1,LY1,LZ1,1)
1847 $ , R1MASK(LX1,LY1,LZ1,1)
1848 $ , R2MASK(LX1,LY1,LZ1,1)
1849 $ , R3MASK(LX1,LY1,LZ1,1)
1850C
1851 NTOT1 = lx1*ly1*lz1*NEL
1852C
1853C (0) Collocate Volume Mask
1854C
1855 CALL COPY (S1,R1,NTOT1)
1856 CALL COPY (S2,R2,NTOT1)
1857 CALL COL2 (R1,R1MASK,NTOT1)
1858 CALL COL2 (R2,R2MASK,NTOT1)
1859 IF (ldim.EQ.3) THEN
1860 CALL COPY (S3,R3,NTOT1)
1861 CALL COL2 (R3,R3MASK,NTOT1)
1862 ENDIF
1863C
1864C (1) Face Mask
1865C
1866 IF (IFLMSF(IFIELD)) THEN
1867 IF (ldim.EQ.2) THEN
1868 CALL FCMSK2 (R1,R2,S1,S2,R1MASK,R2MASK,NEL)
1869 ELSE
1870 CALL FCMSK3 (R1,R2,R3,S1,S2,S3,R1MASK,R2MASK,R3MASK,NEL)
1871 ENDIF
1872 ENDIF
1873C
1874C (2) Edge Mask (3-D only)
1875C
1876 IF (ldim.EQ.3 .AND. IFLMSE(IFIELD))
1877 $ CALL EGMASK (R1,R2,R3,S1,S2,S3,R1MASK,R2MASK,R3MASK,NEL)
1878C
1879C (3) Corner Mask
1880C
1881 IF (IFLMSC(IFIELD)) THEN
1882 IF (ldim.EQ.2) THEN
1883 CALL CRMSK2 (R1,R2,S1,S2,R1MASK,R2MASK,NEL)
1884 ELSE
1885 CALL CRMSK3 (R1,R2,R3,S1,S2,S3,R1MASK,R2MASK,R3MASK,NEL)
1886 ENDIF
1887 ENDIF
1888C
1889 RETURN
1890 END
1891C-----------------------------------------------------------------------
1892 SUBROUTINE FCMSK2 (R1,R2,S1,S2,R1MASK,R2MASK,NEL)
1893C
1894 INCLUDE 'SIZE'
1895 INCLUDE 'GEOM'
1896 INCLUDE 'TSTEP'
1897 DIMENSION R1(LX1,LY1,LZ1,1)
1898 $ , R2(LX1,LY1,LZ1,1)
1899 $ , S1(LX1,LY1,LZ1,1)
1900 $ , S2(LX1,LY1,LZ1,1)
1901 $ , R1MASK(LX1,LY1,LZ1,1)
1902 $ , R2MASK(LX1,LY1,LZ1,1)
1903C
1904 NFACE = 2*ldim
1905C
1906 DO 100 IEL=1,NEL
1907 DO 100 IFC=1,NFACE
1908 IF (.NOT.IFMSFC(IFC,IEL,IFIELD)) GO TO 100
1909 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
1910 DO 120 J2=JS2,JF2,JSKIP2
1911 DO 120 J1=JS1,JF1,JSKIP1
1912 RNOR = ( S1(J1,J2,1,IEL)*VNX(J1,J2,1,IEL) +
1913 $ S2(J1,J2,1,IEL)*VNY(J1,J2,1,IEL) ) *
1914 $ R1MASK(J1,J2,1,IEL)
1915 RTN1 = ( S1(J1,J2,1,IEL)*V1X(J1,J2,1,IEL) +
1916 $ S2(J1,J2,1,IEL)*V1Y(J1,J2,1,IEL) ) *
1917 $ R2MASK(J1,J2,1,IEL)
1918 R1(J1,J2,1,IEL) = RNOR*VNX(J1,J2,1,IEL) +
1919 $ RTN1*V1X(J1,J2,1,IEL)
1920 R2(J1,J2,1,IEL) = RNOR*VNY(J1,J2,1,IEL) +
1921 $ RTN1*V1Y(J1,J2,1,IEL)
1922 120 CONTINUE
1923 100 CONTINUE
1924C
1925 RETURN
1926 END
1927C-----------------------------------------------------------------------
1928 SUBROUTINE FCMSK3 (R1,R2,R3,S1,S2,S3,R1MASK,R2MASK,R3MASK,NEL)
1929C
1930 INCLUDE 'SIZE'
1931 INCLUDE 'GEOM'
1932 INCLUDE 'TSTEP'
1933 DIMENSION R1(LX1,LY1,LZ1,1)
1934 $ , R2(LX1,LY1,LZ1,1)
1935 $ , R3(LX1,LY1,LZ1,1)
1936 $ , S1(LX1,LY1,LZ1,1)
1937 $ , S2(LX1,LY1,LZ1,1)
1938 $ , S3(LX1,LY1,LZ1,1)
1939 $ , R1MASK(LX1,LY1,LZ1,1)
1940 $ , R2MASK(LX1,LY1,LZ1,1)
1941 $ , R3MASK(LX1,LY1,LZ1,1)
1942C
1943 NFACE = 2*ldim
1944C
1945 DO 100 IEL=1,NEL
1946 DO 100 IFC=1,NFACE
1947 IF (.NOT.IFMSFC(IFC,IEL,IFIELD)) GO TO 100
1948 CALL FACIND2 (JS1,JF1,JSKIP1,JS2,JF2,JSKIP2,IFC)
1949 DO 120 J2=JS2,JF2,JSKIP2
1950 DO 120 J1=JS1,JF1,JSKIP1
1951 RNOR = ( S1(J1,J2,1,IEL)*VNX(J1,J2,1,IEL) +
1952 $ S2(J1,J2,1,IEL)*VNY(J1,J2,1,IEL) +
1953 $ S3(J1,J2,1,IEL)*VNZ(J1,J2,1,IEL) ) *
1954 $ R1MASK(J1,J2,1,IEL)
1955 RTN1 = ( S1(J1,J2,1,IEL)*V1X(J1,J2,1,IEL) +
1956 $ S2(J1,J2,1,IEL)*V1Y(J1,J2,1,IEL) +
1957 $ S3(J1,J2,1,IEL)*V1Z(J1,J2,1,IEL) ) *
1958 $ R2MASK(J1,J2,1,IEL)
1959 RTN2 = ( S1(J1,J2,1,IEL)*V2X(J1,J2,1,IEL) +
1960 $ S2(J1,J2,1,IEL)*V2Y(J1,J2,1,IEL) +
1961 $ S3(J1,J2,1,IEL)*V2Z(J1,J2,1,IEL) ) *
1962 $ R3MASK(J1,J2,1,IEL)
1963 R1(J1,J2,1,IEL) = RNOR*VNX(J1,J2,1,IEL) +
1964 $ RTN1*V1X(J1,J2,1,IEL) +
1965 $ RTN2*V2X(J1,J2,1,IEL)
1966 R2(J1,J2,1,IEL) = RNOR*VNY(J1,J2,1,IEL) +
1967 $ RTN1*V1Y(J1,J2,1,IEL) +
1968 $ RTN2*V2Y(J1,J2,1,IEL)
1969 R3(J1,J2,1,IEL) = RNOR*VNZ(J1,J2,1,IEL) +
1970 $ RTN1*V1Z(J1,J2,1,IEL) +
1971 $ RTN2*V2Z(J1,J2,1,IEL)
1972 120 CONTINUE
1973 100 CONTINUE
1974C
1975 RETURN
1976 END
1977C-----------------------------------------------------------------------
1978 SUBROUTINE EGMASK (R1,R2,R3,S1,S2,S3,R1MASK,R2MASK,R3MASK,NEL)
1979C
1980 INCLUDE 'SIZE'
1981 INCLUDE 'GEOM'
1982 INCLUDE 'TSTEP'
1983 DIMENSION R1(LX1,LY1,LZ1,1)
1984 $ , R2(LX1,LY1,LZ1,1)
1985 $ , R3(LX1,LY1,LZ1,1)
1986 $ , S1(LX1,LY1,LZ1,1)
1987 $ , S2(LX1,LY1,LZ1,1)
1988 $ , S3(LX1,LY1,LZ1,1)
1989 $ , R1MASK(LX1,LY1,LZ1,1)
1990 $ , R2MASK(LX1,LY1,LZ1,1)
1991 $ , R3MASK(LX1,LY1,LZ1,1)
1992C
1993 NEDGE = 12
1994C
1995 DO 100 IEL=1,NEL
1996 DO 100 ISD=1,NEDGE
1997 IF (.NOT.IFMSEG(ISD,IEL,IFIELD)) GOTO 100
1998 CALL EDGINDV (LV1,LV2,LVSKIP,ISD)
1999 DO 120 LV=LV1,LV2,LVSKIP
2000 RNOR = ( S1(LV,1,1,IEL)*VNX(LV,1,1,IEL) +
2001 $ S2(LV,1,1,IEL)*VNY(LV,1,1,IEL) +
2002 $ S3(LV,1,1,IEL)*VNZ(LV,1,1,IEL) ) *
2003 $ R1MASK(LV,1,1,IEL)
2004 RTN1 = ( S1(LV,1,1,IEL)*V1X(LV,1,1,IEL) +
2005 $ S2(LV,1,1,IEL)*V1Y(LV,1,1,IEL) +
2006 $ S3(LV,1,1,IEL)*V1Z(LV,1,1,IEL) ) *
2007 $ R2MASK(LV,1,1,IEL)
2008 RTN2 = ( S1(LV,1,1,IEL)*V2X(LV,1,1,IEL) +
2009 $ S2(LV,1,1,IEL)*V2Y(LV,1,1,IEL) +
2010 $ S3(LV,1,1,IEL)*V2Z(LV,1,1,IEL) ) *
2011 $ R3MASK(LV,1,1,IEL)
2012 R1(LV,1,1,IEL) = RNOR*VNX(LV,1,1,IEL) +
2013 $ RTN1*V1X(LV,1,1,IEL) +
2014 $ RTN2*V2X(LV,1,1,IEL)
2015 R2(LV,1,1,IEL) = RNOR*VNY(LV,1,1,IEL) +
2016 $ RTN1*V1Y(LV,1,1,IEL) +
2017 $ RTN2*V2Y(LV,1,1,IEL)
2018 R3(LV,1,1,IEL) = RNOR*VNZ(LV,1,1,IEL) +
2019 $ RTN1*V1Z(LV,1,1,IEL) +
2020 $ RTN2*V2Z(LV,1,1,IEL)
2021 120 CONTINUE
2022 100 CONTINUE
2023C
2024 RETURN
2025 END
2026C-----------------------------------------------------------------------
2027 SUBROUTINE CRMSK2 (R1,R2,S1,S2,R1MASK,R2MASK,NEL)
2028C
2029 INCLUDE 'SIZE'
2030 INCLUDE 'GEOM'
2031 INCLUDE 'TSTEP'
2032 common /indxfc/ mcrfc(4,6)
2033 $ , MFCCR(3,8)
2034 $ , MEGCR(3,8)
2035 $ , MFCEG(2,12)
2036 $ , MCREG(2,12)
2037 $ , MCRRST(3,8)
2038 $ , MIDRST(3,12)
2039 $ , MCRIND(8)
2040 $ , MEDIND(2,4)
2041 $ , NTEFC(2,12)
2042 $ , NTCRF(2,3)
2043 DIMENSION R1(LX1,LY1,LZ1,1)
2044 $ , R2(LX1,LY1,LZ1,1)
2045 $ , S1(LX1,LY1,LZ1,1)
2046 $ , S2(LX1,LY1,LZ1,1)
2047 $ , R1MASK(LX1,LY1,LZ1,1)
2048 $ , R2MASK(LX1,LY1,LZ1,1)
2049C
2050 NCRNR = 2**ldim
2051C
2052 DO 100 IEL=1,NEL
2053 DO 100 ICR=1,NCRNR
2054 IF (.NOT.IFMSCR(ICR,IEL,IFIELD)) GO TO 100
2055 IX = MCRRST(1,ICR)
2056 IY = MCRRST(2,ICR)
2057 IZ = 1
2058 RNOR = ( S1(IX,IY,IZ,IEL)*VNX(IX,IY,IZ,IEL) +
2059 $ S2(IX,IY,IZ,IEL)*VNY(IX,IY,IZ,IEL) ) *
2060 $ R1MASK(IX,IY,IZ,IEL)
2061 RTN1 = ( S1(IX,IY,IZ,IEL)*V1X(IX,IY,IZ,IEL) +
2062 $ S2(IX,IY,IZ,IEL)*V1Y(IX,IY,IZ,IEL) ) *
2063 $ R2MASK(IX,IY,IZ,IEL)
2064 R1(IX,IY,IZ,IEL) = RNOR*VNX(IX,IY,IZ,IEL) +
2065 $ RTN1*V1X(IX,IY,IZ,IEL)
2066 R2(IX,IY,IZ,IEL) = RNOR*VNY(IX,IY,IZ,IEL) +
2067 $ RTN1*V1Y(IX,IY,IZ,IEL)
2068 100 CONTINUE
2069C
2070 RETURN
2071 END
2072C-----------------------------------------------------------------------
2073 SUBROUTINE CRMSK3 (R1,R2,R3,S1,S2,S3,R1MASK,R2MASK,R3MASK,NEL)
2074C
2075 INCLUDE 'SIZE'
2076 INCLUDE 'GEOM'
2077 INCLUDE 'TSTEP'
2078 common /indxfc/ mcrfc(4,6)
2079 $ , MFCCR(3,8)
2080 $ , MEGCR(3,8)
2081 $ , MFCEG(2,12)
2082 $ , MCREG(2,12)
2083 $ , MCRRST(3,8)
2084 $ , MIDRST(3,12)
2085 $ , MCRIND(8)
2086 $ , MEDIND(2,4)
2087 $ , NTEFC(2,12)
2088 $ , NTCRF(2,3)
2089 DIMENSION R1(LX1,LY1,LZ1,1)
2090 $ , R2(LX1,LY1,LZ1,1)
2091 $ , R3(LX1,LY1,LZ1,1)
2092 $ , S1(LX1,LY1,LZ1,1)
2093 $ , S2(LX1,LY1,LZ1,1)
2094 $ , S3(LX1,LY1,LZ1,1)
2095 $ , R1MASK(LX1,LY1,LZ1,1)
2096 $ , R2MASK(LX1,LY1,LZ1,1)
2097 $ , R3MASK(LX1,LY1,LZ1,1)
2098C
2099 NCRNR = 2**ldim
2100C
2101 DO 100 IEL=1,NEL
2102 DO 100 ICR=1,NCRNR
2103 IF (.NOT.IFMSCR(ICR,IEL,IFIELD)) GO TO 100
2104 IX = MCRRST(1,ICR)
2105 IY = MCRRST(2,ICR)
2106 IZ = MCRRST(3,ICR)
2107 RNOR = ( S1(IX,IY,IZ,IEL)*VNX(IX,IY,IZ,IEL) +
2108 $ S2(IX,IY,IZ,IEL)*VNY(IX,IY,IZ,IEL) +
2109 $ S3(IX,IY,IZ,IEL)*VNZ(IX,IY,IZ,IEL) ) *
2110 $ R1MASK(IX,IY,IZ,IEL)
2111 RTN1 = ( S1(IX,IY,IZ,IEL)*V1X(IX,IY,IZ,IEL) +
2112 $ S2(IX,IY,IZ,IEL)*V1Y(IX,IY,IZ,IEL) +
2113 $ S3(IX,IY,IZ,IEL)*V1Z(IX,IY,IZ,IEL) ) *
2114 $ R2MASK(IX,IY,IZ,IEL)
2115 RTN2 = ( S1(IX,IY,IZ,IEL)*V2X(IX,IY,IZ,IEL) +
2116 $ S2(IX,IY,IZ,IEL)*V2Y(IX,IY,IZ,IEL) +
2117 $ S3(IX,IY,IZ,IEL)*V2Z(IX,IY,IZ,IEL) ) *
2118 $ R3MASK(IX,IY,IZ,IEL)
2119 R1(IX,IY,IZ,IEL) = RNOR*VNX(IX,IY,IZ,IEL) +
2120 $ RTN1*V1X(IX,IY,IZ,IEL) +
2121 $ RTN2*V2X(IX,IY,IZ,IEL)
2122 R2(IX,IY,IZ,IEL) = RNOR*VNY(IX,IY,IZ,IEL) +
2123 $ RTN1*V1Y(IX,IY,IZ,IEL) +
2124 $ RTN2*V2Y(IX,IY,IZ,IEL)
2125 R3(IX,IY,IZ,IEL) = RNOR*VNZ(IX,IY,IZ,IEL) +
2126 $ RTN1*V1Z(IX,IY,IZ,IEL) +
2127 $ RTN2*V2Z(IX,IY,IZ,IEL)
2128 100 CONTINUE
2129C
2130 RETURN
2131 END
2132C-----------------------------------------------------------------------
2133 subroutine getSnormal(sn,ix,iy,iz,iside,e)
2134
2135c calculate surface normal
2136
2137 include 'SIZE'
2138 include 'GEOM'
2139 include 'TOPOL'
2140
2141 real sn(3)
2142 integer e,f
2143
2144 f = eface1(iside)
2145
2146 if (1.le.f.and.f.le.2) then ! "r face"
2147 sn(1) = unx(iy,iz,iside,e)
2148 sn(2) = uny(iy,iz,iside,e)
2149 sn(3) = unz(iy,iz,iside,e)
2150 elseif (3.le.f.and.f.le.4) then ! "s face"
2151 sn(1) = unx(ix,iz,iside,e)
2152 sn(2) = uny(ix,iz,iside,e)
2153 sn(3) = unz(ix,iz,iside,e)
2154 elseif (5.le.f.and.f.le.6) then ! "t face"
2155 sn(1) = unx(ix,iy,iside,e)
2156 sn(2) = uny(ix,iy,iside,e)
2157 sn(3) = unz(ix,iy,iside,e)
2158 endif
2159
2160 return
2161 end
2162
2163 subroutine fixmska (c1mask,c2mask,c3mask)
2164
2165c fixes masks for A/SYM face corners
2166
2167 include 'SIZE'
2168 include 'INPUT'
2169
2170 real c1mask(lx1,ly1,lz1,1)
2171 $ ,c2mask(lx1,ly1,lz1,1)
2172 $ ,c3mask(lx1,ly1,lz1,1)
2173
2174 common /ctmp0/ im1(lx1,ly1,lz1),im2(lx1,ly1,lz1)
2175 integer e,f,val,im1,im2
2176
2177 character*3 cb
2178
2179 n = lx1*ly1*lz1
2180
2181 nface = 2*ldim
2182
2183 do e=1,nelv
2184 call izero (im1,n)
2185 call izero (im2,n)
2186 do f=1,nface
2187 cb = cbc (f,e,1)
2188 if (cb.eq.'SYM') call iface_e(im1,f,1,lx1,ly1,lz1)
2189 if (cb.eq.'A ') call iface_e(im2,f,2,lx1,ly1,lz1)
2190 enddo
2191 call icol2(im2,im1,n)
2192
2193 k = 1
2194 do j=1,ly1,ly1-1
2195 do i=1,lx1,lx1-1
2196 if ( im2(i,j,k) .eq. 2) then ! corner of SYM & 'A ' faces
2197 c1mask(i,j,k,e) = 0.
2198 c2mask(i,j,k,e) = 0.
2199 endif
2200 enddo
2201 enddo
2202 enddo
2203
2204 return
2205 end
2206c-----------------------------------------------------------------------
2207 subroutine icol2(a,b,n)
2208 integer a(1),b(1)
2209
2210 do i=1,n
2211 a(i)=a(i)*b(i)
2212 enddo
2213
2214 return
2215 end
2216c-----------------------------------------------------------------------
2217 subroutine iface_e(a,iface,val,nx,ny,nz)
2218
2219C Assign the value VAL to face(IFACE,IE) of array A.
2220C IFACE is the input in the pre-processor ordering scheme.
2221
2222 include 'SIZE'
2223 integer a(nx,ny,nz),val
2224 call facind (kx1,kx2,ky1,ky2,kz1,kz2,nx,ny,nz,iface)
2225 do 100 iz=kz1,kz2
2226 do 100 iy=ky1,ky2
2227 do 100 ix=kx1,kx2
2228 a(ix,iy,iz)=val
2229 100 continue
2230 return
2231 end
2232c-----------------------------------------------------------------------
2233 function op_vlsc2_wt(b1,b2,b3,x1,x2,x3,wt)
2234 include 'SIZE'
2235 include 'INPUT'
2236 include 'TSTEP'
2237 real b1(1),b2(1),b3(1),x1(1),x2(1),x3(1),wt(1)
2238
2239 nel = nelfld(ifield)
2240 n = lx1*ly1*lz1*nel
2241
2242 s = 0
2243 if (if3d) then
2244 do i=1,n
2245 s=s+wt(i)*(b1(i)*x1(i)+b2(i)*x2(i)+b3(i)*x3(i))
2246 enddo
2247 else
2248 do i=1,n
2249 s=s+wt(i)*(b1(i)*x1(i)+b2(i)*x2(i))
2250 enddo
2251 endif
2252 op_vlsc2_wt = s
2253
2254 return
2255 end
2256c-----------------------------------------------------------------------
2257 function op_glsc2_wt(b1,b2,b3,x1,x2,x3,wt)
2258 include 'SIZE'
2259 include 'INPUT'
2260 include 'TSTEP'
2261 real b1(1),b2(1),b3(1),x1(1),x2(1),x3(1),wt(1)
2262
2263 nel = nelfld(ifield)
2264 n = lx1*ly1*lz1*nel
2265
2266 s = 0
2267 if (if3d) then
2268 do i=1,n
2269 s=s+wt(i)*(b1(i)*x1(i)+b2(i)*x2(i)+b3(i)*x3(i))
2270 enddo
2271 else
2272 do i=1,n
2273 s=s+wt(i)*(b1(i)*x1(i)+b2(i)*x2(i))
2274 enddo
2275 endif
2276 op_glsc2_wt = glsum(s,1)
2277
2278 return
2279 end
2280c-----------------------------------------------------------------------
Note: See TracBrowser for help on using the repository browser.