CG  Version 25
consCoeff.h
Go to the documentation of this file.
1 ! ========== This include file contains conservative approximations to coefficient operators ================
2 !
3 ! ajac2d(i1,i2,i3)
4 ! ajac3d(i1,i2,i3
5 ! getCoeffForDxADxPlusDyBDy(au, azmz,amzz,azzz,apzz,azpz, bzmz,bmzz,bzzz,bpzz,bzpz )
6 ! getCoeffForDyADx( au, azmz,amzz,azzz,apzz,azpz )
7 ! getCoeffForDxADy( au, azmz,amzz,azzz,apzz,azpz )
8 ! setDivTensorGradCoeff2d(cmp,eqn,a11ph,a11mh,a22ph,a22mh,a12pzz,a12mzz,a21zpz,a21zmz)
9 ! scaleCoefficients( a11ph,a11mh,a22ph,a22mh,a12pzz,a12mzz,a21zpz,a21zmz )
10 ! getCoeffForDxADxPlusDyBDyPlusDzCDz(au, azzm,azmz,amzz,azzz,apzz,azpz,azzp,...)
11 ! getCoeffForDxADy3d(au, X, Y, azzm,azmz,amzz,azzz,apzz,azpz,azzp )
12 
13 #defineMacro ajac2d(i1,i2,i3) (1./(rx(i1,i2,i3)*sy(i1,i2,i3)-ry(i1,i2,i3)*sx(i1,i2,i3)))
14 #defineMacro ajac3d(i1,i2,i3) (1./((rx(i1,i2,i3)*sy(i1,i2,i3)-ry(i1,i2,i3)*sx(i1,i2,i3))*tz(i1,i2,i3)+\
15  (ry(i1,i2,i3)*sz(i1,i2,i3)-rz(i1,i2,i3)*sy(i1,i2,i3))*tx(i1,i2,i3)+\
16  (rz(i1,i2,i3)*sx(i1,i2,i3)-rx(i1,i2,i3)*sz(i1,i2,i3))*ty(i1,i2,i3)))
17 
18 ! =============================================================================================================
19 ! Declare variables used in compute the coefficients for a non-linear viscosity
20 ! =============================================================================================================
21 #beginMacro declareNonLinearViscosityVariables()
22 
23  real nu0ph,nu0mh,nu1ph,nu1mh,nu2ph,nu2mh
24  real nuzzm,nuzmz,numzz,nuzzz,nupzz,nuzpz,nuzzp
25  real ajzzm,ajzmz,ajmzz,ajzzz,ajpzz,ajzpz,ajzzp
26 
27  real a11ph,a11mh,a22ph,a22mh,a33ph,a33mh,a11mzz,a11zzz,a11pzz,a22zmz,a22zzz,a22zpz,a33zzm,a33zzz,a33zzp,\
28  a12pzz,a12zzz,a12mzz,a13pzz,a13zzz,a13mzz,a21zpz,a21zzz,a21zmz,a23zpz,a23zzz,a23zmz,\
29  a31zzp,a31zzz,a31zzm,a32zzp,a32zzz,a32zzm
30  real b11ph,b11mh,b22ph,b22mh,b33ph,b33mh,b11mzz,b11zzz,b11pzz,b22zmz,b22zzz,b22zpz,b33zzm,b33zzz,b33zzp,\
31  b12pzz,b12zzz,b12mzz,b13pzz,b13zzz,b13mzz,b21zpz,b21zzz,b21zmz,b23zpz,b23zzz,b23zmz,\
32  b31zzp,b31zzz,b31zzm,b32zzp,b32zzz,b32zzm
33  real c11ph,c11mh,c22ph,c22mh,c33ph,c33mh,c11mzz,c11zzz,c11pzz,c22zmz,c22zzz,c22zpz,c33zzm,c33zzz,c33zzp,\
34  c12pzz,c12zzz,c12mzz,c13pzz,c13zzz,c13mzz,c21zpz,c21zzz,c21zmz,c23zpz,c23zzz,c23zmz,\
35  c31zzp,c31zzz,c31zzm,c32zzp,c32zzz,c32zzm
36 
37  real au11ph,au11mh,au22ph,au22mh,au33ph,au33mh,au11mzz,au11zzz,au11pzz,au22zmz,au22zzz,au22zpz,au33zzm,au33zzz,au33zzp,\
38  au12pzz,au12zzz,au12mzz,au13pzz,au13zzz,au13mzz,au21zpz,au21zzz,au21zmz,au23zpz,au23zzz,au23zmz,\
39  au31zzp,au31zzz,au31zzm,au32zzp,au32zzz,au32zzm
40  real av11ph,av11mh,av22ph,av22mh,av33ph,av33mh,av11mzz,av11zzz,av11pzz,av22zmz,av22zzz,av22zpz,av33zzm,av33zzz,av33zzp,\
41  av12pzz,av12zzz,av12mzz,av13pzz,av13zzz,av13mzz,av21zpz,av21zzz,av21zmz,av23zpz,av23zzz,av23zmz,\
42  av31zzp,av31zzz,av31zzm,av32zzp,av32zzz,av32zzm
43  real aw11ph,aw11mh,aw22ph,aw22mh,aw33ph,aw33mh,aw11mzz,aw11zzz,aw11pzz,aw22zmz,aw22zzz,aw22zpz,aw33zzm,aw33zzz,aw33zzp,\
44  aw12pzz,aw12zzz,aw12mzz,aw13pzz,aw13zzz,aw13mzz,aw21zpz,aw21zzz,aw21zmz,aw23zpz,aw23zzz,aw23zmz,\
45  aw31zzp,aw31zzz,aw31zzm,aw32zzp,aw32zzz,aw32zzm
46 
47  real bu11ph,bu11mh,bu22ph,bu22mh,bu33ph,bu33mh,bu11mzz,bu11zzz,bu11pzz,bu22zmz,bu22zzz,bu22zpz,bu33zzm,bu33zzz,bu33zzp,\
48  bu12pzz,bu12zzz,bu12mzz,bu13pzz,bu13zzz,bu13mzz,bu21zpz,bu21zzz,bu21zmz,bu23zpz,bu23zzz,bu23zmz,\
49  bu31zzp,bu31zzz,bu31zzm,bu32zzp,bu32zzz,bu32zzm
50  real bv11ph,bv11mh,bv22ph,bv22mh,bv33ph,bv33mh,bv11mzz,bv11zzz,bv11pzz,bv22zmz,bv22zzz,bv22zpz,bv33zzm,bv33zzz,bv33zzp,\
51  bv12pzz,bv12zzz,bv12mzz,bv13pzz,bv13zzz,bv13mzz,bv21zpz,bv21zzz,bv21zmz,bv23zpz,bv23zzz,bv23zmz,\
52  bv31zzp,bv31zzz,bv31zzm,bv32zzp,bv32zzz,bv32zzm
53  real bw11ph,bw11mh,bw22ph,bw22mh,bw33ph,bw33mh,bw11mzz,bw11zzz,bw11pzz,bw22zmz,bw22zzz,bw22zpz,bw33zzm,bw33zzz,bw33zzp,\
54  bw12pzz,bw12zzz,bw12mzz,bw13pzz,bw13zzz,bw13mzz,bw21zpz,bw21zzz,bw21zmz,bw23zpz,bw23zzz,bw23zmz,\
55  bw31zzp,bw31zzz,bw31zzm,bw32zzp,bw32zzz,bw32zzm
56 
57  real cu11ph,cu11mh,cu22ph,cu22mh,cu33ph,cu33mh,cu11mzz,cu11zzz,cu11pzz,cu22zmz,cu22zzz,cu22zpz,cu33zzm,cu33zzz,cu33zzp,\
58  cu12pzz,cu12zzz,cu12mzz,cu13pzz,cu13zzz,cu13mzz,cu21zpz,cu21zzz,cu21zmz,cu23zpz,cu23zzz,cu23zmz,\
59  cu31zzp,cu31zzz,cu31zzm,cu32zzp,cu32zzz,cu32zzm
60  real cv11ph,cv11mh,cv22ph,cv22mh,cv33ph,cv33mh,cv11mzz,cv11zzz,cv11pzz,cv22zmz,cv22zzz,cv22zpz,cv33zzm,cv33zzz,cv33zzp,\
61  cv12pzz,cv12zzz,cv12mzz,cv13pzz,cv13zzz,cv13mzz,cv21zpz,cv21zzz,cv21zmz,cv23zpz,cv23zzz,cv23zmz,\
62  cv31zzp,cv31zzz,cv31zzm,cv32zzp,cv32zzz,cv32zzm
63  real cw11ph,cw11mh,cw22ph,cw22mh,cw33ph,cw33mh,cw11mzz,cw11zzz,cw11pzz,cw22zmz,cw22zzz,cw22zpz,cw33zzm,cw33zzz,cw33zzp,\
64  cw12pzz,cw12zzz,cw12mzz,cw13pzz,cw13zzz,cw13mzz,cw21zpz,cw21zzz,cw21zmz,cw23zpz,cw23zzz,cw23zmz,\
65  cw31zzp,cw31zzz,cw31zzm,cw32zzp,cw32zzz,cw32zzm
66 
67 
68 #endMacro
69 
70 
71 
72 ! ==================================================================================================
73 ! Define the coefficients in the conservative discretization of:
74 ! L = Dx( a*Dx ) + Dy( b*Dy )
75 !
76 ! L = (1/J)*[ Dr( J*(rx,ry).(aDx,bDy)) + Ds( J*(sx,sy).(aDx,bDy)) ]
77 ! = (1/J)*[ Dr( J*a*rx(rx*Dr + sx*Ds) + J*b*ry*(ry*Dr + sy*Ds)
78 ! +Ds( J*a*sx(rx*Dr + sx*Ds) + J*b*sy*(ry*Dr + sy*Ds) ]
79 ! = (1/J)*[ Dr( a11 Dr) + Dr( a12 Ds) + Ds( a21 Dr) + Ds( a22 Ds) ]
80 ! where
81 ! a11 = J ( a rx^2 + b ry^2 )
82 ! a12 = J ( a rx*sx + b ry*sy )
83 ! a21 = a12
84 ! a22 = J ( a sx^2 + b sy^2 )
85 ! a = a(i1,i2,i3), b=b(i1,i2,i3)
86 !
87 ! Macro Arguments:
88 ! au : prefix of the computed coefficients
89 ! azmz,amzz,azzz,apzz,azpz : a(i1,i2-1,i2),a(i1-1,i2,i3),a(i1,i2,i3),a(i1+1,i2,i3),a(i1,i2+1,i3)
90 ! bzmz,bmzz,bzzz,bpzz,bzpz : b(i1,i2-1,i2),b(i1-1,i2,i3),b(i1,i2,i3),b(i1+1,i2,i3),b(i1,i2+1,i3)
91 ! The following jacobian values should also be defined:
92 ! ajzmz,ajmzz,ajzzz,ajpzz,ajzpz : aj(i1,i2-1,i2),aj(i1-1,i2,i3),aj(i1,i2,i3),aj(i1+1,i2,i3),aj(i1,i2+1,i3)
93 ! ==================================================================================================
94 #beginMacro getCoeffForDxADxPlusDyBDy(au, azmz,amzz,azzz,apzz,azpz, bzmz,bmzz,bzzz,bpzz,bzpz )
95  au ## 11mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*rx(i1-1,i2,i3) + (bmzz)*ry(i1-1,i2,i3)*ry(i1-1,i2,i3) )
96  au ## 11zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*rx(i1 ,i2,i3) + (bzzz)*ry(i1 ,i2,i3)*ry(i1 ,i2,i3) )
97  au ## 11pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*rx(i1+1,i2,i3) + (bpzz)*ry(i1+1,i2,i3)*ry(i1+1,i2,i3) )
98  au ## 11ph = .5*( au ## 11zzz+au ## 11pzz )
99  au ## 11mh = .5*( au ## 11zzz+au ## 11mzz )
100 
101  au ## 22zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*sx(i1,i2-1,i3) + (bzmz)*sy(i1,i2-1,i3)*sy(i1,i2-1,i3) )
102  au ## 22zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*sx(i1,i2 ,i3) + (bzzz)*sy(i1,i2 ,i3)*sy(i1,i2 ,i3) )
103  au ## 22zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*sx(i1,i2+1,i3) + (bzpz)*sy(i1,i2+1,i3)*sy(i1,i2+1,i3) )
104  au ## 22ph = .5*( au ## 22zzz+au ## 22zpz )
105  au ## 22mh = .5*( au ## 22zzz+au ## 22zmz )
106 
107  au ## 12mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*sx(i1-1,i2,i3) + (bmzz)*ry(i1-1,i2,i3)*sy(i1-1,i2,i3) )
108  au ## 12zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*sx(i1 ,i2,i3) + (bzzz)*ry(i1 ,i2,i3)*sy(i1 ,i2,i3) )
109  au ## 12pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*sx(i1+1,i2,i3) + (bpzz)*ry(i1+1,i2,i3)*sy(i1+1,i2,i3) )
110 
111  au ## 21zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*rx(i1,i2-1,i3) + (bzmz)*sy(i1,i2-1,i3)*ry(i1,i2-1,i3) )
112  au ## 21zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*rx(i1,i2 ,i3) + (bzzz)*sy(i1,i2 ,i3)*ry(i1,i2 ,i3) )
113  au ## 21zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*rx(i1,i2+1,i3) + (bzpz)*sy(i1,i2+1,i3)*ry(i1,i2+1,i3) )
114 #endMacro
115 
116 ! ==================================================================================================
117 ! Define the coefficients in the conservative discretization of:
118 ! L = Dy( a*Dx )
119 ! L = div( (0,aDx) )
120 ! L = (1/J)*[ Dr( J*(rx,ry).(0,aDx)) + Ds( J*(sx,sy).(0,aDx)) ]
121 ! = (1/J)*[ Dr( J*a*ry*(rx*Dr + sx*Ds)) + Ds( J*a*sy*(rx*Dr + sx*Ds) ) ]
122 ! = (1/J)*( Dr( a11 Dr) + Dr( a12 Ds) + Ds( a21 Dr) + Ds( a22 Ds) ) u
123 ! where
124 ! a11 = J ( a ry*rx )
125 ! a12 = J ( a ry*sx )
126 ! a21 = J ( a sy*rx )
127 ! a22 = J ( a sy*sx )
128 ! a = a(i1,i2,i3)
129 !
130 ! Macro Arguments:
131 ! au : prefix of the computed coefficients
132 ! azmz,amzz,azzz,apzz,azpz : a(i1,i2-1,i2),a(i1-1,i2,i3),a(i1,i2,i3),a(i1+1,i2,i3),a(i1,i2+1,i3)
133 ! The following jacobian values should also be defined:
134 ! ajzmz,ajmzz,ajzzz,ajpzz,ajzpz : aj(i1,i2-1,i2),aj(i1-1,i2,i3),aj(i1,i2,i3),aj(i1+1,i2,i3),aj(i1,i2+1,i3)
135 ! ==================================================================================================
136 #beginMacro getCoeffForDyADx( au, azmz,amzz,azzz,apzz,azpz )
137 
138  au ## 11mzz = ajmzz*( (amzz)*ry(i1-1,i2,i3)*rx(i1-1,i2,i3) )
139  au ## 11zzz = ajzzz*( (azzz)*ry(i1 ,i2,i3)*rx(i1 ,i2,i3) )
140  au ## 11pzz = ajpzz*( (apzz)*ry(i1+1,i2,i3)*rx(i1+1,i2,i3) )
141  au ## 11ph = .5*( au ## 11zzz+au ## 11pzz )
142  au ## 11mh = .5*( au ## 11zzz+au ## 11mzz )
143 
144  au ## 22zmz = ajzmz*( (azmz)*sy(i1,i2-1,i3)*sx(i1,i2-1,i3) )
145  au ## 22zzz = ajzzz*( (azzz)*sy(i1,i2 ,i3)*sx(i1,i2 ,i3) )
146  au ## 22zpz = ajzpz*( (azpz)*sy(i1,i2+1,i3)*sx(i1,i2+1,i3) )
147  au ## 22ph = .5*( au ## 22zzz+au ## 22zpz )
148  au ## 22mh = .5*( au ## 22zzz+au ## 22zmz )
149 
150  au ## 12mzz = ajmzz*( (amzz)*ry(i1-1,i2,i3)*sx(i1-1,i2,i3) )
151  au ## 12zzz = ajzzz*( (azzz)*ry(i1 ,i2,i3)*sx(i1 ,i2,i3) )
152  au ## 12pzz = ajpzz*( (apzz)*ry(i1+1,i2,i3)*sx(i1+1,i2,i3) )
153 
154  au ## 21zmz = ajzmz*( (azmz)*sy(i1,i2-1,i3)*rx(i1,i2-1,i3) )
155  au ## 21zzz = ajzzz*( (azzz)*sy(i1,i2 ,i3)*rx(i1,i2 ,i3) )
156  au ## 21zpz = ajzpz*( (azpz)*sy(i1,i2+1,i3)*rx(i1,i2+1,i3) )
157 
158 #endMacro
159 
160 ! ==================================================================================================
161 ! Define the coefficients in the conservative discretization of:
162 ! L = Dx( a*Dy )
163 !
164 ! L = div( (aDy,0) )
165 ! L = (1/J)*[ Dr( J*(rx,ry).(aDy,0)) + Ds(J*(sx,sy).(aDy,0)) ]
166 ! = (1/J)*[ Dr( J*a*rx*(ry*Dr + sy*Ds)) + Ds( J*a*sx*(ry*Dr + sy*Ds) ) ]
167 ! = (1/J)*[ Dr( a11 Dr) + Dr( a12 Ds) + Ds( a21 Dr) + Ds( a22 Ds) ]
168 ! where
169 ! a11 = J ( a rx*ry )
170 ! a12 = J ( a rx*sy )
171 ! a21 = J ( a sx*ry )
172 ! a22 = J ( a sx*sy )
173 ! Macro Arguments:
174 ! au : prefix of the computed coefficients
175 ! azmz,amzz,azzz,apzz,azpz : a(i1,i2-1,i2),a(i1-1,i2,i3),a(i1,i2,i3),a(i1+1,i2,i3),a(i1,i2+1,i3)
176 ! The following jacobian values should also be defined:
177 ! ajzmz,ajmzz,ajzzz,ajpzz,ajzpz : aj(i1,i2-1,i2),aj(i1-1,i2,i3),aj(i1,i2,i3),aj(i1+1,i2,i3),aj(i1,i2+1,i3)
178 ! ==================================================================================================
179 #beginMacro getCoeffForDxADy( au, azmz,amzz,azzz,apzz,azpz )
180 
181  au ## 11mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*ry(i1-1,i2,i3) )
182  au ## 11zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*ry(i1 ,i2,i3) )
183  au ## 11pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*ry(i1+1,i2,i3) )
184  au ## 11ph = .5*( au ## 11zzz+au ## 11pzz )
185  au ## 11mh = .5*( au ## 11zzz+au ## 11mzz )
186 
187  au ## 22zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*sy(i1,i2-1,i3) )
188  au ## 22zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*sy(i1,i2 ,i3) )
189  au ## 22zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*sy(i1,i2+1,i3) )
190  au ## 22ph = .5*( au ## 22zzz+au ## 22zpz )
191  au ## 22mh = .5*( au ## 22zzz+au ## 22zmz )
192 
193  au ## 12mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*sy(i1-1,i2,i3) )
194  au ## 12zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*sy(i1 ,i2,i3) )
195  au ## 12pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*sy(i1+1,i2,i3) )
196 
197  au ## 21zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*ry(i1,i2-1,i3) )
198  au ## 21zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*ry(i1,i2 ,i3) )
199  au ## 21zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*ry(i1,i2+1,i3) )
200 
201 #endMacro
202 
203 
204 ! =============================================================================================================
205 ! Assign the coefficients for a component of the conservative discretization of the div(tensor grad) operator
206 ! =============================================================================================================
207 #beginMacro setDivTensorGradCoeff2d(cmp,eqn,a11ph,a11mh,a22ph,a22mh,a12pzz,a12mzz,a21zpz,a21zmz)
208  coeff(MCE(-1,-1,0,cmp,eqn),i1,i2,i3)= a12mzz+a21zmz
209  coeff(MCE( 0,-1,0,cmp,eqn),i1,i2,i3)= a22mh
210  coeff(MCE( 1,-1,0,cmp,eqn),i1,i2,i3)= -a12pzz-a21zmz
211  coeff(MCE(-1, 0,0,cmp,eqn),i1,i2,i3)= a11mh
212  coeff(MCE( 0, 0,0,cmp,eqn),i1,i2,i3)= -a11ph-a11mh -a22ph -a22mh
213  coeff(MCE( 1, 0,0,cmp,eqn),i1,i2,i3)= a11ph
214  coeff(MCE(-1, 1,0,cmp,eqn),i1,i2,i3)= -a12mzz-a21zpz
215  coeff(MCE( 0, 1,0,cmp,eqn),i1,i2,i3)= a22ph
216  coeff(MCE( 1, 1,0,cmp,eqn),i1,i2,i3)= a12pzz+a21zpz
217 #endMacro
218 
219 
220 ! =======================================================================================================
221 ! This macro scaled the coefficients that appear in the discretization of the div(tensor grad) operator
222 ! =======================================================================================================
223 #beginMacro scaleCoefficients( a11ph,a11mh,a22ph,a22mh,a12pzz,a12mzz,a21zpz,a21zmz )
224  a11ph=a11ph*dr0i
225  a11mh=a11mh*dr0i
226  a22ph=a22ph*dr1i
227  a22mh=a22mh*dr1i
228  a12pzz=a12pzz*dr0dr1
229  a12mzz=a12mzz*dr0dr1
230  a21zpz=a21zpz*dr0dr1
231  a21zmz=a21zmz*dr0dr1
232 #endMacro
233 
234 
235 
236 ! ****************************************************************
237 ! ****************** THREE DIMENSIONS ****************************
238 ! ****************************************************************
239 
240 
241 ! ==================================================================================================
242 ! Define the coefficients in the conservative discretization of:
243 ! L = Dx( a*Dx ) + Dy( b*Dy ) + Dz( c*Dz )
244 !
245 ! L = (1/J)*[ Dr( J*(rx,ry,rz).(aDx,bDy,cDz)) + Ds( J*(sx,sy,sz).(aDx,bDy,cDz)) + Dt( J*(tx,ty,tz).(aDx,bDy,cDz)) ]
246 ! = (1/J)*[ Dr( J*a*rx(rx*Dr + sx*Ds+ tx*Dt) + J*b*ry*(ry*Dr + sy*Ds+ ty*Dt) + J*c*rz*(rz*Dr + sz*Ds+ tz*Dt)
247 ! +Ds( J*a*sx(rx*Dr + sx*Ds+ tx*Dt) + J*b*sy*(ry*Dr + sy*Ds+ ty*Dt) + J*c*sz*(rz*Dr + sz*Ds+ tz*Dt)
248 ! +Dt( J*a*tx(rx*Dr + sx*Ds+ tx*Dt) + J*b*ty*(ry*Dr + sy*Ds+ ty*Dt) + J*c*tz*(rz*Dr + sz*Ds+ tz*Dt)
249 ! = (1/J)*[ Dr( a11 Dr) + Dr( a12 Ds) + Dr( a13 Dt) + Ds( a21 Dr) + Ds( a22 Ds) + Ds( a23 Dt) + Dt( a31 Dr) + Dt( a32 Ds) + Dt( a33 Dt) ]
250 ! where
251 ! a11 = J ( a rx^2 + b ry^2 + c rz^2 )
252 ! a12 = J ( a rx*sx + b ry*sy + b rz*sz )
253 ! a13 = J ( a rx*tx + b ry*ty + b rz*tz )
254 ! a21 = a12
255 ! a22 = J ( a sx^2 + b sy^2 + c sz^2 )
256 ! a23 = J ( a sx*tx + b sy*ty + b sz*tz )
257 ! a31 = a13
258 ! a32 = a23
259 ! a33 = J ( a tx^2 + b ty^2 + c tz^2 )
260 ! a = a(i1,i2,i3), b=b(i1,i2,i3)
261 !
262 ! Macro Arguments:
263 ! au : prefix of the computed coefficients
264 ! azzm,azmz,amzz,azzz,apzz,azpz,azzp : a(i1,i2,i3-1),a(i1,i2-1,i3),a(i1-1,i2,i3),a(i1,i2,i3),a(i1+1,i2,i3),a(i1,i2+1,i3),a(i1,i2,i3+1)
265 ! bzzm,bzmz,bmzz,bzzz,bpzz,bzpz,bzzp : b(i1,i2,i3-1),b(i1,i2-1,i3),b(i1-1,i2,i3),b(i1,i2,i3),b(i1+1,i2,i3),b(i1,i2+1,i3),b(i1,i2,i3+1)
266 ! czzm,czmz,cmzz,czzz,cpzz,czpz,czzp : c(i1,i2,i3-1),c(i1,i2-1,i3),c(i1-1,i2,i3),c(i1,i2,i3),c(i1+1,i2,i3),c(i1,i2+1,i3),c(i1,i2,i3+1)
267 ! The following jacobian values should also be defined:
268 ! ajzmz,ajzzm,ajmzz,ajzzz,ajpzz,ajzpz,ajzzp : aj(i1,i2,i3-1),aj(i1,i2-1,i3),aj(i1-1,i2,i3),aj(i1,i2,i3),aj(i1+1,i2,i3),aj(i1,i2+1,i3),aj(i1,i2,i3+1)
269 ! ==================================================================================================
270 #beginMacro getCoeffForDxADxPlusDyBDyPlusDzCDz(au, azzm,azmz,amzz,azzz,apzz,azpz,azzp,\
271  bzzm,bzmz,bmzz,bzzz,bpzz,bzpz,bzzp,\
272  czzm,czmz,cmzz,czzz,cpzz,czpz,czzp )
273  au ## 11mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*rx(i1-1,i2,i3) + (bmzz)*ry(i1-1,i2,i3)*ry(i1-1,i2,i3) + (cmzz)*rz(i1-1,i2,i3)*rz(i1-1,i2,i3) )
274  au ## 11zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*rx(i1 ,i2,i3) + (bzzz)*ry(i1 ,i2,i3)*ry(i1 ,i2,i3) + (czzz)*rz(i1 ,i2,i3)*rz(i1 ,i2,i3) )
275  au ## 11pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*rx(i1+1,i2,i3) + (bpzz)*ry(i1+1,i2,i3)*ry(i1+1,i2,i3) + (cpzz)*rz(i1+1,i2,i3)*rz(i1+1,i2,i3) )
276  au ## 11ph = .5*( au ## 11zzz+au ## 11pzz )
277  au ## 11mh = .5*( au ## 11zzz+au ## 11mzz )
278 
279  au ## 22zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*sx(i1,i2-1,i3) + (bzmz)*sy(i1,i2-1,i3)*sy(i1,i2-1,i3) + (czmz)*sz(i1,i2-1,i3)*sz(i1,i2-1,i3) )
280  au ## 22zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*sx(i1,i2 ,i3) + (bzzz)*sy(i1,i2 ,i3)*sy(i1,i2 ,i3) + (czzz)*sz(i1,i2 ,i3)*sz(i1,i2 ,i3) )
281  au ## 22zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*sx(i1,i2+1,i3) + (bzpz)*sy(i1,i2+1,i3)*sy(i1,i2+1,i3) + (czpz)*sz(i1,i2+1,i3)*sz(i1,i2+1,i3) )
282  au ## 22ph = .5*( au ## 22zzz+au ## 22zpz )
283  au ## 22mh = .5*( au ## 22zzz+au ## 22zmz )
284 
285  au ## 33zzm = ajzzm*( (azzm)*tx(i1,i2,i3-1)*tx(i1,i2,i3-1) + (bzzm)*ty(i1,i2,i3-1)*ty(i1,i2,i3-1) + (czzm)*tz(i1,i2,i3-1)*tz(i1,i2,i3-1) )
286  au ## 33zzz = ajzzz*( (azzz)*tx(i1,i2,i3 )*tx(i1,i2,i3 ) + (bzzz)*ty(i1,i2,i3 )*ty(i1,i2,i3 ) + (czzz)*tz(i1,i2,i3 )*tz(i1,i2,i3 ) )
287  au ## 33zzp = ajzzp*( (azzp)*tx(i1,i2,i3+1)*tx(i1,i2,i3+1) + (bzzp)*ty(i1,i2,i3+1)*ty(i1,i2,i3+1) + (czzp)*tz(i1,i2,i3+1)*tz(i1,i2,i3+1) )
288  au ## 33ph = .5*( au ## 33zzz+au ## 33zzp )
289  au ## 33mh = .5*( au ## 33zzz+au ## 33zzm )
290 
291 
292  au ## 12mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*sx(i1-1,i2,i3) + (bmzz)*ry(i1-1,i2,i3)*sy(i1-1,i2,i3) + (cmzz)*rz(i1-1,i2,i3)*sz(i1-1,i2,i3) )
293  au ## 12zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*sx(i1 ,i2,i3) + (bzzz)*ry(i1 ,i2,i3)*sy(i1 ,i2,i3) + (czzz)*rz(i1 ,i2,i3)*sz(i1 ,i2,i3) )
294  au ## 12pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*sx(i1+1,i2,i3) + (bpzz)*ry(i1+1,i2,i3)*sy(i1+1,i2,i3) + (cpzz)*rz(i1+1,i2,i3)*sz(i1+1,i2,i3) )
295 
296  au ## 13mzz = ajmzz*( (amzz)*rx(i1-1,i2,i3)*tx(i1-1,i2,i3) + (bmzz)*ry(i1-1,i2,i3)*ty(i1-1,i2,i3) + (cmzz)*rz(i1-1,i2,i3)*tz(i1-1,i2,i3) )
297  au ## 13zzz = ajzzz*( (azzz)*rx(i1 ,i2,i3)*tx(i1 ,i2,i3) + (bzzz)*ry(i1 ,i2,i3)*ty(i1 ,i2,i3) + (czzz)*rz(i1 ,i2,i3)*tz(i1 ,i2,i3) )
298  au ## 13pzz = ajpzz*( (apzz)*rx(i1+1,i2,i3)*tx(i1+1,i2,i3) + (bpzz)*ry(i1+1,i2,i3)*ty(i1+1,i2,i3) + (cpzz)*rz(i1+1,i2,i3)*tz(i1+1,i2,i3) )
299 
300  au ## 21zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*rx(i1,i2-1,i3) + (bzmz)*sy(i1,i2-1,i3)*ry(i1,i2-1,i3) + (czmz)*sz(i1,i2-1,i3)*rz(i1,i2-1,i3) )
301  au ## 21zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*rx(i1,i2 ,i3) + (bzzz)*sy(i1,i2 ,i3)*ry(i1,i2 ,i3) + (czzz)*sz(i1,i2 ,i3)*rz(i1,i2 ,i3) )
302  au ## 21zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*rx(i1,i2+1,i3) + (bzpz)*sy(i1,i2+1,i3)*ry(i1,i2+1,i3) + (czpz)*sz(i1,i2+1,i3)*rz(i1,i2+1,i3) )
303 
304  au ## 23zmz = ajzmz*( (azmz)*sx(i1,i2-1,i3)*tx(i1,i2-1,i3) + (bzmz)*sy(i1,i2-1,i3)*ty(i1,i2-1,i3) + (czmz)*sz(i1,i2-1,i3)*tz(i1,i2-1,i3) )
305  au ## 23zzz = ajzzz*( (azzz)*sx(i1,i2 ,i3)*tx(i1,i2 ,i3) + (bzzz)*sy(i1,i2 ,i3)*ty(i1,i2 ,i3) + (czzz)*sz(i1,i2 ,i3)*tz(i1,i2 ,i3) )
306  au ## 23zpz = ajzpz*( (azpz)*sx(i1,i2+1,i3)*tx(i1,i2+1,i3) + (bzpz)*sy(i1,i2+1,i3)*ty(i1,i2+1,i3) + (czpz)*sz(i1,i2+1,i3)*tz(i1,i2+1,i3) )
307 
308 
309  au ## 31zzm = ajzzm*( (azzm)*tx(i1,i2,i3-1)*rx(i1,i2,i3-1) + (bzzm)*ty(i1,i2,i3-1)*ry(i1,i2,i3-1) + (czzm)*tz(i1,i2,i3-1)*rz(i1,i2,i3-1) )
310  au ## 31zzz = ajzzz*( (azzz)*tx(i1,i2,i3 )*rx(i1,i2,i3 ) + (bzzz)*ty(i1,i2,i3 )*ry(i1,i2,i3 ) + (czzz)*tz(i1,i2,i3 )*rz(i1,i2,i3 ) )
311  au ## 31zzp = ajzzp*( (azzp)*tx(i1,i2,i3+1)*rx(i1,i2,i3+1) + (bzzp)*ty(i1,i2,i3+1)*ry(i1,i2,i3+1) + (czzp)*tz(i1,i2,i3+1)*rz(i1,i2,i3+1) )
312 
313  au ## 32zzm = ajzzm*( (azzm)*tx(i1,i2,i3-1)*sx(i1,i2,i3-1) + (bzzm)*ty(i1,i2,i3-1)*sy(i1,i2,i3-1) + (czzm)*tz(i1,i2,i3-1)*sz(i1,i2,i3-1) )
314  au ## 32zzz = ajzzz*( (azzz)*tx(i1,i2,i3 )*sx(i1,i2,i3 ) + (bzzz)*ty(i1,i2,i3 )*sy(i1,i2,i3 ) + (czzz)*tz(i1,i2,i3 )*sz(i1,i2,i3 ) )
315  au ## 32zzp = ajzzp*( (azzp)*tx(i1,i2,i3+1)*sx(i1,i2,i3+1) + (bzzp)*ty(i1,i2,i3+1)*sy(i1,i2,i3+1) + (czzp)*tz(i1,i2,i3+1)*sz(i1,i2,i3+1) )
316 
317 #endMacro
318 
319 ! ==========================================================================================================================================
320 !
321 ! Define the coefficients in the conservative discretization of any mixed or non-mixed derivative:
322 !
323 ! L = D_X( a*D_Y ) where X=x, y, or z and Y=x, y, or z
324 !
325 ! Example:
326 ! L = div( (aDy,0,0) )
327 ! L = (1/J)*[ Dr( J*(rx,ry,rx).(aDy,0,0)) + Ds(J*(sx,sy,sz).(aDy,0,0))+ Dt(J*(tx,ty,tz).(aDy,0,0)) ]
328 ! = (1/J)*[ Dr( J*a*rx(ry*Dr + sy*Ds+ ty*Dt)
329 ! +Ds( J*a*sx(ry*Dr + sy*Ds+ ty*Dt)
330 ! +Dt( J*a*tx(ry*Dr + sy*Ds+ ty*Dt) ]
331 ! = (1/J)*[ Dr( a11 Dr) + Dr( a12 Ds) + Dr( a13 Dt) + Ds( a21 Dr) + Ds( a22 Ds) + Ds( a23 Dt) + Dt( a31 Dr) + Dt( a32 Ds) + Dt( a33 Dt) ]
332 ! where
333 ! a11 = J ( a rx*ry )
334 ! a12 = J ( a rx*sy )
335 ! a13 = J ( a rx*ty )
336 ! a21 = J ( a sx*ry )
337 ! a22 = J ( a sx*sy )
338 ! a23 = J ( a sx*ty )
339 ! a31 = J ( a tx*ry )
340 ! a32 = J ( a tx*sy )
341 ! a33 = J ( a tx*ty )
342 ! a = a(i1,i2,i3)
343 !
344 ! Macro Arguments:
345 ! au : prefix of the computed coefficients
346 ! X,Y : X=[x,y,z] and Y=[x,y,z] to compute the coeffcients of D_X( a D_Y )
347 ! azzm,azmz,amzz,azzz,apzz,azpz,azzp : a(i1,i2,i3-1),a(i1,i2-1,i3),a(i1-1,i2,i3),a(i1,i2,i3),a(i1+1,i2,i3),a(i1,i2+1,i3),a(i1,i2,i3+1)
348 ! The following jacobian values should also be defined:
349 ! ajzmz,ajzzm,ajmzz,ajzzz,ajpzz,ajzpz,ajzzp : aj(i1,i2,i3-1),aj(i1,i2-1,i3),aj(i1-1,i2,i3),aj(i1,i2,i3),aj(i1+1,i2,i3),aj(i1,i2+1,i3),aj(i1,i2,i3+1)
350 ! ===========================================================================================================================================
351 #beginMacro getCoeffForDxADy3d(au, X, Y, azzm,azmz,amzz,azzz,apzz,azpz,azzp )
352 
353  au ## 11mzz = ajmzz*( (amzz)*r ## X(i1-1,i2,i3)*r ## Y(i1-1,i2,i3) )
354  au ## 11zzz = ajzzz*( (azzz)*r ## X(i1 ,i2,i3)*r ## Y(i1 ,i2,i3) )
355  au ## 11pzz = ajpzz*( (apzz)*r ## X(i1+1,i2,i3)*r ## Y(i1+1,i2,i3) )
356  au ## 11ph = .5*( au ## 11zzz+au ## 11pzz )
357  au ## 11mh = .5*( au ## 11zzz+au ## 11mzz )
358 
359  au ## 22zmz = ajzmz*( (azmz)*s ## X(i1,i2-1,i3)*s ## Y(i1,i2-1,i3) )
360  au ## 22zzz = ajzzz*( (azzz)*s ## X(i1,i2 ,i3)*s ## Y(i1,i2 ,i3) )
361  au ## 22zpz = ajzpz*( (azpz)*s ## X(i1,i2+1,i3)*s ## Y(i1,i2+1,i3) )
362  au ## 22ph = .5*( au ## 22zzz+au ## 22zpz )
363  au ## 22mh = .5*( au ## 22zzz+au ## 22zmz )
364 
365  au ## 33zzm = ajzzm*( (azzm)*t ## X(i1,i2,i3-1)*t ## Y(i1,i2,i3-1) )
366  au ## 33zzz = ajzzz*( (azzz)*t ## X(i1,i2,i3 )*t ## Y(i1,i2,i3 ) )
367  au ## 33zzp = ajzzp*( (azzp)*t ## X(i1,i2,i3+1)*t ## Y(i1,i2,i3+1) )
368  au ## 33ph = .5*( au ## 33zzz+au ## 33zzp )
369  au ## 33mh = .5*( au ## 33zzz+au ## 33zzm )
370 
371 
372  au ## 12mzz = ajmzz*( (amzz)*r ## X(i1-1,i2,i3)*s ## Y(i1-1,i2,i3) )
373  au ## 12zzz = ajzzz*( (azzz)*r ## X(i1 ,i2,i3)*s ## Y(i1 ,i2,i3) )
374  au ## 12pzz = ajpzz*( (apzz)*r ## X(i1+1,i2,i3)*s ## Y(i1+1,i2,i3) )
375 
376  au ## 13mzz = ajmzz*( (amzz)*r ## X(i1-1,i2,i3)*t ## Y(i1-1,i2,i3) )
377  au ## 13zzz = ajzzz*( (azzz)*r ## X(i1 ,i2,i3)*t ## Y(i1 ,i2,i3) )
378  au ## 13pzz = ajpzz*( (apzz)*r ## X(i1+1,i2,i3)*t ## Y(i1+1,i2,i3) )
379 
380  au ## 21zmz = ajzmz*( (azmz)*s ## X(i1,i2-1,i3)*r ## Y(i1,i2-1,i3) )
381  au ## 21zzz = ajzzz*( (azzz)*s ## X(i1,i2 ,i3)*r ## Y(i1,i2 ,i3) )
382  au ## 21zpz = ajzpz*( (azpz)*s ## X(i1,i2+1,i3)*r ## Y(i1,i2+1,i3) )
383 
384  au ## 23zmz = ajzmz*( (azmz)*s ## X(i1,i2-1,i3)*t ## Y(i1,i2-1,i3) )
385  au ## 23zzz = ajzzz*( (azzz)*s ## X(i1,i2 ,i3)*t ## Y(i1,i2 ,i3) )
386  au ## 23zpz = ajzpz*( (azpz)*s ## X(i1,i2+1,i3)*t ## Y(i1,i2+1,i3) )
387 
388 
389  au ## 31zzm = ajzzm*( (azzm)*t ## X(i1,i2,i3-1)*r ## Y(i1,i2,i3-1) )
390  au ## 31zzz = ajzzz*( (azzz)*t ## X(i1,i2,i3 )*r ## Y(i1,i2,i3 ) )
391  au ## 31zzp = ajzzp*( (azzp)*t ## X(i1,i2,i3+1)*r ## Y(i1,i2,i3+1) )
392 
393  au ## 32zzm = ajzzm*( (azzm)*t ## X(i1,i2,i3-1)*s ## Y(i1,i2,i3-1) )
394  au ## 32zzz = ajzzz*( (azzz)*t ## X(i1,i2,i3 )*s ## Y(i1,i2,i3 ) )
395  au ## 32zzp = ajzzp*( (azzp)*t ## X(i1,i2,i3+1)*s ## Y(i1,i2,i3+1) )
396 
397 #endMacro
398 
399 
400 
401 ! ================================================================================================================================================
402 ! Assign the coefficients for a component of the conservative 3D discretization of the div(tensor grad) operator
403 !
404 ! L = (1/J)*[ Dr( a11 Dr) + Dr( a12 Ds) + Dr( a13 Dt) + Ds( a21 Dr) + Ds( a22 Ds) + Ds( a23 Dt) + Dt( a31 Dr) + Dt( a32 Ds) + Dt( a33 Dt) ]
405 ! Dr( a11 Dr)u = a11pzz*( u(i1+1)-u(i1) ) - a11mzz*( u(i1)-u(i1-1))
406 ! = a11pzz*u(i1+1) -(a11pzz+a11mzz)*u(i1) +a11mzz*u(i1-1)
407 ! Dr( a13 Dt ) = D0r( a13 D0t ) = a13pzz*( u(i1+1,i2,i3+1)-u(i1+1,i2,i3-1)) - a13mzz*( u(i1-1,i2,i3+1)-u(i1-1,i2,i3-1) )
408 ! Ds( a23 Dt ) = D0s( a23 D0t ) = a23zpz*( u(i1,i2+1,i3+1)-u(i1,i2+1,i3-1)) - a23zmz*( u(i1,i2-1,i3+1)-u(i1,i2-1,i3-1) )
409 !
410 ! Dt( a31 Dr ) = D0t( a31 D0r ) = a31zzp*( u(i1+1,i2,i3+1)-u(i1-1,i2,i3+1)) - a31zzm*( u(i1+1,i2,i3-1)-u(i1-1,i2,i3-1) )
411 ! Dt( a32 Ds ) = D0t( a32 D0s ) = a32zzp*( u(i1,i2+1,i3+1)-u(i1,i2-1,i3+1)) - a32zzm*( u(i1,i2+1,i3-1)-u(i1,i2-1,i3-1) )
412 !
413 ! Macro args:
414 ! cmp,eqn : component and equation
415 ! A : generic name for coefficients
416 ! =================================================================================================================================================
417 #beginMacro setDivTensorGradCoeff3d(cmp,eqn,A)
418  coeff(MCE(-1,-1,-1,cmp,eqn),i1,i2,i3)= 0.
419  coeff(MCE( 0,-1,-1,cmp,eqn),i1,i2,i3)= A ## 23zmz+A ## 32zzm
420  coeff(MCE( 1,-1,-1,cmp,eqn),i1,i2,i3)= 0.
421  coeff(MCE(-1, 0,-1,cmp,eqn),i1,i2,i3)= A ## 13mzz+A ## 31zzm
422  coeff(MCE( 0, 0,-1,cmp,eqn),i1,i2,i3)= A ## 33mh
423  coeff(MCE( 1, 0,-1,cmp,eqn),i1,i2,i3)= -A ## 13pzz-A ## 31zzm
424  coeff(MCE(-1, 1,-1,cmp,eqn),i1,i2,i3)= 0.
425  coeff(MCE( 0, 1,-1,cmp,eqn),i1,i2,i3)= -A ## 23zpz-A ## 32zzm
426  coeff(MCE( 1, 1,-1,cmp,eqn),i1,i2,i3)= 0.
427 
428  coeff(MCE(-1,-1, 0,cmp,eqn),i1,i2,i3)= A ## 12mzz+A ## 21zmz
429  coeff(MCE( 0,-1, 0,cmp,eqn),i1,i2,i3)= A ## 22mh
430  coeff(MCE( 1,-1, 0,cmp,eqn),i1,i2,i3)= -A ## 12pzz-A ## 21zmz
431  coeff(MCE(-1, 0, 0,cmp,eqn),i1,i2,i3)= A ## 11mh
432  coeff(MCE( 0, 0, 0,cmp,eqn),i1,i2,i3)= -A ## 11ph-A ## 11mh -A ## 22ph -A ## 22mh -A ## 33ph -A ## 33mh
433  coeff(MCE( 1, 0, 0,cmp,eqn),i1,i2,i3)= A ## 11ph
434  coeff(MCE(-1, 1, 0,cmp,eqn),i1,i2,i3)= -A ## 12mzz-A ## 21zpz
435  coeff(MCE( 0, 1, 0,cmp,eqn),i1,i2,i3)= A ## 22ph
436  coeff(MCE( 1, 1, 0,cmp,eqn),i1,i2,i3)= A ## 12pzz+A ## 21zpz
437 
438  coeff(MCE(-1,-1, 1,cmp,eqn),i1,i2,i3)= 0.
439  coeff(MCE( 0,-1, 1,cmp,eqn),i1,i2,i3)= -A ## 23zmz-A ## 32zzp
440  coeff(MCE( 1,-1, 1,cmp,eqn),i1,i2,i3)= 0.
441  coeff(MCE(-1, 0, 1,cmp,eqn),i1,i2,i3)= -A ## 13mzz-A ## 31zzp
442  coeff(MCE( 0, 0, 1,cmp,eqn),i1,i2,i3)= A ## 33ph
443  coeff(MCE( 1, 0, 1,cmp,eqn),i1,i2,i3)= A ## 13pzz+A ## 31zzp
444  coeff(MCE(-1, 1, 1,cmp,eqn),i1,i2,i3)= 0.
445  coeff(MCE( 0, 1, 1,cmp,eqn),i1,i2,i3)= A ## 23zpz+A ## 32zzp
446  coeff(MCE( 1, 1, 1,cmp,eqn),i1,i2,i3)= 0.
447 #endMacro
448 
449 
450 ! =======================================================================================================
451 ! This macro scaled the coefficients that appear in the discretization of the div(tensor grad) operator
452 ! =======================================================================================================
453 #beginMacro scaleCoefficients3d( a11ph,a11mh,a22ph,a22mh,a33ph,a33mh,a12pzz,a12mzz,a13pzz,a13mzz,a21zpz,a21zmz,a23zpz,a23zmz,a31zzp,a31zzm,a32zzp,a32zzm )
454  a11ph=a11ph*dr0i
455  a11mh=a11mh*dr0i
456  a22ph=a22ph*dr1i
457  a22mh=a22mh*dr1i
458  a33ph=a33ph*dr2i
459  a33mh=a33mh*dr2i
460 
461  a12pzz=a12pzz*dr0dr1
462  a12mzz=a12mzz*dr0dr1
463  a13pzz=a13pzz*dr0dr2
464  a13mzz=a13mzz*dr0dr2
465 
466  a21zpz=a21zpz*dr0dr1
467  a21zmz=a21zmz*dr0dr1
468  a23zpz=a23zpz*dr1dr2
469  a23zmz=a23zmz*dr1dr2
470 
471  a31zzp=a31zzp*dr0dr2
472  a31zzm=a31zzm*dr0dr2
473  a32zzp=a32zzp*dr1dr2
474  a32zzm=a32zzm*dr1dr2
475 
476 #endMacro
Generated on Fri Jan 4 2013 10:19:17 for CG by   doxygen 1.8.3