CG  Version 25
bcOptSmFOS3D.h
Go to the documentation of this file.
1 c*******
2 c******* Fill in forcing arrays if they are not provided ***********
3 c*******
4  beginLoopOverSides( numGhost,numGhost )
5  if( boundaryCondition(side,axis).eq.displacementBC )then
6  if( addBoundaryForcing(side,axis).eq.0 ) then
7  beginGhostLoopsMask3d() ! *wdh* we can assign dirichlet like conditions at ghost/imter points too
8  ! given displacements
9  bcfa(side,axis,i1,i2,i3,uc) = 0.0
10  bcfa(side,axis,i1,i2,i3,vc) = 0.0
11  bcfa(side,axis,i1,i2,i3,wc) = 0.0
12  ! given velocities
13  bcfa(side,axis,i1,i2,i3,v1c) = 0.0
14  bcfa(side,axis,i1,i2,i3,v2c) = 0.0
15  bcfa(side,axis,i1,i2,i3,v3c) = 0.0
16  ! given acceleration
17  bcfa(side,axis,i1,i2,i3,s11c) = 0.0
18  bcfa(side,axis,i1,i2,i3,s12c) = 0.0
19  bcfa(side,axis,i1,i2,i3,s13c) = 0.0
20  endLoopsMask3d()
21  else if( twilightZone.ne.0 ) then
22  beginGhostLoopsMask3d()
23  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,ue )
24  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,ve )
25  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,we )
26 
27  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1e )
28  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2e )
29  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3e )
30 
31  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11x )
32  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21y )
33  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31z )
34 
35  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12x )
36  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22y )
37  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32z )
38 
39  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13x )
40  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23y )
41  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33z )
42 
43  bcfa(side,axis,i1,i2,i3,uc) = ue
44  bcfa(side,axis,i1,i2,i3,vc) = ve
45  bcfa(side,axis,i1,i2,i3,wc) = we
46 
47  bcfa(side,axis,i1,i2,i3,v1c) = v1e
48  bcfa(side,axis,i1,i2,i3,v2c) = v2e
49  bcfa(side,axis,i1,i2,i3,v3c) = v3e
50 
51  bcfa(side,axis,i1,i2,i3,s11c) = tau11x+tau21y+tau31z
52  bcfa(side,axis,i1,i2,i3,s12c) = tau12x+tau22y+tau32z
53  bcfa(side,axis,i1,i2,i3,s13c) = tau13x+tau23y+tau33z
54  endLoopsMask3d()
55  end if
56  else if( boundaryCondition(side,axis).eq.tractionBC ) then
57  if( addBoundaryForcing(side,axis).eq.0 )then
58  beginGhostLoopsMask3d()
59  ! given traction (for the traction BC)
60  bcfa(side,axis,i1,i2,i3,s11c) = 0.0
61  bcfa(side,axis,i1,i2,i3,s12c) = 0.0
62  bcfa(side,axis,i1,i2,i3,s13c) = 0.0
63 
64  ! given traction (for determining displacements). Normally this is equal to the above
65  ! traction values except when using twilight-zone
66  bcfa(side,axis,i1,i2,i3,uc) = 0.0
67  bcfa(side,axis,i1,i2,i3,vc) = 0.0
68  bcfa(side,axis,i1,i2,i3,wc) = 0.0
69 
70  ! given rate of change of traction (for determining the velocity)
71  bcfa(side,axis,i1,i2,i3,v1c) = 0.0
72  bcfa(side,axis,i1,i2,i3,v2c) = 0.0
73  bcfa(side,axis,i1,i2,i3,v3c) = 0.0
74  endLoopsMask3d()
75  else if( twilightZone.ne.0 )then
76  beginGhostLoopsMask3d()
77  ! (an1,an2,an3) = outward normal
78  if( gridType.eq.rectangular )then
79  if( axis.eq.0 )then
80  an1 = -is
81  an2 = 0.0
82  an3 = 0.0
83  else if( axis.eq.1 ) then
84  an1 = 0.0
85  an2 = -is
86  an3 = 0.0
87  else
88  an1 = 0.0
89  an2 = 0.0
90  an3 = -is
91  end if
92  else
93  aNormi = 1.0/max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
94  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
95  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
96  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
97  end if
98  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
99  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
100  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
101 
102  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
103  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
104  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
105 
106  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
107  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
108  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
109 
110  bcfa(side,axis,i1,i2,i3,uc) = an1*( kappa*u1x+lambda*(u2y+u3z) ) + an2*( mu*(u2x+u1y) ) + an3*( mu*(u3x+u1z) )
111  bcfa(side,axis,i1,i2,i3,vc) = an1*( mu*(u2x+u1y) ) + an2*( kappa*u2y+lambda*(u1x+u3z) ) + an3*( mu*(u3y+u2z) )
112  bcfa(side,axis,i1,i2,i3,wc) = an1*( mu*(u3x+u1z) ) + an2*( mu*(u3y+u2z) ) + an3*( kappa*u3z+lambda*(u1x+u2y) )
113 
114  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1x )
115  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1y )
116  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1z )
117 
118  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2x )
119  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2y )
120  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2z )
121 
122  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3x )
123  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3y )
124  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3z )
125 
126  bcfa(side,axis,i1,i2,i3,v1c) = an1*( kappa*v1x+lambda*(v2y+v3z) ) + an2*( mu*(v2x+v1y) ) + an3*( mu*(v3x+v1z) )
127  bcfa(side,axis,i1,i2,i3,v2c) = an1*( mu*(v2x+v1y) ) + an2*( kappa*v2y+lambda*(v1x+v3z) ) + an3*( mu*(v3y+v2z) )
128  bcfa(side,axis,i1,i2,i3,v3c) = an1*( mu*(v3x+v1z) ) + an2*( mu*(v3y+v2z) ) + an3*( kappa*v3z+lambda*(v1x+v2y) )
129 
130  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11 )
131  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21 )
132  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31 )
133  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12 )
134  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22 )
135  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32 )
136  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13 )
137  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23 )
138  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33 )
139 
140  ! note : n_j sigma_ji : sum over first index
141  bcfa(side,axis,i1,i2,i3,s11c) = an1*tau11+an2*tau21+an3*tau31
142  bcfa(side,axis,i1,i2,i3,s12c) = an1*tau12+an2*tau22+an3*tau32
143  bcfa(side,axis,i1,i2,i3,s13c) = an1*tau13+an2*tau23+an3*tau33
144 
145  endLoopsMask3d()
146 
147  else
148  ! fill in the traction BC into the stress components
149  ! (this is needed since for TZ flow these values are different)
150  beginGhostLoopsMask3d()
151  bcfa(side,axis,i1,i2,i3,s11c) = bcfa(side,axis,i1,i2,i3,uc)
152  bcfa(side,axis,i1,i2,i3,s12c) = bcfa(side,axis,i1,i2,i3,vc)
153  bcfa(side,axis,i1,i2,i3,s13c) = bcfa(side,axis,i1,i2,i3,wc)
154  endLoopsMask3d()
155  end if
156 
157  else if( boundaryCondition(side,axis).eq.slipWall ) then
158  if( addBoundaryForcing(side,axis).eq.0 ) then
159  beginGhostLoopsMask3d()
160  !! check these components with Bill ... FIX ME!! ...
161  ! given tangential stresses (often zero)
162  bcfa(side,axis,i1,i2,i3,s11c) = 0.0
163  bcfa(side,axis,i1,i2,i3,s12c) = 0.0
164 
165  ! time rate of change of tangential stresses
166  bcfa(side,axis,i1,i2,i3,s21c) = 0.0
167  bcfa(side,axis,i1,i2,i3,s22c) = 0.0
168 
169  ! given normal displacement
170  bcfa(side,axis,i1,i2,i3,uc) = 0.0
171 
172  ! time rate of change of normal displacement
173  bcfa(side,axis,i1,i2,i3,v1c) = 0.0
174  endLoopsMask3d()
175  else if( twilightZone.ne.0 ) then
176  beginGhostLoopsMask3d()
177  ! (an1,an2,an3) = outward normal
178  if( gridType.eq.rectangular ) then
179  an1 = 0.0
180  an2 = 0.0
181  an3 = 0.0
182 
183  sn1 = 0.0
184  sn2 = 0.0
185  sn3 = 0.0
186 
187  tn1 = 0.0
188  tn2 = 0.0
189  tn3 = 0.0
190  if( axis.eq.0 ) then
191  an1 = -is
192  sn2 = -is
193  tn3 = -is
194  else if( axis.eq.1 ) then
195  an2 = -is
196  sn1 = -is
197  tn3 = -is
198  else
199  an3 = -is
200  sn1 = -is
201  tn2 = -is
202  end if
203  else
204  if( axis.eq.0 ) then
205  tan1c = 1
206  tan2c = 2
207  else if( axis.eq.1 ) then
208  tan1c = 0
209  tan2c = 2
210  else
211  tan1c = 0
212  tan2c = 1
213  end if
214  aNormi = 1.0/max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
215  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
216  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
217  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
218 
219  sn1 = rx(i1,i2,i3,tan1c,0)
220  sn2 = rx(i1,i2,i3,tan1c,1)
221  sn3 = rx(i1,i2,i3,tan1c,2)
222 
223  tn1 = rx(i1,i2,i3,tan2c,0)
224  tn2 = rx(i1,i2,i3,tan2c,1)
225  tn3 = rx(i1,i2,i3,tan2c,2)
226 
227  ! set sn to be part of sn which is orthogonal to an
228  alpha = an1*sn1+an2*sn2+an3*sn3
229  sn1 = sn1-alpha*an1
230  sn2 = sn2-alpha*an2
231  sn3 = sn3-alpha*an3
232  ! normalize sn
233  aNormi = 1.0/max(epsx,sqrt(sn1**2+sn2**2+sn3**2))
234  sn1 = sn1*aNormi
235  sn2 = sn2*aNormi
236  sn3 = sn3*aNormi
237 
238  ! set tn to be part of tn which is orthogonal to an and sn
239  alpha = an1*tn1+an2*tn2+an3*tn3
240  tn1 = tn1-alpha*an1
241  tn2 = tn2-alpha*an2
242  tn3 = tn3-alpha*an3
243  alpha = sn1*tn1+sn2*tn2+sn3*tn3
244  tn1 = tn1-alpha*sn1
245  tn2 = tn2-alpha*sn2
246  tn3 = tn3-alpha*sn3
247  ! normalize tn
248  aNormi = 1.0/max(epsx,sqrt(tn1**2+tn2**2+tn3**2))
249  tn1 = tn1*aNormi
250  tn2 = tn2*aNormi
251  tn3 = tn3*aNormi
252  end if
253 
254  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,ue)
255  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,ve)
256  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,we)
257 
258  bcfa(side,axis,i1,i2,i3,uc) = an1*ue+an2*ve+an3*we
259 
260  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,ue)
261  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,ve)
262  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,we)
263 
264  bcfa(side,axis,i1,i2,i3,v1c) = an1*ue+an2*ve+an3*we
265 
266  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11)
267  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21)
268  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31)
269  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12)
270  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22)
271  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32)
272  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13)
273  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23)
274  call ogDeriv(ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33)
275 
276  ! check indicies ... FIX ME!! ...
277  bcfa(side,axis,i1,i2,i3,s11c) = an1*(sn1*tau11+sn2*tau12+sn3*tau13)+ \
278  an2*(sn1*tau21+sn2*tau22+sn3*tau23)+ \
279  an3*(sn1*tau31+sn2*tau32+sn3*tau33)
280  bcfa(side,axis,i1,i2,i3,s12c) = an1*(tn1*tau11+tn2*tau12+tn3*tau13)+ \
281  an2*(tn1*tau21+tn2*tau22+tn3*tau23)+ \
282  an3*(tn1*tau31+tn2*tau32+tn3*tau33)
283 
284  call ogDeriv(ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x)
285  call ogDeriv(ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y)
286  call ogDeriv(ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z)
287  call ogDeriv(ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x)
288  call ogDeriv(ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y)
289  call ogDeriv(ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z)
290  call ogDeriv(ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x)
291  call ogDeriv(ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y)
292  call ogDeriv(ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z)
293 
294  tau11 = kappa*u1x+lambda*(u2y+u3z)
295  tau21 = mu*(u2x+u1y)
296  tau31 = mu*(u3x+u1z)
297  tau12 = tau21
298  tau22 = kappa*u2y+lambda*(u1x+u3z)
299  tau32 = mu*(u3y+u2z)
300  tau13 = tau31
301  tau23 = tau32
302  tau33 = kappa*u3z+lambda*(u1x+u2y)
303 ccccccccccccccccccccccccccccccccccccccccccccccccccc
304 
305  bcfa(side,axis,i1,i2,i3,s21c) = an1*(-an2*tau11+an1*tau12)+an2*(-an2*tau21+an1*tau22)
306 
307  call ogDeriv(ep,1,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11)
308  call ogDeriv(ep,1,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21)
309  call ogDeriv(ep,1,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12)
310  call ogDeriv(ep,1,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22)
311 
312  bcfa(side,axis,i1,i2,i3,s12c) = an1*(-an2*tau11+an1*tau12)+an2*(-an2*tau21+an1*tau22)
313 
314  call ogDeriv(ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1ex)
315  call ogDeriv(ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1ey)
316  call ogDeriv(ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2ex)
317  call ogDeriv(ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2ey)
318 
319  tau11=(lambda+2.*mu)*v1ex+lambda*v2ey
320  tau12=mu*(v1ey+v2ex)
321  tau21=tau12
322  tau22=lambda*v1ex+(lambda+2.*mu)*v2ey
323  bcfa(side,axis,i1,i2,i3,s22c) = an1*(-an2*tau11+an1*tau12)+an2*(-an2*tau21+an1*tau22)
324 
325  endLoopsMask3d()
326  else
327  ! fill in the traction BC into the stress components *wdh* 081109
328  ! (this is needed since for TZ flow these values are different)
329  beginGhostLoopsMask3d()
330  bcfa(side,axis,i1,i2,i3,s11c) = bcfa(side,axis,i1,i2,i3,uc)
331  bcfa(side,axis,i1,i2,i3,s12c) = bcfa(side,axis,i1,i2,i3,vc)
332  endLoopsMask3d()
333  end if
334 
335  else if( boundaryCondition(side,axis).gt.0 .and. boundaryCondition(side,axis).ne.dirichletBoundaryCondition ) then
336  write(*,'("smg3d:BC: unknown BC: side,axis,grid, boundaryCondition=",i2,i2,i4,i8)') side,axis,grid,boundaryCondition(side,axis)
337  end if
338  endLoopOverSides()
339 
340 c*******
341 c******* Primary Dirichlet boundary conditions ***********
342 c*******
343  beginLoopOverSides( numGhost,numGhost )
344  if( boundaryCondition(side,axis).eq.displacementBC )then
345  ! ..step 0: Dirichlet bcs for displacement and velocity
346  beginGhostLoopsMask3d()
347  ! given displacements
348  u(i1,i2,i3,uc) = bcf(side,axis,i1,i2,i3,uc)
349  u(i1,i2,i3,vc) = bcf(side,axis,i1,i2,i3,vc)
350  u(i1,i2,i3,wc) = bcf(side,axis,i1,i2,i3,wc)
351 
352  ! given velocities
353  u(i1,i2,i3,v1c) = bcf(side,axis,i1,i2,i3,v1c)
354  u(i1,i2,i3,v2c) = bcf(side,axis,i1,i2,i3,v2c)
355  u(i1,i2,i3,v3c) = bcf(side,axis,i1,i2,i3,v3c)
356  endLoopsMask3d()
357  else if( boundaryCondition(side,axis).eq.tractionBC )then
358  ! dirichlet portion of traction BC
359  if( gridType.eq.rectangular ) then
360  if( axis.eq.0 )then
361  beginGhostLoopsMask3d()
362  ! set normal components of the stress, n=(-is,0,0)
363  u(i1,i2,i3,s11c) = -is*bcf(side,axis,i1,i2,i3,s11c)
364  u(i1,i2,i3,s12c) = -is*bcf(side,axis,i1,i2,i3,s12c)
365  u(i1,i2,i3,s13c) = -is*bcf(side,axis,i1,i2,i3,s13c)
366  endLoopsMask3d()
367  else if( axis.eq.1 ) then
368  beginGhostLoopsMask3d()
369  ! set normal components of the stress, n=(0,-is,0)
370  u(i1,i2,i3,s21c) = -is*bcf(side,axis,i1,i2,i3,s11c)
371  u(i1,i2,i3,s22c) = -is*bcf(side,axis,i1,i2,i3,s12c)
372  u(i1,i2,i3,s23c) = -is*bcf(side,axis,i1,i2,i3,s13c)
373  endLoopsMask3d()
374  else
375  beginGhostLoopsMask3d()
376  ! set normal components of the stress, n=(0,0,-is)
377  u(i1,i2,i3,s31c) = -is*bcf(side,axis,i1,i2,i3,s11c)
378  u(i1,i2,i3,s32c) = -is*bcf(side,axis,i1,i2,i3,s12c)
379  u(i1,i2,i3,s33c) = -is*bcf(side,axis,i1,i2,i3,s13c)
380  endLoopsMask3d()
381  end if
382 
383  else ! curvilinear
384  beginGhostLoopsMask3d()
385  f1 = bcf(side,axis,i1,i2,i3,s11c)
386  f2 = bcf(side,axis,i1,i2,i3,s12c)
387  f3 = bcf(side,axis,i1,i2,i3,s13c)
388 
389  ! (an1,an2,an3) = outward normal
390  aNormi = 1.0/max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
391  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
392  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
393  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
394 
395  b1 = f1-(an1*u(i1,i2,i3,s11c)+an2*u(i1,i2,i3,s21c)+an3*u(i1,i2,i3,s31c))
396  b2 = f2-(an1*u(i1,i2,i3,s12c)+an2*u(i1,i2,i3,s22c)+an3*u(i1,i2,i3,s32c))
397  b3 = f3-(an1*u(i1,i2,i3,s13c)+an2*u(i1,i2,i3,s23c)+an3*u(i1,i2,i3,s33c))
398 
399  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)+an1*b1
400  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)+an1*b2
401  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)+an1*b3
402 
403  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)+an2*b1
404  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)+an2*b2
405  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)+an2*b3
406 
407  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)+an3*b1
408  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)+an3*b2
409  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)+an3*b3
410  endLoopsMask3d()
411  end if ! grid type
412 
413  else if( boundaryCondition(side,axis).eq.slipWall ) then
414  ! ********* SlipWall BC ********
415  ! set "dirichlet" parts of the slipwall BC
416  if( gridType.eq.rectangular ) then
417  if( axis.eq.0 ) then
418  beginGhostLoopsMask3d()
419  ! set n.tau.t and the normal component of displacement, n=(-is,0,0)
420  u(i1,i2,i3,s12c) = bcf(side,axis,i1,i2,i3,s11c)
421  u(i1,i2,i3,s13c) = bcf(side,axis,i1,i2,i3,s12c)
422  u(i1,i2,i3,uc) = -is*bcf(side,axis,i1,i2,i3,uc)
423  u(i1,i2,i3,v1c) = -is*bcf(side,axis,i1,i2,i3,v1c)
424  endLoopsMask3d()
425  else if( axis.eq.1 ) then
426  beginGhostLoopsMask3d()
427  ! set n.tau.t and the normal component of displacement, n=(-is,0,0)
428  u(i1,i2,i3,s21c) = bcf(side,axis,i1,i2,i3,s11c)
429  u(i1,i2,i3,s23c) = bcf(side,axis,i1,i2,i3,s12c)
430  u(i1,i2,i3,vc) = -is*bcf(side,axis,i1,i2,i3,uc)
431  u(i1,i2,i3,v2c) = -is*bcf(side,axis,i1,i2,i3,v1c)
432  endLoopsMask3d()
433  else
434  beginGhostLoopsMask3d()
435  ! set n.tau.t and the normal component of displacement, n=(-is,0,0)
436  u(i1,i2,i3,s31c) = bcf(side,axis,i1,i2,i3,s11c)
437  u(i1,i2,i3,s32c) = bcf(side,axis,i1,i2,i3,s12c)
438  u(i1,i2,i3,wc) = -is*bcf(side,axis,i1,i2,i3,uc)
439  u(i1,i2,i3,v3c) = -is*bcf(side,axis,i1,i2,i3,v1c)
440  endLoopsMask3d()
441  end if
442 
443  else ! curvilinear
444  beginGhostLoopsMask3d()
445  ! given tangential traction forces
446  f1 = bcf(side,axis,i1,i2,i3,s11c)
447  f2 = bcf(side,axis,i1,i2,i3,s12c)
448 
449  ! (an1,an2,an3) = outward normal
450  if( axis.eq.0 ) then
451  tan1c = 1
452  tan2c = 2
453  else if( axis.eq.1 ) then
454  tan1c = 0
455  tan2c = 2
456  else
457  tan1c = 0
458  tan2c = 1
459  end if
460  aNormi = 1.0/max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
461  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
462  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
463  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
464 
465  sn1 = rx(i1,i2,i3,tan1c,0)
466  sn2 = rx(i1,i2,i3,tan1c,1)
467  sn3 = rx(i1,i2,i3,tan1c,2)
468 
469  tn1 = rx(i1,i2,i3,tan2c,0)
470  tn2 = rx(i1,i2,i3,tan2c,1)
471  tn3 = rx(i1,i2,i3,tan2c,2)
472 
473  ! set sn to be part of sn which is orthogonal to an
474  alpha = an1*sn1+an2*sn2+an3*sn3
475  sn1 = sn1-alpha*an1
476  sn2 = sn2-alpha*an2
477  sn3 = sn3-alpha*an3
478  ! normalize sn
479  aNormi = 1.0/max(epsx,sqrt(sn1**2+sn2**2+sn3**2))
480  sn1 = sn1*aNormi
481  sn2 = sn2*aNormi
482  sn3 = sn3*aNormi
483 
484  ! set tn to be part of tn which is orthogonal to an and sn
485  alpha = an1*tn1+an2*tn2+an3*tn3
486  tn1 = tn1-alpha*an1
487  tn2 = tn2-alpha*an2
488  tn3 = tn3-alpha*an3
489  alpha = sn1*tn1+sn2*tn2+sn3*tn3
490  tn1 = tn1-alpha*sn1
491  tn2 = tn2-alpha*sn2
492  tn3 = tn3-alpha*sn3
493  ! normalize tn
494  aNormi = 1.0/max(epsx,sqrt(tn1**2+tn2**2+tn3**2))
495  tn1 = tn1*aNormi
496  tn2 = tn2*aNormi
497  tn3 = tn3*aNormi
498 
499  b1 = f1-an1*(u(i1,i2,i3,s11c)*sn1+u(i1,i2,i3,s12c)*sn2+u(i1,i2,i3,s13c)*sn3)- \
500  an2*(u(i1,i2,i3,s21c)*sn1+u(i1,i2,i3,s22c)*sn2+u(i1,i2,i3,s23c)*sn3)- \
501  an3*(u(i1,i2,i3,s31c)*sn1+u(i1,i2,i3,s32c)*sn2+u(i1,i2,i3,s33c)*sn3)
502  b2 = f2-an1*(u(i1,i2,i3,s11c)*tn1+u(i1,i2,i3,s12c)*tn2+u(i1,i2,i3,s13c)*tn3)- \
503  an2*(u(i1,i2,i3,s21c)*tn1+u(i1,i2,i3,s22c)*tn2+u(i1,i2,i3,s23c)*tn3)- \
504  an3*(u(i1,i2,i3,s31c)*tn1+u(i1,i2,i3,s32c)*tn2+u(i1,i2,i3,s33c)*tn3)
505 
506 
507  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)+an1*b1*sn1+an1*b2*tn1
508  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)+an1*b1*sn2+an1*b2*tn2
509  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)+an1*b1*sn3+an1*b2*tn3
510 
511  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)+an2*b1*sn1+an2*b2*tn1
512  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)+an2*b1*sn2+an2*b2*tn2
513  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)+an2*b1*sn3+an2*b2*tn3
514 
515  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)+an3*b1*sn1+an3*b2*tn1
516  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)+an3*b1*sn2+an3*b2*tn2
517  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)+an3*b1*sn3+an3*b2*tn3
518 
519 
520  ! given normal displacement
521  f1 = bcf(side,axis,i1,i2,i3,uc)
522 
523  ! given normal velocity
524  f2 = bcf(side,axis,i1,i2,i3,v1c)
525 
526  b1 = f1-an1*u(i1,i2,i3,uc)-an2*u(i1,i2,i3,vc)-an3*u(i1,i2,i3,wc)
527  b2 = f2-an1*u(i1,i2,i3,v1c)-an2*u(i1,i2,i3,v2c)-an3*u(i1,i2,i3,v3c)
528 
529  u(i1,i2,i3,uc) = u(i1,i2,i3,uc)+an1*b1
530  u(i1,i2,i3,vc) = u(i1,i2,i3,vc)+an2*b1
531  u(i1,i2,i3,wc) = u(i1,i2,i3,wc)+an3*b1
532 
533  u(i1,i2,i3,v1c) = u(i1,i2,i3,v1c)+an1*b2
534  u(i1,i2,i3,v2c) = u(i1,i2,i3,v2c)+an2*b2
535  u(i1,i2,i3,v3c) = u(i1,i2,i3,v3c)+an3*b2
536  endLoopsMask3d()
537 
538  end if ! end gridType
539 
540  else if( boundaryCondition(side,axis).gt.0 .and. boundaryCondition(side,axis).ne.dirichletBoundaryCondition ) then
541  write(*,'("smg3d:BC: unknown BC: side,axis,grid, boundaryCondition=",i2,i2,i4,i8)') side,axis,grid,boundaryCondition(side,axis)
542 
543  end if ! bc type
544  endLoopOverSides()
545 
546 c*******
547 c******* Extrapolate to the first ghost cells (only for physical sides) ********
548 c*******
549  ! *wdh* For now assign 2 ghost lines and points outside edges and corners
550  extra = 2
551  numGhostExtrap=2
552  beginLoopOverSides( extra,numGhostExtrap )
553  if( boundaryCondition(side,axis).gt.0.and.boundaryCondition(side,axis).ne.dirichletBoundaryCondition )then
554 
555  if( .false. )then
556  write(*,'(" bcOpt: Extrap ghost: grid,side,axis=",3i3,", \
557  loop bounds: nn1a,nn1b,nn2a,nn2b,nn3a,nn3b=",6i3)') grid,side,axis,\
558  nn1a,nn1b,nn2a,nn2b,nn3a,nn3b
559 
560  end if
561 
562  beginGhostLoops3d()
563  if( mask(i1,i2,i3).ne.0 ) then
564  do n=0,numberOfComponents-1
565  u(i1-is1,i2-is2,i3-is3,n)=extrap3(u,i1,i2,i3,n,is1,is2,is3)
566  u(i1-2*is1,i2-2*is2,i3-2*is3,n)=extrap3(u,i1-is1,i2-is2,i3-is3,n,is1,is2,is3)
567  end do
568  end if
569  endLoops3d()
570  end if
571  endLoopOverSides()
572 
573 c*******
574 c******* Fix up components of stress along the edges
575 c*******
576  #Include 'bcOptSmFOS3DEdge.h'
577 
578 cccccccccccccccccccccccccccccccccccccccccccccccccc
579 c .. set exact solution for corners for now
580  if( setCornersWithTZ .and. twilightZone.ne.0 ) then ! *wdh* 090909
581  write(*,'(" bcOptSmFOS3D: INFO set exact values on corners")')
582  beginLoopOverCorners3d(0)
583  if( mask(i1,i2,i3).gt.0 ) then
584  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,ue )
585  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,ve )
586  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,we )
587 
588  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1e )
589  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2e )
590  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3e )
591 
592  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11e )
593  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21e )
594  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31e )
595 
596  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12e )
597  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22e )
598  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32e )
599 
600  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13e )
601  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23e )
602  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33e )
603 
604  u(i1,i2,i3,uc) = ue
605  u(i1,i2,i3,vc) = ve
606  u(i1,i2,i3,wc) = we
607 
608  u(i1,i2,i3,v1c) = v1e
609  u(i1,i2,i3,v2c) = v2e
610  u(i1,i2,i3,v3c) = v3e
611 
612  u(i1,i2,i3,s11c) = tau11e
613  u(i1,i2,i3,s21c) = tau21e
614  u(i1,i2,i3,s31c) = tau31e
615 
616  u(i1,i2,i3,s12c) = tau12e
617  u(i1,i2,i3,s22c) = tau22e
618  u(i1,i2,i3,s32c) = tau32e
619 
620  u(i1,i2,i3,s13c) = tau13e
621  u(i1,i2,i3,s23c) = tau23e
622  u(i1,i2,i3,s33c) = tau33e
623  end if
624  endLoopOverCorners3d()
625  end if
626 cccccccccccccccccccccccccccccccccccccccccccccccccc
627 
628 c*******
629 c******* Fix up components of stress in the corners
630 c*******
631 c.. for now we are going to ignore this too
632 
633 c*******
634 c******* Secondary Neumann boundary conditions (compatibility conditions) ********
635 c*******
636  beginLoopOverSides(numGhost,numGhost)
637 
638  if( boundaryCondition(side,axis).eq.displacementBC ) then
639 
640  if( gridType.eq.rectangular ) then
641 
642  ! ********* DISPLACEMENT : Cartesian Grid **********
643 
644  ! Use momentum equations ...
645  ! s11_x + s21_y + s31_z = rho * u_tt
646  ! s12_x + s22_y + s32_z = rho * v_tt
647  ! s13_x + s23_y + s33_z = rho * w_tt
648  ! *wdh* 090909 -- only assign pts where mask > 0 since we assume values at adjacent points.
649  beginLoopsMask3d()
650  accel(1) = rho*bcf(side,axis,i1,i2,i3,s11c)
651  accel(2) = rho*bcf(side,axis,i1,i2,i3,s12c)
652  accel(3) = rho*bcf(side,axis,i1,i2,i3,s13c)
653 
654 c write(6,*)is1,is2,is3,is,axis
655 c write(6,*)accel(1),accel(2),accel(3)
656 
657  do isc = 1,3
658  u(i1-is1,i2-is2,i3-is3,sc(axis+1,isc)) = u(i1+is1,i2+is2,i3+is3,sc(axis+1,isc))-2.0*is*dx(axis)*(accel(isc)- \
659  (1.0-delta(axis+1,1))*(u(i1+1,i2,i3,sc(1,isc))-u(i1-1,i2,i3,sc(1,isc)))/(2.0*dx(0))- \
660  (1.0-delta(axis+1,2))*(u(i1,i2+1,i3,sc(2,isc))-u(i1,i2-1,i3,sc(2,isc)))/(2.0*dx(1))- \
661  (1.0-delta(axis+1,3))*(u(i1,i2,i3+1,sc(3,isc))-u(i1,i2,i3-1,sc(3,isc)))/(2.0*dx(2)))
662  end do
663  endLoopsMask3d()
664  else
665 
666  if( .false. ) then ! non-free stream preserving method
667  ! *********** DISPLACEMENT : Curvilinear Grid (not free stream preserving) ****************
668 
669  ! Use momentum equations to get J*(rx,ry,rz).(s11,s21,s31)(-1) = s11tilde
670  ! (1) D_r1[ J*(rx,ry,rz).(s11,s21,s31)] + D_r2[J*(sx,sy,sz).(s11,s21,s31)] + D_r3[J*(tx,ty,tz).(s11,s21,s31)] = J * rho * u_tt
671  ! s11tilde s21tilde s31tilde
672  ! (2) Use extrapolated values to get J*(sx,sy,sz).(s11,s21,s31)(-1) = s21tilde
673  ! (3) Use extrapolated values to get J*(tx,ty,tz).(s11,s21,s31)(-1) = s31tilde
674  ! To give 3 equations for (s11,s21,s31) on the ghost point:
675  ! (J rx) s11(-1) + (J ry) s21(-1) + (J rz) s31(-1) = f1 = s11tilde (from momentum eqn)
676  ! (J sx) s11(-1) + (J sy) s21(-1) + (J sz) s31(-1) = f2 = s21tilde (from extrapolated values)
677  ! (J tx) s11(-1) + (J ty) s21(-1) + (J tz) s31(-1) = f3 = s31tilde (from extrapolated values)
678  ! Solve: (note the Jacobian cancels when the matrix inversion is determined)
679  ! s11(-1) = (sy*tz-sz*ty)*f1 + (rz*ty-ry*tz)*f2 + (ry*sz-rz*sy)*f3
680  ! s21(-1) = (tx*sz-sx*tz)*f1 + (rx*tz-rz*tx)*f2 + (rz*sx-rx*sz)*f3
681  ! s31(-1) = (sx*ty-sy*tx)*f1 + (ry*tx-rx*ty)*f2 + (rx*sy-ry*sx)*f3
682  !
683  ! (A similar expression holds for other stresses)
684  beginLoopsMask3d()
685  accel(1) = rho*bcf(side,axis,i1,i2,i3,s11c)
686  accel(2) = rho*bcf(side,axis,i1,i2,i3,s12c)
687  accel(3) = rho*bcf(side,axis,i1,i2,i3,s13c)
688 
689  met(0,0) = rx(i1,i2,i3,0,0)
690  met(0,1) = rx(i1,i2,i3,0,1)
691  met(0,2) = rx(i1,i2,i3,0,2)
692  met(1,0) = rx(i1,i2,i3,1,0)
693  met(1,1) = rx(i1,i2,i3,1,1)
694  met(1,2) = rx(i1,i2,i3,1,2)
695  met(2,0) = rx(i1,i2,i3,2,0)
696  met(2,1) = rx(i1,i2,i3,2,1)
697  met(2,2) = rx(i1,i2,i3,2,2)
698 
699  ! loop over stress components
700  do isc = 1,3
701  do idot = 1,3
702  ! these are extrapolated components
703  stilde(idot) = det(i1-is1,i2-is2,i3-is3)*(rx(i1-is1,i2-is2,i3-is3,idot-1,0)*u(i1-is1,i2-is2,i3-is3,sc(1,isc))+ \
704  rx(i1-is1,i2-is2,i3-is3,idot-1,1)*u(i1-is1,i2-is2,i3-is3,sc(2,isc))+ \
705  rx(i1-is1,i2-is2,i3-is3,idot-1,2)*u(i1-is1,i2-is2,i3-is3,sc(3,isc)))
706  end do
707  ! now override in the direction we are currently looking
708  stilde(axis+1) = det(i1+is1,i2+is2,i3+is3)*(rx(i1+is1,i2+is2,i3+is3,axis,0)*u(i1+is1,i2+is2,i3+is3,sc(1,isc))+ \
709  rx(i1+is1,i2+is2,i3+is3,axis,1)*u(i1+is1,i2+is2,i3+is3,sc(2,isc))+ \
710  rx(i1+is1,i2+is2,i3+is3,axis,2)*u(i1+is1,i2+is2,i3+is3,sc(3,isc)))- \
711  2.0*dr(axis)*is*(det(i1,i2,i3)*accel(isc)- \
712  (1.0-delta(axis+1,1))* \
713  (det(i1+1,i2,i3)*(rx(i1+1,i2,i3,0,0)*u(i1+1,i2,i3,sc(1,isc))+ \
714  rx(i1+1,i2,i3,0,1)*u(i1+1,i2,i3,sc(2,isc))+ \
715  rx(i1+1,i2,i3,0,2)*u(i1+1,i2,i3,sc(3,isc)))- \
716  det(i1-1,i2,i3)*(rx(i1-1,i2,i3,0,0)*u(i1-1,i2,i3,sc(1,isc))+ \
717  rx(i1-1,i2,i3,0,1)*u(i1-1,i2,i3,sc(2,isc))+ \
718  rx(i1-1,i2,i3,0,2)*u(i1-1,i2,i3,sc(3,isc))))/(2.0*dr(0))- \
719  (1.0-delta(axis+1,2))* \
720  (det(i1,i2+1,i3)*(rx(i1,i2+1,i3,1,0)*u(i1,i2+1,i3,sc(1,isc))+ \
721  rx(i1,i2+1,i3,1,1)*u(i1,i2+1,i3,sc(2,isc))+ \
722  rx(i1,i2+1,i3,1,2)*u(i1,i2+1,i3,sc(3,isc)))- \
723  det(i1,i2-1,i3)*(rx(i1,i2-1,i3,1,0)*u(i1,i2-1,i3,sc(1,isc))+ \
724  rx(i1,i2-1,i3,1,1)*u(i1,i2-1,i3,sc(2,isc))+ \
725  rx(i1,i2-1,i3,1,2)*u(i1,i2-1,i3,sc(3,isc))))/(2.0*dr(1))- \
726  (1.0-delta(axis+1,3))* \
727  (det(i1,i2,i3+1)*(rx(i1,i2,i3+1,2,0)*u(i1,i2,i3+1,sc(1,isc))+ \
728  rx(i1,i2,i3+1,2,1)*u(i1,i2,i3+1,sc(2,isc))+ \
729  rx(i1,i2,i3+1,2,2)*u(i1,i2,i3+1,sc(3,isc)))- \
730  det(i1,i2,i3-1)*(rx(i1,i2,i3-1,2,0)*u(i1,i2,i3-1,sc(1,isc))+ \
731  rx(i1,i2,i3-1,2,1)*u(i1,i2,i3-1,sc(2,isc))+ \
732  rx(i1,i2,i3-1,2,2)*u(i1,i2,i3-1,sc(3,isc))))/(2.0*dr(2)))
733 
734  u(i1-is1,i2-is2,i3-is3,sc(1,isc)) = (met(1,1)*met(2,2)-met(1,2)*met(2,1))*stilde(1)+ \
735  (met(0,2)*met(2,1)-met(0,1)*met(2,2))*stilde(2)+ \
736  (met(0,1)*met(1,2)-met(0,2)*met(1,1))*stilde(3)
737  u(i1-is1,i2-is2,i3-is3,sc(2,isc)) = (met(1,2)*met(2,0)-met(1,0)*met(2,2))*stilde(1)+ \
738  (met(0,0)*met(2,2)-met(0,2)*met(2,0))*stilde(2)+ \
739  (met(0,2)*met(1,0)-met(0,0)*met(1,2))*stilde(3)
740  u(i1-is1,i2-is2,i3-is3,sc(3,isc)) = (met(1,0)*met(2,1)-met(1,1)*met(2,0))*stilde(1)+ \
741  (met(0,1)*met(2,0)-met(0,0)*met(2,1))*stilde(2)+ \
742  (met(0,0)*met(1,1)-met(0,1)*met(1,0))*stilde(3)
743  end do
744  endLoopsMask3d()
745  else ! free stream preserving method
746 cccccccccccc
747  ! *********** DISPLACEMENT : Curvilinear Grid (free stream preserving) ****************
748 
749  ! Use momentum equations to get (rx,ry,rz)(0).(s11,s21,s31)(-1) = s11tilde
750  ! (1) (rx,ry,rz).D_r1[(s11,s21,s31)] + (sx,sy,sz).D_r2[(s11,s21,s31)] + (tx,ty,tz).D_r3[(s11,s21,s31)] = rho * u_tt
751  ! s11tilde s21tilde s31tilde
752  ! (2) Use extrapolated values to get (sx,sy,sz)(0).(s11,s21,s31)(-1) = s21tilde
753  ! (3) Use extrapolated values to get (tx,ty,tz)(0).(s11,s21,s31)(-1) = s31tilde
754  ! To give 3 equations for (s11,s21,s31) on the ghost point:
755  ! (rx)(0) s11(-1) + (ry)(0) s21(-1) + (rz)(0) s31(-1) = f1 = s11tilde (from momentum eqn)
756  ! (sx)(0) s11(-1) + (sy)(0) s21(-1) + (sz)(0) s31(-1) = f2 = s21tilde (from extrapolated values)
757  ! (tx)(0) s11(-1) + (ty)(0) s21(-1) + (tz)(0) s31(-1) = f3 = s31tilde (from extrapolated values)
758  ! Solve: (note that det is the inverse of the determinant det[rx,ry,rz; sx,sy,sz; tx,ty,tz])
759  ! s11(-1) = ((sy*tz-sz*ty)*f1 + (rz*ty-ry*tz)*f2 + (ry*sz-rz*sy)*f3)*det
760  ! s21(-1) = ((tx*sz-sx*tz)*f1 + (rx*tz-rz*tx)*f2 + (rz*sx-rx*sz)*f3)*det
761  ! s31(-1) = ((sx*ty-sy*tx)*f1 + (ry*tx-rx*ty)*f2 + (rx*sy-ry*sx)*f3)*det
762  !
763  ! (A similar expression holds for other stresses)
764  beginLoopsMask3d()
765  accel(1) = rho*bcf(side,axis,i1,i2,i3,s11c)
766  accel(2) = rho*bcf(side,axis,i1,i2,i3,s12c)
767  accel(3) = rho*bcf(side,axis,i1,i2,i3,s13c)
768 
769  met(0,0) = rx(i1,i2,i3,0,0)
770  met(0,1) = rx(i1,i2,i3,0,1)
771  met(0,2) = rx(i1,i2,i3,0,2)
772  met(1,0) = rx(i1,i2,i3,1,0)
773  met(1,1) = rx(i1,i2,i3,1,1)
774  met(1,2) = rx(i1,i2,i3,1,2)
775  met(2,0) = rx(i1,i2,i3,2,0)
776  met(2,1) = rx(i1,i2,i3,2,1)
777  met(2,2) = rx(i1,i2,i3,2,2)
778 
779  ! loop over stress components
780  do isc = 1,3
781  do idot = 1,3
782  ! these are extrapolated components
783  stilde(idot) = rx(i1,i2,i3,idot-1,0)*u(i1-is1,i2-is2,i3-is3,sc(1,isc))+ \
784  rx(i1,i2,i3,idot-1,1)*u(i1-is1,i2-is2,i3-is3,sc(2,isc))+ \
785  rx(i1,i2,i3,idot-1,2)*u(i1-is1,i2-is2,i3-is3,sc(3,isc))
786  end do
787  ! now override in the direction we are currently looking
788  stilde(axis+1) = (rx(i1,i2,i3,axis,0)*u(i1+is1,i2+is2,i3+is3,sc(1,isc))+ \
789  rx(i1,i2,i3,axis,1)*u(i1+is1,i2+is2,i3+is3,sc(2,isc))+ \
790  rx(i1,i2,i3,axis,2)*u(i1+is1,i2+is2,i3+is3,sc(3,isc)))- \
791  2.0*dr(axis)*is*(accel(isc)- \
792  (1.0-delta(axis+1,1))* \
793  (rx(i1,i2,i3,0,0)*(u(i1+1,i2,i3,sc(1,isc))-u(i1-1,i2,i3,sc(1,isc)))+ \
794  rx(i1,i2,i3,0,1)*(u(i1+1,i2,i3,sc(2,isc))-u(i1-1,i2,i3,sc(2,isc)))+ \
795  rx(i1,i2,i3,0,2)*(u(i1+1,i2,i3,sc(3,isc))-u(i1-1,i2,i3,sc(3,isc))))/(2.0*dr(0))- \
796  (1.0-delta(axis+1,2))* \
797  (rx(i1,i2,i3,1,0)*(u(i1,i2+1,i3,sc(1,isc))-u(i1,i2-1,i3,sc(1,isc)))+ \
798  rx(i1,i2,i3,1,1)*(u(i1,i2+1,i3,sc(2,isc))-u(i1,i2-1,i3,sc(2,isc)))+ \
799  rx(i1,i2,i3,1,2)*(u(i1,i2+1,i3,sc(3,isc))-u(i1,i2-1,i3,sc(3,isc))))/(2.0*dr(1))- \
800  (1.0-delta(axis+1,3))* \
801  (rx(i1,i2,i3,2,0)*(u(i1,i2,i3+1,sc(1,isc))-u(i1,i2,i3-1,sc(1,isc)))+ \
802  rx(i1,i2,i3,2,1)*(u(i1,i2,i3+1,sc(2,isc))-u(i1,i2,i3-1,sc(2,isc)))+ \
803  rx(i1,i2,i3,2,2)*(u(i1,i2,i3+1,sc(3,isc))-u(i1,i2,i3-1,sc(3,isc))))/(2.0*dr(2)))
804 
805  u(i1-is1,i2-is2,i3-is3,sc(1,isc)) = ((met(1,1)*met(2,2)-met(1,2)*met(2,1))*stilde(1)+ \
806  (met(0,2)*met(2,1)-met(0,1)*met(2,2))*stilde(2)+ \
807  (met(0,1)*met(1,2)-met(0,2)*met(1,1))*stilde(3))*det(i1,i2,i3)
808  u(i1-is1,i2-is2,i3-is3,sc(2,isc)) = ((met(1,2)*met(2,0)-met(1,0)*met(2,2))*stilde(1)+ \
809  (met(0,0)*met(2,2)-met(0,2)*met(2,0))*stilde(2)+ \
810  (met(0,2)*met(1,0)-met(0,0)*met(1,2))*stilde(3))*det(i1,i2,i3)
811  u(i1-is1,i2-is2,i3-is3,sc(3,isc)) = ((met(1,0)*met(2,1)-met(1,1)*met(2,0))*stilde(1)+ \
812  (met(0,1)*met(2,0)-met(0,0)*met(2,1))*stilde(2)+ \
813  (met(0,0)*met(1,1)-met(0,1)*met(1,0))*stilde(3))*det(i1,i2,i3)
814  end do
815  endLoopsMask3d()
816 cccccccccccc
817  end if
818  end if ! end gridType
819 
820  else if( boundaryCondition(side,axis).eq.tractionBC ) then
821  ! **************** TRACTION : Neumann type conditions ******************
822  if( gridType.eq.rectangular )then
823 
824  ! ********* TRACTION : Cartesian Grid **********
825 
826  ! Assign displacements on the ghost points from given tractions on the boundary
827  ! s11 = kappa*u.x + lambda*( v.y + w.z )
828  ! s22 = kappa*v.y + lambda*( u.x + w.z )
829  ! s33 = kappa*w.z + lambda*( u.x + v.y )
830  ! s12 = s21 = mu*( u.y + v.x )
831  ! s13 = s31 = mu*( w.x + u.z )
832  !
833  ! an1*s11 + an2*s21 + an3*s31 = f1
834  ! an1*s12 + an2*s22 + an3*s32 = f2
835  ! an1*s13 + an2*s23 + an3*s33 = f3
836 
837  ! Assign velocities on the ghost points from given time derivatives of the tractions on the boundary
838  beginLoopsMask3d()
839  f1 = bcf(side,axis,i1,i2,i3,uc)
840  f2 = bcf(side,axis,i1,i2,i3,vc)
841  f3 = bcf(side,axis,i1,i2,i3,wc)
842 
843  if( axis.eq.0 )then
844  an1 = -is
845  an2 = 0.0
846  an3 = 0.0
847  else if( axis.eq.1 ) then
848  an1 = 0.0
849  an2 = -is
850  an3 = 0.0
851  else
852  an1 = 0.0
853  an2 = 0.0
854  an3 = -is
855  end if
856 
857  dux = diffr1(uc,dx(0))*(1.0-delta(axis+1,1))
858  duy = diffr2(uc,dx(1))*(1.0-delta(axis+1,2))
859  duz = diffr3(uc,dx(2))*(1.0-delta(axis+1,3))
860 
861  dvx = diffr1(vc,dx(0))*(1.0-delta(axis+1,1))
862  dvy = diffr2(vc,dx(1))*(1.0-delta(axis+1,2))
863  dvz = diffr3(vc,dx(2))*(1.0-delta(axis+1,3))
864 
865  dwx = diffr1(wc,dx(0))*(1.0-delta(axis+1,1))
866  dwy = diffr2(wc,dx(1))*(1.0-delta(axis+1,2))
867  dwz = diffr3(wc,dx(2))*(1.0-delta(axis+1,3))
868 
869  f1 = f1-an1*(kappa*dux+lambda*(dvy+dwz))- \
870  an2*(mu*(dvx+duy))- \
871  an3*(mu*(dwx+duz))
872  f2 = f2-an1*(mu*(dvx+duy))- \
873  an2*(kappa*dvy+lambda*(dux+dwz))- \
874  an3*(mu*(dwy+dvz))
875  f3 = f3-an1*(mu*(dwx+duz))- \
876  an2*(mu*(dwy+dvz))- \
877  an3*(kappa*dwz+lambda*(dux+dvy))
878 
879  ! in the Cartesian case all that survives in the Matrix are the diagonal terms
880  f1 = f1/(an1*kappa +an2*mu +an3*mu)
881  f2 = f2/(an1*mu +an2*kappa +an3*mu)
882  f3 = f3/(an1*mu +an2*mu +an3*kappa)
883 
884  u(i1-is1,i2-is2,i3-is3,uc) = -is*2.0*dx(axis)*f1+u(i1+is1,i2+is2,i3+is3,uc)
885  u(i1-is1,i2-is2,i3-is3,vc) = -is*2.0*dx(axis)*f2+u(i1+is1,i2+is2,i3+is3,vc)
886  u(i1-is1,i2-is2,i3-is3,wc) = -is*2.0*dx(axis)*f3+u(i1+is1,i2+is2,i3+is3,wc)
887 
888  !!! now do velocities
889  fdot1 = bcf(side,axis,i1,i2,i3,v1c)
890  fdot2 = bcf(side,axis,i1,i2,i3,v2c)
891  fdot3 = bcf(side,axis,i1,i2,i3,v3c)
892 
893  dux = diffr1(v1c,dx(0))*(1.0-delta(axis+1,1))
894  duy = diffr2(v1c,dx(1))*(1.0-delta(axis+1,2))
895  duz = diffr3(v1c,dx(2))*(1.0-delta(axis+1,3))
896 
897  dvx = diffr1(v2c,dx(0))*(1.0-delta(axis+1,1))
898  dvy = diffr2(v2c,dx(1))*(1.0-delta(axis+1,2))
899  dvz = diffr3(v2c,dx(2))*(1.0-delta(axis+1,3))
900 
901  dwx = diffr1(v3c,dx(0))*(1.0-delta(axis+1,1))
902  dwy = diffr2(v3c,dx(1))*(1.0-delta(axis+1,2))
903  dwz = diffr3(v3c,dx(2))*(1.0-delta(axis+1,3))
904 
905  fdot1 = fdot1-an1*(kappa*dux+lambda*(dvy+dwz))- \
906  an2*(mu*(dvx+duy))- \
907  an3*(mu*(dwx+duz))
908  fdot2 = fdot2-an1*(mu*(dvx+duy))- \
909  an2*(kappa*dvy+lambda*(dux+dwz))- \
910  an3*(mu*(dwy+dvz))
911  fdot3 = fdot3-an1*(mu*(dwx+duz))- \
912  an2*(mu*(dwy+dvz))- \
913  an3*(kappa*dwz+lambda*(dux+dvy))
914 
915  ! in the Cartesian case all that survives in the Matrix are the diagonal terms
916  fdot1 = fdot1/(an1*kappa +an2*mu +an3*mu)
917  fdot2 = fdot2/(an1*mu +an2*kappa +an3*mu)
918  fdot3 = fdot3/(an1*mu +an2*mu +an3*kappa)
919 
920  u(i1-is1,i2-is2,i3-is3,v1c) = -is*2.0*dx(axis)*fdot1+u(i1+is1,i2+is2,i3+is3,v1c)
921  u(i1-is1,i2-is2,i3-is3,v2c) = -is*2.0*dx(axis)*fdot2+u(i1+is1,i2+is2,i3+is3,v2c)
922  u(i1-is1,i2-is2,i3-is3,v3c) = -is*2.0*dx(axis)*fdot3+u(i1+is1,i2+is2,i3+is3,v3c)
923 
924  endLoopsMask3d()
925  else
926  ! *********** TRACTION : Curvilinear Grid ****************
927  beginLoopsMask3d()
928  f1 = bcf(side,axis,i1,i2,i3,uc)
929  f2 = bcf(side,axis,i1,i2,i3,vc)
930  f3 = bcf(side,axis,i1,i2,i3,wc)
931 
932  fdot1 = bcf(side,axis,i1,i2,i3,v1c)
933  fdot2 = bcf(side,axis,i1,i2,i3,v2c)
934  fdot3 = bcf(side,axis,i1,i2,i3,v3c)
935 
936  met(0,0) = rx(i1,i2,i3,0,0)
937  met(0,1) = rx(i1,i2,i3,0,1)
938  met(0,2) = rx(i1,i2,i3,0,2)
939  met(1,0) = rx(i1,i2,i3,1,0)
940  met(1,1) = rx(i1,i2,i3,1,1)
941  met(1,2) = rx(i1,i2,i3,1,2)
942  met(2,0) = rx(i1,i2,i3,2,0)
943  met(2,1) = rx(i1,i2,i3,2,1)
944  met(2,2) = rx(i1,i2,i3,2,2)
945 
946  ! (an1,an2,an3) = outward normal
947  aNormi = 1.0/max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
948  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
949  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
950  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
951 
952  dur(0) = diffr1(uc,dr(0))*(1.0-delta(axis+1,1))
953  dur(1) = diffr2(uc,dr(1))*(1.0-delta(axis+1,2))
954  dur(2) = diffr3(uc,dr(2))*(1.0-delta(axis+1,3))
955 
956  dvr(0) = diffr1(vc,dr(0))*(1.0-delta(axis+1,1))
957  dvr(1) = diffr2(vc,dr(1))*(1.0-delta(axis+1,2))
958  dvr(2) = diffr3(vc,dr(2))*(1.0-delta(axis+1,3))
959 
960  dwr(0) = diffr1(wc,dr(0))*(1.0-delta(axis+1,1))
961  dwr(1) = diffr2(wc,dr(1))*(1.0-delta(axis+1,2))
962  dwr(2) = diffr3(wc,dr(2))*(1.0-delta(axis+1,3))
963 
964  dv1r(0) = diffr1(v1c,dr(0))*(1.0-delta(axis+1,1))
965  dv1r(1) = diffr2(v1c,dr(1))*(1.0-delta(axis+1,2))
966  dv1r(2) = diffr3(v1c,dr(2))*(1.0-delta(axis+1,3))
967 
968  dv2r(0) = diffr1(v2c,dr(0))*(1.0-delta(axis+1,1))
969  dv2r(1) = diffr2(v2c,dr(1))*(1.0-delta(axis+1,2))
970  dv2r(2) = diffr3(v2c,dr(2))*(1.0-delta(axis+1,3))
971 
972  dv3r(0) = diffr1(v3c,dr(0))*(1.0-delta(axis+1,1))
973  dv3r(1) = diffr2(v3c,dr(1))*(1.0-delta(axis+1,2))
974  dv3r(2) = diffr3(v3c,dr(2))*(1.0-delta(axis+1,3))
975 
976  do isc = 0,2
977  mat(0,0) = an1*kappa*met(isc,0) +an2*mu*met(isc,1) +an3*mu*met(isc,2)
978  mat(0,1) = an1*lambda*met(isc,1)+an2*mu*met(isc,0)
979  mat(0,2) = an1*lambda*met(isc,2) +an3*mu*met(isc,0)
980  mat(1,0) = an1*mu*met(isc,1) +an2*lambda*met(isc,0)
981  mat(1,1) = an1*mu*met(isc,0) +an2*kappa*met(isc,1) +an3*mu*met(isc,2)
982  mat(1,2) = an2*lambda*met(isc,2)+an3*mu*met(isc,1)
983  mat(2,0) = an1*mu*met(isc,2) +an3*lambda*met(isc,0)
984  mat(2,1) = an2*mu*met(isc,2) +an3*lambda*met(isc,1)
985  mat(2,2) = an1*mu*met(isc,0) +an2*mu*met(isc,1) +an3*kappa*met(isc,2)
986 
987  f1 = f1-(mat(0,0)*dur(isc)+mat(0,1)*dvr(isc)+mat(0,2)*dwr(isc))
988  f2 = f2-(mat(1,0)*dur(isc)+mat(1,1)*dvr(isc)+mat(1,2)*dwr(isc))
989  f3 = f3-(mat(2,0)*dur(isc)+mat(2,1)*dvr(isc)+mat(2,2)*dwr(isc))
990 
991  fdot1 = fdot1-(mat(0,0)*dv1r(isc)+mat(0,1)*dv2r(isc)+mat(0,2)*dv3r(isc))
992  fdot2 = fdot2-(mat(1,0)*dv1r(isc)+mat(1,1)*dv2r(isc)+mat(1,2)*dv3r(isc))
993  fdot3 = fdot3-(mat(2,0)*dv1r(isc)+mat(2,1)*dv2r(isc)+mat(2,2)*dv3r(isc))
994 
995  if( axis.eq.isc ) then
996  lhs(0,0) = mat(0,0)
997  lhs(0,1) = mat(0,1)
998  lhs(0,2) = mat(0,2)
999  lhs(1,0) = mat(1,0)
1000  lhs(1,1) = mat(1,1)
1001  lhs(1,2) = mat(1,2)
1002  lhs(2,0) = mat(2,0)
1003  lhs(2,1) = mat(2,1)
1004  lhs(2,2) = mat(2,2)
1005  end if
1006  end do
1007 
1008 
1009  !! solve linear systems to get the solution (grid derivatives)
1010  rhs(0,0) = f1
1011  rhs(1,0) = f2
1012  rhs(2,0) = f3
1013  rhs(0,1) = fdot1
1014  rhs(1,1) = fdot2
1015  rhs(2,1) = fdot3
1016 
1017  call dgesv( 3,2,lhs,3,ipiv,rhs,3,info )
1018  if( info.ne.0 ) then
1019  write(6,*)'Error (compat3D) : error in linear system'
1020  stop
1021  end if
1022 
1023  u(i1-is1,i2-is2,i3-is3,uc) = -is*2.0*rhs(0,0)*dr(axis)+u(i1+is1,i2+is2,i3+is3,uc)
1024  u(i1-is1,i2-is2,i3-is3,vc) = -is*2.0*rhs(1,0)*dr(axis)+u(i1+is1,i2+is2,i3+is3,vc)
1025  u(i1-is1,i2-is2,i3-is3,wc) = -is*2.0*rhs(2,0)*dr(axis)+u(i1+is1,i2+is2,i3+is3,wc)
1026 
1027  u(i1-is1,i2-is2,i3-is3,v1c) = -is*2.0*rhs(0,1)*dr(axis)+u(i1+is1,i2+is2,i3+is3,v1c)
1028  u(i1-is1,i2-is2,i3-is3,v2c) = -is*2.0*rhs(1,1)*dr(axis)+u(i1+is1,i2+is2,i3+is3,v2c)
1029  u(i1-is1,i2-is2,i3-is3,v3c) = -is*2.0*rhs(2,1)*dr(axis)+u(i1+is1,i2+is2,i3+is3,v3c)
1030 
1031  endLoopsMask3d()
1032  end if ! gridType
1033  else if( boundaryCondition(side,axis).eq.slipWall ) then
1034  ! **************** SLIPWALL : Neumann type conditions ******************
1035  if( gridType.eq.rectangular ) then
1036  ! ********* SLIPWALL : Cartesian Grid **********
1037  else
1038  ! ********* SLIPWALL : Curvilinear Grid **********
1039  end if ! gridType
1040 
1041  else if( boundaryCondition(side,axis).gt.0 .and. boundaryCondition(side,axis).ne.dirichletBoundaryCondition ) then
1042  write(*,'("smg3d:BC: unknown BC: side,axis,grid, boundaryCondition=",i2,i2,i4,i8)') side,axis,grid,boundaryCondition(side,axis)
1043 
1044  end if ! bc
1045  endLoopOverSides()
1046 
1047 c*******
1048 c******* Secondary Dirichlet conditions for the tangential components of stress (tractionBC only) ********
1049 c*******
1050 c if( .false. ) then
1051  beginLoopOverSides(numGhost,numGhost)
1052  if( boundaryCondition(side,axis).eq.tractionBC )then
1053  if( gridType.eq.rectangular )then
1054  if( axis.eq.0 )then
1055  beginLoopsMask3d()
1056  u1x = diffr1(uc,dx(0))
1057  u1y = diffr2(uc,dx(1))
1058  u1z = diffr3(uc,dx(2))
1059 
1060  u2x = diffr1(vc,dx(0))
1061  u2y = diffr2(vc,dx(1))
1062  u2z = diffr3(vc,dx(2))
1063 
1064  u3x = diffr1(wc,dx(0))
1065  u3y = diffr2(wc,dx(1))
1066  u3z = diffr3(wc,dx(2))
1067 
1068  u(i1,i2,i3,s21c) = mu*(u1y+u2x)
1069  u(i1,i2,i3,s22c) = kappa*u2y+lambda*(u1x+u3z)
1070  u(i1,i2,i3,s23c) = mu*(u3y+u2z)
1071 
1072  u(i1,i2,i3,s31c) = mu*(u3x+u1z)
1073  u(i1,i2,i3,s32c) = mu*(u3y+u2z)
1074  u(i1,i2,i3,s33c) = kappa*u3z+lambda*(u1x+u2y)
1075  endLoopsMask3d()
1076  else if( axis.eq.1 ) then
1077  beginLoopsMask3d()
1078  u1x = diffr1(uc,dx(0))
1079  u1y = diffr2(uc,dx(1))
1080  u1z = diffr3(uc,dx(2))
1081 
1082  u2x = diffr1(vc,dx(0))
1083  u2y = diffr2(vc,dx(1))
1084  u2z = diffr3(vc,dx(2))
1085 
1086  u3x = diffr1(wc,dx(0))
1087  u3y = diffr2(wc,dx(1))
1088  u3z = diffr3(wc,dx(2))
1089 
1090  u(i1,i2,i3,s11c) = kappa*u1x+lambda*(u2y+u3z)
1091  u(i1,i2,i3,s12c) = mu*(u2x+u1y)
1092  u(i1,i2,i3,s13c) = mu*(u3x+u1z)
1093 
1094  u(i1,i2,i3,s31c) = mu*(u3x+u1z)
1095  u(i1,i2,i3,s32c) = mu*(u3y+u2z)
1096  u(i1,i2,i3,s33c) = kappa*u3z+lambda*(u1x+u2y)
1097  endLoopsMask3d()
1098  else
1099  beginLoopsMask3d()
1100  u1x = diffr1(uc,dx(0))
1101  u1y = diffr2(uc,dx(1))
1102  u1z = diffr3(uc,dx(2))
1103 
1104  u2x = diffr1(vc,dx(0))
1105  u2y = diffr2(vc,dx(1))
1106  u2z = diffr3(vc,dx(2))
1107 
1108  u3x = diffr1(wc,dx(0))
1109  u3y = diffr2(wc,dx(1))
1110  u3z = diffr3(wc,dx(2))
1111 
1112  u(i1,i2,i3,s11c) = kappa*u1x+lambda*(u2y+u3z)
1113  u(i1,i2,i3,s12c) = mu*(u2x+u1y)
1114  u(i1,i2,i3,s13c) = mu*(u3x+u1z)
1115 
1116  u(i1,i2,i3,s21c) = mu*(u1y+u2x)
1117  u(i1,i2,i3,s22c) = kappa*u2y+lambda*(u1x+u3z)
1118  u(i1,i2,i3,s23c) = mu*(u3y+u2z)
1119  endLoopsMask3d()
1120  end if
1121 
1122  else ! curvilinear
1123  if( axis.eq.0 ) then
1124  tan1c = 1
1125  tan2c = 2
1126  else if( axis.eq.1 ) then
1127  tan1c = 2
1128  tan2c = 0
1129  else
1130  tan1c = 0
1131  tan2c = 1
1132  end if
1133  beginLoopsMask3d()
1134  u1r = diffr1(uc,dr(0))
1135  u1s = diffr2(uc,dr(1))
1136  u1t = diffr3(uc,dr(2))
1137 
1138  u2r = diffr1(vc,dr(0))
1139  u2s = diffr2(vc,dr(1))
1140  u2t = diffr3(vc,dr(2))
1141 
1142  u3r = diffr1(wc,dr(0))
1143  u3s = diffr2(wc,dr(1))
1144  u3t = diffr3(wc,dr(2))
1145 
1146  u1x = u1r*rx(i1,i2,i3,0,0)+u1s*rx(i1,i2,i3,1,0)+u1t*rx(i1,i2,i3,2,0)
1147  u1y = u1r*rx(i1,i2,i3,0,1)+u1s*rx(i1,i2,i3,1,1)+u1t*rx(i1,i2,i3,2,1)
1148  u1z = u1r*rx(i1,i2,i3,0,2)+u1s*rx(i1,i2,i3,1,2)+u1t*rx(i1,i2,i3,2,2)
1149 
1150  u2x = u2r*rx(i1,i2,i3,0,0)+u2s*rx(i1,i2,i3,1,0)+u2t*rx(i1,i2,i3,2,0)
1151  u2y = u2r*rx(i1,i2,i3,0,1)+u2s*rx(i1,i2,i3,1,1)+u2t*rx(i1,i2,i3,2,1)
1152  u2z = u2r*rx(i1,i2,i3,0,2)+u2s*rx(i1,i2,i3,1,2)+u2t*rx(i1,i2,i3,2,2)
1153 
1154  u3x = u3r*rx(i1,i2,i3,0,0)+u3s*rx(i1,i2,i3,1,0)+u3t*rx(i1,i2,i3,2,0)
1155  u3y = u3r*rx(i1,i2,i3,0,1)+u3s*rx(i1,i2,i3,1,1)+u3t*rx(i1,i2,i3,2,1)
1156  u3z = u3r*rx(i1,i2,i3,0,2)+u3s*rx(i1,i2,i3,1,2)+u3t*rx(i1,i2,i3,2,2)
1157 
1158  s11t = kappa*u1x+lambda*(u2y+u3z)
1159  s12t = mu*(u2x+u1y)
1160  s13t = mu*(u3x+u1z)
1161 
1162  s21t = mu*(u2x+u1y)
1163  s22t = kappa*u2y+lambda*(u1x+u3z)
1164  s23t = mu*(u3y+u2z)
1165 
1166  s31t = mu*(u3x+u1z)
1167  s32t = mu*(u3y+u2z)
1168  s33t = kappa*u3z+lambda*(u1x+u2y)
1169 
1170  ! (an1,an2,3) = outward unit normal
1171  aNormi = 1.0/max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
1172  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
1173  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
1174  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
1175 
1176  !! do the signs of sn and tn matter?? ... I do not think so
1177  sn1 = rx(i1,i2,i3,tan1c,0)
1178  sn2 = rx(i1,i2,i3,tan1c,1)
1179  sn3 = rx(i1,i2,i3,tan1c,2)
1180 
1181  tn1 = rx(i1,i2,i3,tan2c,0)
1182  tn2 = rx(i1,i2,i3,tan2c,1)
1183  tn3 = rx(i1,i2,i3,tan2c,2)
1184 
1185  ! set sn to be part of sn which is orthogonal to an
1186  alpha = an1*sn1+an2*sn2+an3*sn3
1187  sn1 = sn1-alpha*an1
1188  sn2 = sn2-alpha*an2
1189  sn3 = sn3-alpha*an3
1190  ! normalize sn
1191  aNormi = 1.0/max(epsx,sqrt(sn1**2+sn2**2+sn3**2))
1192  sn1 = sn1*aNormi
1193  sn2 = sn2*aNormi
1194  sn3 = sn3*aNormi
1195 
1196  ! set tn to be part of tn which is orthogonal to an and sn
1197  alpha = an1*tn1+an2*tn2+an3*tn3
1198  tn1 = tn1-alpha*an1
1199  tn2 = tn2-alpha*an2
1200  tn3 = tn3-alpha*an3
1201  alpha = sn1*tn1+sn2*tn2+sn3*tn3
1202  tn1 = tn1-alpha*sn1
1203  tn2 = tn2-alpha*sn2
1204  tn3 = tn3-alpha*sn3
1205  ! normalize tn
1206  aNormi = 1.0/max(epsx,sqrt(tn1**2+tn2**2+tn3**2))
1207  tn1 = tn1*aNormi
1208  tn2 = tn2*aNormi
1209  tn3 = tn3*aNormi
1210 
1211  ! compute components of stress in normal direction (primary condition)
1212  ns1 = an1*u(i1,i2,i3,s11c)+an2*u(i1,i2,i3,s21c)+an3*u(i1,i2,i3,s31c)
1213  ns2 = an1*u(i1,i2,i3,s12c)+an2*u(i1,i2,i3,s22c)+an3*u(i1,i2,i3,s32c)
1214  ns3 = an1*u(i1,i2,i3,s13c)+an2*u(i1,i2,i3,s23c)+an3*u(i1,i2,i3,s33c)
1215 
1216  ! compute componenets of stress in 1st tangential direction (secondary condition)
1217  ss1 = sn1*s11t+sn2*s21t+sn3*s31t
1218  ss2 = sn1*s12t+sn2*s22t+sn3*s32t
1219  ss3 = sn1*s13t+sn2*s23t+sn3*s33t
1220 
1221  ! compute componenets of stress in 2nd tangential direction (secondary condition)
1222  ts1 = tn1*s11t+tn2*s21t+tn3*s31t
1223  ts2 = tn1*s12t+tn2*s22t+tn3*s32t
1224  ts3 = tn1*s13t+tn2*s23t+tn3*s33t
1225 
1226  u(i1,i2,i3,s11c) = an1*ns1+sn1*ss1+tn1*ts1
1227  u(i1,i2,i3,s12c) = an1*ns2+sn1*ss2+tn1*ts2
1228  u(i1,i2,i3,s13c) = an1*ns3+sn1*ss3+tn1*ts3
1229 
1230  u(i1,i2,i3,s21c) = an2*ns1+sn2*ss1+tn2*ts1
1231  u(i1,i2,i3,s22c) = an2*ns2+sn2*ss2+tn2*ts2
1232  u(i1,i2,i3,s23c) = an2*ns3+sn2*ss3+tn2*ts3
1233 
1234  u(i1,i2,i3,s31c) = an3*ns1+sn3*ss1+tn3*ts1
1235  u(i1,i2,i3,s32c) = an3*ns2+sn3*ss2+tn3*ts2
1236  u(i1,i2,i3,s33c) = an3*ns3+sn3*ss3+tn3*ts3
1237  endLoopsMask3d()
1238 
1239  end if ! end gridType
1240 
1241  end if ! bc
1242  endLoopOverSides()
1243 
1244 c set tangential components of stress on the boundary (TZ forcing, if necessary)
1245  if( twilightZone.ne.0 ) then
1246 c if( .false. ) then
1247 
1248  beginLoopOverSides(numGhost,numGhost)
1249  if( boundaryCondition(side,axis).eq.tractionBC )then
1250 
1251  if( gridType.eq.rectangular )then
1252  if( axis.eq.0 )then
1253  beginLoopsMask3d()
1254  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
1255  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
1256  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
1257 
1258  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
1259  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
1260  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
1261 
1262  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
1263  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
1264  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
1265 
1266  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,s21e )
1267  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,s22e )
1268  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,s23e )
1269 
1270  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,s31e )
1271  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,s32e )
1272  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,s33e )
1273 
1274  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)+s21e-mu*(u1y+u2x)
1275  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)+s22e-(kappa*u2y+lambda*(u1x+u3z))
1276  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)+s23e-mu*(u3y+u2z)
1277 
1278  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)+s31e-mu*(u3x+u1z)
1279  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)+s32e-mu*(u3y+u2z)
1280  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)+s33e-(kappa*u3z+lambda*(u1x+u2y))
1281  endLoopsMask3d()
1282  else if( axis.eq.1 ) then
1283  beginLoopsMask3d()
1284  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
1285  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
1286  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
1287 
1288  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
1289  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
1290  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
1291 
1292  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
1293  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
1294  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
1295 
1296  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,s11e )
1297  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,s12e )
1298  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,s13e )
1299 
1300  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,s31e )
1301  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,s32e )
1302  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,s33e )
1303 
1304  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)+s11e-(kappa*u1x+lambda*(u2y+u3z))
1305  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)+s12e-mu*(u2x+u1y)
1306  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)+s13e-mu*(u3x+u1z)
1307 
1308  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)+s31e-mu*(u3x+u1z)
1309  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)+s32e-mu*(u3y+u2z)
1310  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)+s33e-(kappa*u3z+lambda*(u1x+u2y))
1311  endLoopsMask3d()
1312  else
1313  beginLoopsMask3d()
1314  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
1315  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
1316  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
1317 
1318  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
1319  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
1320  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
1321 
1322  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
1323  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
1324  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
1325 
1326  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,s11e )
1327  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,s12e )
1328  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,s13e )
1329 
1330  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,s21e )
1331  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,s22e )
1332  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,s23e )
1333 
1334  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)+s11e-(kappa*u1x+lambda*(u2y+u3z))
1335  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)+s12e-mu*(u2x+u1y)
1336  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)+s13e-mu*(u3x+u1z)
1337 
1338  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)+s21e-mu*(u1y+u2x)
1339  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)+s22e-(kappa*u2y+lambda*(u1x+u3z))
1340  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)+s23e-mu*(u3y+u2z)
1341  endLoopsMask3d()
1342  end if
1343 
1344  else ! curvilinear
1345  if( axis.eq.0 ) then
1346  tan1c = 1
1347  tan2c = 2
1348  else if( axis.eq.1 ) then
1349  tan1c = 2
1350  tan2c = 0
1351  else
1352  tan1c = 0
1353  tan2c = 1
1354  end if
1355  beginLoopsMask3d()
1356  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
1357  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
1358  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
1359 
1360  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
1361  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
1362  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
1363 
1364  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
1365  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
1366  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
1367 
1368  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,s11e )
1369  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,s12e )
1370  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,s13e )
1371 
1372  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,s21e )
1373  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,s22e )
1374  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,s23e )
1375 
1376  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,s31e )
1377  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,s32e )
1378  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,s33e )
1379 
1380  s11t = kappa*u1x+lambda*(u2y+u3z)
1381  s12t = mu*(u2x+u1y)
1382  s13t = mu*(u3x+u1z)
1383 
1384  s21t = mu*(u2x+u1y)
1385  s22t = kappa*u2y+lambda*(u1x+u3z)
1386  s23t = mu*(u3y+u2z)
1387 
1388  s31t = mu*(u3x+u1z)
1389  s32t = mu*(u3y+u2z)
1390  s33t = kappa*u3z+lambda*(u1x+u2y)
1391 
1392  ! (an1,an2,3) = outward unit normal
1393  aNormi = 1./max(epsx,sqrt(rx(i1,i2,i3,axis,0)**2+rx(i1,i2,i3,axis,1)**2+rx(i1,i2,i3,axis,2)**2))
1394  an1 = -is*rx(i1,i2,i3,axis,0)*aNormi
1395  an2 = -is*rx(i1,i2,i3,axis,1)*aNormi
1396  an3 = -is*rx(i1,i2,i3,axis,2)*aNormi
1397 
1398  sn1 = rx(i1,i2,i3,tan1c,0)
1399  sn2 = rx(i1,i2,i3,tan1c,1)
1400  sn3 = rx(i1,i2,i3,tan1c,2)
1401 
1402  tn1 = rx(i1,i2,i3,tan2c,0)
1403  tn2 = rx(i1,i2,i3,tan2c,1)
1404  tn3 = rx(i1,i2,i3,tan2c,2)
1405 
1406  ! set sn to be part of sn which is orthogonal to an
1407  alpha = an1*sn1+an2*sn2+an3*sn3
1408  sn1 = sn1-alpha*an1
1409  sn2 = sn2-alpha*an2
1410  sn3 = sn3-alpha*an3
1411  ! normalize sn
1412  aNormi = 1./max(epsx,sqrt(sn1**2+sn2**2+sn3**2))
1413  sn1 = sn1*aNormi
1414  sn2 = sn2*aNormi
1415  sn3 = sn3*aNormi
1416 
1417  ! set tn to be part of tn which is orthogonal to an and sn
1418  alpha = an1*tn1+an2*tn2+an3*tn3
1419  tn1 = tn1-alpha*an1
1420  tn2 = tn2-alpha*an2
1421  tn3 = tn3-alpha*an3
1422  alpha = sn1*tn1+sn2*tn2+sn3*tn3
1423  tn1 = tn1-alpha*sn1
1424  tn2 = tn2-alpha*sn2
1425  tn3 = tn3-alpha*sn3
1426  ! normalize tn
1427  aNormi = 1./max(epsx,sqrt(tn1**2+tn2**2+tn3**2))
1428  tn1 = tn1*aNormi
1429  tn2 = tn2*aNormi
1430  tn3 = tn3*aNormi
1431 
1432  ! compute components of stress in normal direction (leave these alone)
1433  ns1 = an1*u(i1,i2,i3,s11c)+an2*u(i1,i2,i3,s21c)+an3*u(i1,i2,i3,s31c)
1434  ns2 = an1*u(i1,i2,i3,s12c)+an2*u(i1,i2,i3,s22c)+an3*u(i1,i2,i3,s32c)
1435  ns3 = an1*u(i1,i2,i3,s13c)+an2*u(i1,i2,i3,s23c)+an3*u(i1,i2,i3,s33c)
1436 
1437  ! compute componenets of stress in tangential directions (add forcing to these)
1438  ss1 = sn1*u(i1,i2,i3,s11c)+sn2*u(i1,i2,i3,s21c)+sn3*u(i1,i2,i3,s31c)
1439  ss2 = sn1*u(i1,i2,i3,s12c)+sn2*u(i1,i2,i3,s22c)+sn3*u(i1,i2,i3,s32c)
1440  ss3 = sn1*u(i1,i2,i3,s13c)+sn2*u(i1,i2,i3,s23c)+sn3*u(i1,i2,i3,s33c)
1441  ts1 = tn1*u(i1,i2,i3,s11c)+tn2*u(i1,i2,i3,s21c)+tn3*u(i1,i2,i3,s31c)
1442  ts2 = tn1*u(i1,i2,i3,s12c)+tn2*u(i1,i2,i3,s22c)+tn3*u(i1,i2,i3,s32c)
1443  ts3 = tn1*u(i1,i2,i3,s13c)+tn2*u(i1,i2,i3,s23c)+tn3*u(i1,i2,i3,s33c)
1444 
1445  ! compute componenets of derived stress in tangential directions
1446  ss1d = sn1*s11t+sn2*s21t+sn3*s31t
1447  ss2d = sn1*s12t+sn2*s22t+sn3*s32t
1448  ss3d = sn1*s13t+sn2*s23t+sn3*s33t
1449  ts1d = tn1*s11t+tn2*s21t+tn3*s31t
1450  ts2d = tn1*s12t+tn2*s22t+tn3*s32t
1451  ts3d = tn1*s13t+tn2*s23t+tn3*s33t
1452 
1453  ! compute componenets of exact stress in tangential directions
1454  ss1e = sn1*s11e+sn2*s21e+sn3*s31e
1455  ss2e = sn1*s12e+sn2*s22e+sn3*s32e
1456  ss3e = sn1*s13e+sn2*s23e+sn3*s33e
1457  ts1e = tn1*s11e+tn2*s21e+tn3*s31e
1458  ts2e = tn1*s12e+tn2*s22e+tn3*s32e
1459  ts3e = tn1*s13e+tn2*s23e+tn3*s33e
1460 
1461  ss1 = ss1+ss1e-ss1d
1462  ss2 = ss2+ss2e-ss2d
1463  ss3 = ss3+ss3e-ss3d
1464  ts1 = ts1+ts1e-ts1d
1465  ts2 = ts2+ts2e-ts2d
1466  ts3 = ts3+ts3e-ts3d
1467 
1468  u(i1,i2,i3,s11c) = an1*ns1+sn1*ss1+tn1*ts1
1469  u(i1,i2,i3,s12c) = an1*ns2+sn1*ss2+tn1*ts2
1470  u(i1,i2,i3,s13c) = an1*ns3+sn1*ss3+tn1*ts3
1471 
1472  u(i1,i2,i3,s21c) = an2*ns1+sn2*ss1+tn2*ts1
1473  u(i1,i2,i3,s22c) = an2*ns2+sn2*ss2+tn2*ts2
1474  u(i1,i2,i3,s23c) = an2*ns3+sn2*ss3+tn2*ts3
1475 
1476  u(i1,i2,i3,s31c) = an3*ns1+sn3*ss1+tn3*ts1
1477  u(i1,i2,i3,s32c) = an3*ns2+sn3*ss2+tn3*ts2
1478  u(i1,i2,i3,s33c) = an3*ns3+sn3*ss3+tn3*ts3
1479 
1480 c u(i1,i2,i3,s11c) = s11e
1481 c u(i1,i2,i3,s12c) = s12e
1482 c u(i1,i2,i3,s13c) = s13e
1483 c
1484 c u(i1,i2,i3,s21c) = s21e
1485 c u(i1,i2,i3,s22c) = s22e
1486 c u(i1,i2,i3,s23c) = s23e
1487 c
1488 c u(i1,i2,i3,s31c) = s31e
1489 c u(i1,i2,i3,s32c) = s32e
1490 c u(i1,i2,i3,s33c) = s33e
1491  endLoopsMask3d()
1492 
1493  end if ! end gridType
1494 
1495  end if ! bc
1496  endLoopOverSides()
1497 c
1498 c.. substract off components of TZ force that were added twice
1499  beginEdgeMacro()
1500  if( bc1.eq.tractionBC .and. bc2.eq.tractionBC ) then
1501  if( gridType.eq.rectangular ) then
1502  beginLoopsMask3d()
1503  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
1504  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
1505  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
1506 
1507  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
1508  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
1509  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
1510 
1511  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
1512  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
1513  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
1514 
1515  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,s11e )
1516  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,s12e )
1517  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,s13e )
1518 
1519  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,s21e )
1520  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,s22e )
1521  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,s23e )
1522 
1523  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,s31e )
1524  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,s32e )
1525  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,s33e )
1526 
1527  if( edgeDirection.eq.0 ) then
1528  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)-(s11e-(kappa*u1x+lambda*(u2y+u3z)))
1529  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)-(s12e-mu*(u2x+u1y))
1530  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)-(s13e-mu*(u3x+u1z))
1531  else if( edgeDirection.eq.1 ) then
1532  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)-(s21e-mu*(u1y+u2x))
1533  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)-(s22e-(kappa*u2y+lambda*(u1x+u3z)))
1534  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)-(s23e-mu*(u3y+u2z))
1535  else
1536  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)-(s31e-mu*(u3x+u1z))
1537  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)-(s32e-mu*(u3y+u2z))
1538  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)-(s33e-(kappa*u3z+lambda*(u1x+u2y)))
1539  end if
1540  endLoopsMask3d()
1541  else
1542  beginLoopsMask3d()
1543  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1x )
1544  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1y )
1545  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,u1z )
1546 
1547  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2x )
1548  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2y )
1549  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,u2z )
1550 
1551  call ogDeriv( ep,0,1,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3x )
1552  call ogDeriv( ep,0,0,1,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3y )
1553  call ogDeriv( ep,0,0,0,1,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,u3z )
1554 
1555  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,s11e )
1556  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,s12e )
1557  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,s13e )
1558 
1559  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,s21e )
1560  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,s22e )
1561  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,s23e )
1562 
1563  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,s31e )
1564  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,s32e )
1565  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,s33e )
1566 
1567  tn1 = rx(i1,i2,i3,edgeDirection,0)
1568  tn2 = rx(i1,i2,i3,edgeDirection,1)
1569  tn3 = rx(i1,i2,i3,edgeDirection,2)
1570  aNormi = 1.0/max(epsx,sqrt(tn1**2+tn2**2+tn3**2))
1571  tn1 = tn1*aNormi
1572  tn2 = tn2*aNormi
1573  tn3 = tn3*aNormi
1574 
1575  f1 = tn1*(s11e-(kappa*u1x+lambda*(u2y+u3z)))+tn2*(s21e-mu*(u1y+u2x)) +tn3*(s31e-mu*(u3x+u1z))
1576  f2 = tn1*(s12e-mu*(u2x+u1y)) +tn2*(s22e-(kappa*u2y+lambda*(u1x+u3z)))+tn3*(s32e-mu*(u3y+u2z))
1577  f3 = tn1*(s13e-mu*(u3x+u1z)) +tn2*(s23e-mu*(u3y+u2z)) +tn3*(s33e-(kappa*u3z+lambda*(u1x+u2y)))
1578 
1579 c u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)-tn1*f1
1580 c u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)-tn1*f2
1581 c u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)-tn1*f3
1582 c
1583 c u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)-tn2*f1
1584 c u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)-tn2*f2
1585 c u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)-tn2*f3
1586 c
1587 c u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)-tn3*f1
1588 c u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)-tn3*f2
1589 c u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)-tn3*f3
1590 ccc
1591  u(i1,i2,i3,s11c) = s11e
1592  u(i1,i2,i3,s12c) = s12e
1593  u(i1,i2,i3,s13c) = s13e
1594 c
1595  u(i1,i2,i3,s21c) = s21e
1596  u(i1,i2,i3,s22c) = s22e
1597  u(i1,i2,i3,s23c) = s23e
1598 c
1599  u(i1,i2,i3,s31c) = s31e
1600  u(i1,i2,i3,s32c) = s32e
1601  u(i1,i2,i3,s33c) = s33e
1602  endLoopsMask3d()
1603  end if ! gridType
1604  end if ! bcTypes
1605  endEdgeMacro()
1606  end if
1607 c
1608 c.. re-compute stress at traction-traction edges
1609  beginEdgeMacro()
1610  if( gridType.eq.rectangular ) then
1611  ! do nothing ...
1612  else
1613  if( bc1.eq.tractionBC .and. bc2.eq.tractionBC ) then
1614  beginLoopsMask3d()
1615  if( edgeDirection.eq.0 ) then
1616  norm1(0) = -is2*rx(i1,i2,i3,1,0)
1617  norm1(1) = -is2*rx(i1,i2,i3,1,1)
1618  norm1(2) = -is2*rx(i1,i2,i3,1,2)
1619  f11 = bcf(side2,1,i1,i2,i3,s11c)
1620  f21 = bcf(side2,1,i1,i2,i3,s12c)
1621  f31 = bcf(side2,1,i1,i2,i3,s13c)
1622 
1623  norm2(0) = -is3*rx(i1,i2,i3,2,0)
1624  norm2(1) = -is3*rx(i1,i2,i3,2,1)
1625  norm2(2) = -is3*rx(i1,i2,i3,2,2)
1626  f12 = bcf(side3,2,i1,i2,i3,s11c)
1627  f22 = bcf(side3,2,i1,i2,i3,s12c)
1628  f32 = bcf(side3,2,i1,i2,i3,s13c)
1629  else if( edgeDirection.eq.1 ) then
1630  norm1(0) = -is3*rx(i1,i2,i3,2,0)
1631  norm1(1) = -is3*rx(i1,i2,i3,2,1)
1632  norm1(2) = -is3*rx(i1,i2,i3,2,2)
1633  f11 = bcf(side3,2,i1,i2,i3,s11c)
1634  f21 = bcf(side3,2,i1,i2,i3,s12c)
1635  f31 = bcf(side3,2,i1,i2,i3,s13c)
1636 
1637  norm2(0) = -is1*rx(i1,i2,i3,0,0)
1638  norm2(1) = -is1*rx(i1,i2,i3,0,1)
1639  norm2(2) = -is1*rx(i1,i2,i3,0,2)
1640  f12 = bcf(side1,0,i1,i2,i3,s11c)
1641  f22 = bcf(side1,0,i1,i2,i3,s12c)
1642  f32 = bcf(side1,0,i1,i2,i3,s13c)
1643  else
1644  norm1(0) = -is1*rx(i1,i2,i3,0,0)
1645  norm1(1) = -is1*rx(i1,i2,i3,0,1)
1646  norm1(2) = -is1*rx(i1,i2,i3,0,2)
1647  f11 = bcf(side1,0,i1,i2,i3,s11c)
1648  f21 = bcf(side1,0,i1,i2,i3,s12c)
1649  f31 = bcf(side1,0,i1,i2,i3,s13c)
1650 
1651  norm2(0) = -is2*rx(i1,i2,i3,1,0)
1652  norm2(1) = -is2*rx(i1,i2,i3,1,1)
1653  norm2(2) = -is2*rx(i1,i2,i3,1,2)
1654  f12 = bcf(side2,1,i1,i2,i3,s11c)
1655  f22 = bcf(side2,1,i1,i2,i3,s12c)
1656  f32 = bcf(side2,1,i1,i2,i3,s13c)
1657  end if
1658 
1659  aNormi = 1.0/max(epsx,sqrt(norm1(0)**2+norm1(1)**2+norm1(2)**2))
1660  norm1(0) = norm1(0)*aNormi
1661  norm1(1) = norm1(1)*aNormi
1662  norm1(2) = norm1(2)*aNormi
1663 
1664  aNormi = 1.0/max(epsx,sqrt(norm2(0)**2+norm2(1)**2+norm2(2)**2))
1665  norm2(0) = norm2(0)*aNormi
1666  norm2(1) = norm2(1)*aNormi
1667  norm2(2) = norm2(2)*aNormi
1668 
1669  b11 = f11-(norm1(0)*u(i1,i2,i3,s11c)+norm1(1)*u(i1,i2,i3,s21c)+norm1(2)*u(i1,i2,i3,s31c))
1670  b21 = f21-(norm1(0)*u(i1,i2,i3,s12c)+norm1(1)*u(i1,i2,i3,s22c)+norm1(2)*u(i1,i2,i3,s32c))
1671  b31 = f31-(norm1(0)*u(i1,i2,i3,s13c)+norm1(1)*u(i1,i2,i3,s23c)+norm1(2)*u(i1,i2,i3,s33c))
1672 
1673  dot1 = norm1(0)*norm2(0)+norm1(1)*norm2(1)+norm1(2)*norm2(2)
1674  dot2 = -sin(acos(dot1))
1675 
1676  b12 = (f12-(norm2(0)*u(i1,i2,i3,s11c)+norm2(1)*u(i1,i2,i3,s21c)+norm2(2)*u(i1,i2,i3,s31c))-dot1*b11)/dot2
1677  b22 = (f22-(norm2(0)*u(i1,i2,i3,s12c)+norm2(1)*u(i1,i2,i3,s22c)+norm2(2)*u(i1,i2,i3,s32c))-dot1*b21)/dot2
1678  b32 = (f32-(norm2(0)*u(i1,i2,i3,s13c)+norm2(1)*u(i1,i2,i3,s23c)+norm2(2)*u(i1,i2,i3,s33c))-dot1*b31)/dot2
1679 
1680  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)+norm1(0)*b11+norm2(0)*b12
1681  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)+norm1(0)*b21+norm2(0)*b22
1682  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)+norm1(0)*b31+norm2(0)*b32
1683 
1684  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)+norm1(1)*b11+norm2(1)*b12
1685  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)+norm1(1)*b21+norm2(1)*b22
1686  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)+norm1(1)*b31+norm2(1)*b32
1687 
1688  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)+norm1(2)*b11+norm2(2)*b12
1689  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)+norm1(2)*b21+norm2(2)*b22
1690  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)+norm1(2)*b31+norm2(2)*b32
1691 
1692  endLoopsMask3d()
1693  end if ! bcTypes
1694  end if ! gridType
1695  endEdgeMacro()
1696 c end if
1697 
1698 c.. set exact corner conditions for now
1699  if( setCornersWithTZ .and.twilightZone.ne.0 ) then ! *wdh* 090909
1700  write(*,'(" bcOptSmFOS3D: INFO set exact values on corners")')
1701  beginLoopOverCorners3d(0)
1702  ! *wdh* Need to check the boundaryConditions on the adjacent faces before applying these values:
1703  if( mask(i1,i2,i3).gt.0 ) then
1704  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,ue )
1705  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,ve )
1706  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,we )
1707 
1708  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1e )
1709  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2e )
1710  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3e )
1711 
1712  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11e )
1713  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21e )
1714  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31e )
1715 
1716  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12e )
1717  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22e )
1718  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32e )
1719 
1720  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13e )
1721  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23e )
1722  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33e )
1723 
1724  u(i1,i2,i3,uc) = ue
1725  u(i1,i2,i3,vc) = ve
1726  u(i1,i2,i3,wc) = we
1727 
1728  u(i1,i2,i3,v1c) = v1e
1729  u(i1,i2,i3,v2c) = v2e
1730  u(i1,i2,i3,v3c) = v3e
1731 
1732  u(i1,i2,i3,s11c) = tau11e
1733  u(i1,i2,i3,s21c) = tau21e
1734  u(i1,i2,i3,s31c) = tau31e
1735 
1736  u(i1,i2,i3,s12c) = tau12e
1737  u(i1,i2,i3,s22c) = tau22e
1738  u(i1,i2,i3,s32c) = tau32e
1739 
1740  u(i1,i2,i3,s13c) = tau13e
1741  u(i1,i2,i3,s23c) = tau23e
1742  u(i1,i2,i3,s33c) = tau33e
1743  end if
1744  endLoopOverCorners3d()
1745  end if
1746 c*******
1747 c******* re-extrapolation components of stress to first ghost line ********
1748 c*******
1749 
1750  beginLoopOverSides(numGhost,numGhost)
1751  if( boundaryCondition(side,axis).eq.tractionBC.or.boundaryCondition(side,axis).eq.slipWall ) then
1752  beginLoops3d()
1753  if( mask(i1,i2,i3).ne.0 ) then
1754  u(i1-is1,i2-is2,i3-is3,s11c) = extrap3(u,i1,i2,i3,s11c,is1,is2,is3)
1755  u(i1-is1,i2-is2,i3-is3,s12c) = extrap3(u,i1,i2,i3,s12c,is1,is2,is3)
1756  u(i1-is1,i2-is2,i3-is3,s13c) = extrap3(u,i1,i2,i3,s13c,is1,is2,is3)
1757 
1758  u(i1-is1,i2-is2,i3-is3,s21c) = extrap3(u,i1,i2,i3,s21c,is1,is2,is3)
1759  u(i1-is1,i2-is2,i3-is3,s22c) = extrap3(u,i1,i2,i3,s22c,is1,is2,is3)
1760  u(i1-is1,i2-is2,i3-is3,s23c) = extrap3(u,i1,i2,i3,s23c,is1,is2,is3)
1761 
1762  u(i1-is1,i2-is2,i3-is3,s31c) = extrap3(u,i1,i2,i3,s31c,is1,is2,is3)
1763  u(i1-is1,i2-is2,i3-is3,s32c) = extrap3(u,i1,i2,i3,s32c,is1,is2,is3)
1764  u(i1-is1,i2-is2,i3-is3,s33c) = extrap3(u,i1,i2,i3,s33c,is1,is2,is3)
1765  end if
1766  endLoops3d()
1767 
1768  end if ! bc
1769  endLoopOverSides()
1770 
1771 c.. set the corners to the exact twilight zone function for testing ... CHANGE ME!! ...
1772  if( .false. ) then
1773  beginLoopOverEdges3d(1)
1774  if( mask(i1,i2,i3).gt.0 ) then
1775  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,ue )
1776  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,ve )
1777  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,we )
1778 
1779  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1e )
1780  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2e )
1781  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3e )
1782 
1783  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11e )
1784  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21e )
1785  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31e )
1786 
1787  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12e )
1788  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22e )
1789  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32e )
1790 
1791  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13e )
1792  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23e )
1793  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33e )
1794 
1795  u(i1,i2,i3,uc) = ue
1796  u(i1,i2,i3,vc) = ve
1797  u(i1,i2,i3,wc) = we
1798 
1799  u(i1,i2,i3,v1c) = v1e
1800  u(i1,i2,i3,v2c) = v2e
1801  u(i1,i2,i3,v3c) = v3e
1802 
1803  u(i1,i2,i3,s11c) = tau11e
1804  u(i1,i2,i3,s21c) = tau21e
1805  u(i1,i2,i3,s31c) = tau31e
1806 
1807  u(i1,i2,i3,s12c) = tau12e
1808  u(i1,i2,i3,s22c) = tau22e
1809  u(i1,i2,i3,s32c) = tau32e
1810 
1811  u(i1,i2,i3,s13c) = tau13e
1812  u(i1,i2,i3,s23c) = tau23e
1813  u(i1,i2,i3,s33c) = tau33e
1814  end if
1815  endLoopOverEdges3d()
1816  end if
1817 
1818 c*******
1819 c******* Extrapolation to the second ghost line ********
1820 c*******
1821 
1822  beginLoopOverSides(numGhost,numGhost)
1823  if( boundaryCondition(side,axis).gt.0 ) then
1824  beginLoops3d()
1825  if( mask(i1,i2,i3).ne.0 ) then
1826  do n=0,numberOfComponents-1
1827  u(i1-2*is1,i2-2*is2,i3-2*is3,n)=extrap3(u,i1-is1,i2-is2,i3-is3,n,is1,is2,is3)
1828  end do
1829  end if
1830  endLoops3d()
1831  end if ! bc
1832  endLoopOverSides()
1833 
1834 c..extrapolate the 2nd ghost line near the corners.
1835  do side1=0,1
1836  i1 = gridIndexRange(side1,axis1)
1837  is1 = 1-2*side1
1838  do side2=0,1
1839  i2 = gridIndexRange(side2,axis2)
1840  is2 = 1-2*side2
1841  do side3=0,1
1842  i3 = gridIndexRange(side3,axis3)
1843  is3 = 1-2*side3
1844 
1845 c extrapolate in the i1 direction
1846  if( boundaryCondition(side1,axis1).gt.0 ) then
1847  if( mask(i1,i2,i3).ne.0 ) then
1848  do n=0,numberOfComponents-1
1849  u(i1-2*is1,i2-is2,i3-is3,n)=extrap3(u,i1-is1,i2-is2,i3-is3,n,is1,0,0)
1850  end do
1851  end if
1852  end if
1853 
1854 c extrapolate in the i2 direction
1855  if( boundaryCondition(side2,axis2).gt.0 ) then
1856  if( mask(i1,i2,i3).ne.0 ) then
1857  do n=0,numberOfComponents-1
1858  u(i1-is1,i2-2*is2,i3-is3,n)=extrap3(u,i1-is1,i2-is2,i3-is3,n,0,is2,0)
1859  end do
1860  end if
1861  end if
1862 
1863 c extrapolate in the i3 direction
1864  if( boundaryCondition(side3,axis3).gt.0 ) then
1865  if( mask(i1,i2,i3).ne.0 ) then
1866  do n=0,numberOfComponents-1
1867  u(i1-is1,i2-is2,i3-2*is3,n)=extrap3(u,i1-is1,i2-is2,i3-is3,n,0,0,is3)
1868  end do
1869  end if
1870  end if
1871 
1872 c extrapolate in the diagonal direction
1873  if( boundaryCondition(side1,axis1).gt.0.and.boundaryCondition(side2,axis2).gt.0.and.boundaryCondition(side3,axis3).gt.0) then
1874  if( mask(i1,i2,i3).ne.0 ) then
1875  do n=0,numberOfComponents-1
1876  u(i1-2*is1,i2-2*is2,i3-2*is3,n)=extrap3(u,i1-is1,i2-is2,i3-is3,n,is1,is2,is3)
1877  end do
1878  end if
1879  end if
1880  end do
1881  end do
1882  end do
1883 
1884  if( .false. ) then
1885 c.. set the corners to the exact twilight zone function for testing ... CHANGE ME!! ...
1886  beginLoopOverEdges3d(1)
1887  if( mask(i1,i2,i3).gt.0 ) then
1888  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,uc,ue )
1889  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,vc,ve )
1890  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,wc,we )
1891 
1892  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v1c,v1e )
1893  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v2c,v2e )
1894  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,v3c,v3e )
1895 
1896  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s11c,tau11e )
1897  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s21c,tau21e )
1898  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s31c,tau31e )
1899 
1900  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s12c,tau12e )
1901  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s22c,tau22e )
1902  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s32c,tau32e )
1903 
1904  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s13c,tau13e )
1905  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s23c,tau23e )
1906  call ogDeriv( ep,0,0,0,0,xy(i1,i2,i3,0),xy(i1,i2,i3,1),xy(i1,i2,i3,2),t,s33c,tau33e )
1907 
1908  u(i1,i2,i3,uc) = ue
1909  u(i1,i2,i3,vc) = ve
1910  u(i1,i2,i3,wc) = we
1911 
1912  u(i1,i2,i3,v1c) = v1e
1913  u(i1,i2,i3,v2c) = v2e
1914  u(i1,i2,i3,v3c) = v3e
1915 
1916  u(i1,i2,i3,s11c) = tau11e
1917  u(i1,i2,i3,s21c) = tau21e
1918  u(i1,i2,i3,s31c) = tau31e
1919 
1920  u(i1,i2,i3,s12c) = tau12e
1921  u(i1,i2,i3,s22c) = tau22e
1922  u(i1,i2,i3,s32c) = tau32e
1923 
1924  u(i1,i2,i3,s13c) = tau13e
1925  u(i1,i2,i3,s23c) = tau23e
1926  u(i1,i2,i3,s33c) = tau33e
1927  end if
1928  endLoopOverEdges3d()
1929  end if
1930 
Generated on Fri Jan 4 2013 10:19:15 for CG by   doxygen 1.8.3