CG  Version 25
bcOptSmFOS3DEdge.h
Go to the documentation of this file.
1 c*******
2 c******* Fix up components of stress along the edges
3 c*******
4  if( .false. )then ! *wdh* 090910 -- turn this off for now --------------------
5 
6  write(*,'(" bcOptSmFOS3DEdge: do NOT apply edge fixup, grid,gridType=",i4,i2)') grid,gridType
7 
8  else
9 
10  ! write(*,'(" bcOptSmFOS3DEdge: DO apply edge fixup, grid,gridType=",i4,i2)') grid,gridType
11 
12  beginEdgeMacro()
13  if( bc1.eq.displacementBC .and. bc2.eq.displacementBC ) then
14  if( gridType.eq.rectangular ) then
15  beginLoopsMask3d()
16 c if( mask(i1,i2,i3).gt.0 ) then
17 c write(6,*)i1,i2,i3
18  u1x = diffr1(uc,dx(0))
19  u1y = diffr2(uc,dx(1))
20  u1z = diffr3(uc,dx(2))
21 
22  u2x = diffr1(vc,dx(0))
23  u2y = diffr2(vc,dx(1))
24  u2z = diffr3(vc,dx(2))
25 
26  u3x = diffr1(wc,dx(0))
27  u3y = diffr2(wc,dx(1))
28  u3z = diffr3(wc,dx(2))
29 
30  u(i1,i2,i3,s11c) = kappa*u1x+lambda*(u2y+u3z)
31  u(i1,i2,i3,s21c) = mu*(u2x+u1y)
32  u(i1,i2,i3,s31c) = mu*(u3x+u1z)
33  u(i1,i2,i3,s12c) = mu*(u2x+u1y)
34  u(i1,i2,i3,s22c) = kappa*u2y+lambda*(u1x+u3z)
35  u(i1,i2,i3,s32c) = mu*(u3y+u2z)
36  u(i1,i2,i3,s13c) = mu*(u3x+u1z)
37  u(i1,i2,i3,s23c) = mu*(u3y+u2z)
38  u(i1,i2,i3,s33c) = kappa*u3z+lambda*(u1x+u2y)
39 c end if
40  endLoopsMask3d()
41  else
42  beginLoopsMask3d()
43 c if( mask(i1,i2,i3).gt.0 ) then
44  u1r = diffr1(uc,dr(0))
45  u1s = diffr2(uc,dr(1))
46  u1t = diffr3(uc,dr(2))
47 
48  u2r = diffr1(vc,dr(0))
49  u2s = diffr2(vc,dr(1))
50  u2t = diffr3(vc,dr(2))
51 
52  u3r = diffr1(wc,dr(0))
53  u3s = diffr2(wc,dr(1))
54  u3t = diffr3(wc,dr(2))
55 
56  u1x = u1r*rx(i1,i2,i3,0,0)+u1s*rx(i1,i2,i3,1,0)+u1t*rx(i1,i2,i3,2,0)
57  u1y = u1r*rx(i1,i2,i3,0,1)+u1s*rx(i1,i2,i3,1,1)+u1t*rx(i1,i2,i3,2,1)
58  u1z = u1r*rx(i1,i2,i3,0,2)+u1s*rx(i1,i2,i3,1,2)+u1t*rx(i1,i2,i3,2,2)
59 
60  u2x = u2r*rx(i1,i2,i3,0,0)+u2s*rx(i1,i2,i3,1,0)+u2t*rx(i1,i2,i3,2,0)
61  u2y = u2r*rx(i1,i2,i3,0,1)+u2s*rx(i1,i2,i3,1,1)+u2t*rx(i1,i2,i3,2,1)
62  u2z = u2r*rx(i1,i2,i3,0,2)+u2s*rx(i1,i2,i3,1,2)+u2t*rx(i1,i2,i3,2,2)
63 
64  u3x = u3r*rx(i1,i2,i3,0,0)+u3s*rx(i1,i2,i3,1,0)+u3t*rx(i1,i2,i3,2,0)
65  u3y = u3r*rx(i1,i2,i3,0,1)+u3s*rx(i1,i2,i3,1,1)+u3t*rx(i1,i2,i3,2,1)
66  u3z = u3r*rx(i1,i2,i3,0,2)+u3s*rx(i1,i2,i3,1,2)+u3t*rx(i1,i2,i3,2,2)
67 
68  u(i1,i2,i3,s11c) = kappa*u1x+lambda*(u2y+u3z)
69  u(i1,i2,i3,s21c) = mu*(u2x+u1y)
70  u(i1,i2,i3,s31c) = mu*(u3x+u1z)
71  u(i1,i2,i3,s12c) = mu*(u2x+u1y)
72  u(i1,i2,i3,s22c) = kappa*u2y+lambda*(u1x+u3z)
73  u(i1,i2,i3,s32c) = mu*(u3y+u2z)
74  u(i1,i2,i3,s13c) = mu*(u3x+u1z)
75  u(i1,i2,i3,s23c) = mu*(u3y+u2z)
76  u(i1,i2,i3,s33c) = kappa*u3z+lambda*(u1x+u2y)
77 c end if
78  endLoopsMask3d()
79  end if ! gridType
80 
81  ! *wdh* we need to adjust for TZ here (not in a separate edgeMacro loop for then
82  ! corner points are corrected twice)
83  if( twilightZone.ne.0 ) then
84  beginLoopsMask3d()
85 c if( mask(i1,i2,i3).gt.0 ) then
86  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 )
87  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 )
88  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 )
89 
90  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 )
91  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 )
92  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 )
93 
94  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 )
95  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 )
96  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 )
97 
98  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 )
99  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 )
100  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 )
101 
102  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 )
103  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 )
104  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 )
105 
106  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 )
107  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 )
108  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 )
109 
110  tau11 = kappa*u1x+lambda*(u2y+u3z)
111  tau21 = mu*(u2x+u1y)
112  tau31 = mu*(u3x+u1z)
113  tau12 = tau21
114  tau22 = kappa*u2y+lambda*(u1x+u3z)
115  tau32 = mu*(u3y+u2z)
116  tau13 = tau31
117  tau23 = tau32
118  tau33 = kappa*u3z+lambda*(u1x+u2y)
119 
120  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)-tau11+s11e
121  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)-tau21+s21e
122  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)-tau31+s31e
123 
124  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)-tau12+s12e
125  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)-tau22+s22e
126  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)-tau32+s32e
127 
128  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)-tau13+s13e
129  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)-tau23+s23e
130  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)-tau33+s33e
131 
132 c u(i1,i2,i3,s11c) = s11e
133 c u(i1,i2,i3,s21c) = s21e
134 c u(i1,i2,i3,s31c) = s31e
135 c
136 c u(i1,i2,i3,s12c) = s12e
137 c u(i1,i2,i3,s22c) = s22e
138 c u(i1,i2,i3,s32c) = s32e
139 c
140 c u(i1,i2,i3,s13c) = s13e
141 c u(i1,i2,i3,s23c) = s23e
142 c u(i1,i2,i3,s33c) = s33e
143 c end if ! end if mask
144  endLoopsMask3d()
145  end if ! end if TZ
146 
147 
148  else if( bc1.eq.tractionBC .and. bc2.eq.tractionBC ) then
149  if( gridType.eq.rectangular ) then
150  ! do nothing because normals are perpendicular and so no part of the force is counted twice
151  else
152  ! do stuff because normals are not perpendicular and some part of the stress might be counted twice
153  beginLoopsMask3d()
154  if( edgeDirection.eq.0 ) then
155  norm1(0) = -is2*rx(i1,i2,i3,1,0)
156  norm1(1) = -is2*rx(i1,i2,i3,1,1)
157  norm1(2) = -is2*rx(i1,i2,i3,1,2)
158  f11 = bcf(side2,1,i1,i2,i3,s11c)
159  f21 = bcf(side2,1,i1,i2,i3,s12c)
160  f31 = bcf(side2,1,i1,i2,i3,s13c)
161 
162  norm2(0) = -is3*rx(i1,i2,i3,2,0)
163  norm2(1) = -is3*rx(i1,i2,i3,2,1)
164  norm2(2) = -is3*rx(i1,i2,i3,2,2)
165  f12 = bcf(side3,2,i1,i2,i3,s11c)
166  f22 = bcf(side3,2,i1,i2,i3,s12c)
167  f32 = bcf(side3,2,i1,i2,i3,s13c)
168  else if( edgeDirection.eq.1 ) then
169  norm1(0) = -is3*rx(i1,i2,i3,2,0)
170  norm1(1) = -is3*rx(i1,i2,i3,2,1)
171  norm1(2) = -is3*rx(i1,i2,i3,2,2)
172  f11 = bcf(side3,2,i1,i2,i3,s11c)
173  f21 = bcf(side3,2,i1,i2,i3,s12c)
174  f31 = bcf(side3,2,i1,i2,i3,s13c)
175 
176  norm2(0) = -is1*rx(i1,i2,i3,0,0)
177  norm2(1) = -is1*rx(i1,i2,i3,0,1)
178  norm2(2) = -is1*rx(i1,i2,i3,0,2)
179  f12 = bcf(side1,0,i1,i2,i3,s11c)
180  f22 = bcf(side1,0,i1,i2,i3,s12c)
181  f32 = bcf(side1,0,i1,i2,i3,s13c)
182  else
183  norm1(0) = -is1*rx(i1,i2,i3,0,0)
184  norm1(1) = -is1*rx(i1,i2,i3,0,1)
185  norm1(2) = -is1*rx(i1,i2,i3,0,2)
186  f11 = bcf(side1,0,i1,i2,i3,s11c)
187  f21 = bcf(side1,0,i1,i2,i3,s12c)
188  f31 = bcf(side1,0,i1,i2,i3,s13c)
189 
190  norm2(0) = -is2*rx(i1,i2,i3,1,0)
191  norm2(1) = -is2*rx(i1,i2,i3,1,1)
192  norm2(2) = -is2*rx(i1,i2,i3,1,2)
193  f12 = bcf(side2,1,i1,i2,i3,s11c)
194  f22 = bcf(side2,1,i1,i2,i3,s12c)
195  f32 = bcf(side2,1,i1,i2,i3,s13c)
196  end if
197 
198  aNormi = 1.0/max(epsx,sqrt(norm1(0)**2+norm1(1)**2+norm1(2)**2))
199  norm1(0) = norm1(0)*aNormi
200  norm1(1) = norm1(1)*aNormi
201  norm1(2) = norm1(2)*aNormi
202 
203  aNormi = 1.0/max(epsx,sqrt(norm2(0)**2+norm2(1)**2+norm2(2)**2))
204  norm2(0) = norm2(0)*aNormi
205  norm2(1) = norm2(1)*aNormi
206  norm2(2) = norm2(2)*aNormi
207 
208  b11 = f11-(norm1(0)*u(i1,i2,i3,s11c)+norm1(1)*u(i1,i2,i3,s21c)+norm1(2)*u(i1,i2,i3,s31c))
209  b21 = f21-(norm1(0)*u(i1,i2,i3,s12c)+norm1(1)*u(i1,i2,i3,s22c)+norm1(2)*u(i1,i2,i3,s32c))
210  b31 = f31-(norm1(0)*u(i1,i2,i3,s13c)+norm1(1)*u(i1,i2,i3,s23c)+norm1(2)*u(i1,i2,i3,s33c))
211 
212  dot1 = norm1(0)*norm2(0)+norm1(1)*norm2(1)+norm1(2)*norm2(2)
213  dot2 = -sin(acos(dot1))
214 
215  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
216  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
217  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
218 
219  u(i1,i2,i3,s11c) = u(i1,i2,i3,s11c)+norm1(0)*b11+norm2(0)*b12
220  u(i1,i2,i3,s12c) = u(i1,i2,i3,s12c)+norm1(0)*b21+norm2(0)*b22
221  u(i1,i2,i3,s13c) = u(i1,i2,i3,s13c)+norm1(0)*b31+norm2(0)*b32
222 
223  u(i1,i2,i3,s21c) = u(i1,i2,i3,s21c)+norm1(1)*b11+norm2(1)*b12
224  u(i1,i2,i3,s22c) = u(i1,i2,i3,s22c)+norm1(1)*b21+norm2(1)*b22
225  u(i1,i2,i3,s23c) = u(i1,i2,i3,s23c)+norm1(1)*b31+norm2(1)*b32
226 
227  u(i1,i2,i3,s31c) = u(i1,i2,i3,s31c)+norm1(2)*b11+norm2(2)*b12
228  u(i1,i2,i3,s32c) = u(i1,i2,i3,s32c)+norm1(2)*b21+norm2(2)*b22
229  u(i1,i2,i3,s33c) = u(i1,i2,i3,s33c)+norm1(2)*b31+norm2(2)*b32
230  endLoopsMask3d()
231  end if ! gridType
232 
233  end if ! bcType
234  endEdgeMacro()
235 
236  end if
Generated on Fri Jan 4 2013 10:19:15 for CG by   doxygen 1.8.3