CG  Version 25
bcExtended3d4.h
Go to the documentation of this file.
1  ! results from bcExtended3d4.maple
2  ! Assign values on the extended boundary next to two PEC boundaries
3  !
4  ! Here we assume the following are defined
5  ! c11,c22,c33,c1,c2,c3
6  ! urr,uss,utt,ur,us,ut (also for v and w)
7  ! deltaFu,deltaFv,deltaFw = RHS for Delta(u,v,w)
8  ! g1f,g2f = RHS for extrapolation, a1.D+2^4u(i1,i2-2)=g1f, a2.D+2^4u(i1-2,i2)=g2f,
9  !
10  DeltaU = c11*urr+c22*uss+c33*utt+c1*ur+c2*us+c3*ut - deltaFu
11  DeltaV = c11*vrr+c22*vss+c33*vtt+c1*vr+c2*vs+c3*vt - deltaFv
12  DeltaW = c11*wrr+c22*wss+c33*wtt+c1*wr+c2*ws+c3*wt - deltaFw
13 
14 ! ** decompose point u(i1-is1,i2-is2,i3-is3) into components along a1,a2,a3 **
15  a11c=A11D3(i1-is1,i2-is2,i3-is3)
16  a12c=A12D3(i1-is1,i2-is2,i3-is3)
17  a13c=A13D3(i1-is1,i2-is2,i3-is3)
18  a21c=A21D3(i1-is1,i2-is2,i3-is3)
19  a22c=A22D3(i1-is1,i2-is2,i3-is3)
20  a23c=A23D3(i1-is1,i2-is2,i3-is3)
21  a31c=A31D3(i1-is1,i2-is2,i3-is3)
22  a32c=A32D3(i1-is1,i2-is2,i3-is3)
23  a33c=A33D3(i1-is1,i2-is2,i3-is3)
24 
25  a1a1=a11c*a11c+a12c*a12c+a13c*a13c
26  a1a2=a11c*a21c+a12c*a22c+a13c*a23c
27  a1a3=a11c*a31c+a12c*a32c+a13c*a33c
28  a2a3=a21c*a31c+a22c*a32c+a23c*a33c
29  a3a3=a31c*a31c+a32c*a32c+a33c*a33c
30 ! The tangential component is assumed valid:
31  a1Dotu1=a11c*u(i1-is1,i2-is2,i3-is3,ex)+a12c*u(i1-is1,i2-is2,i3-is3,ey)+a13c*u(i1-is1,i2-is2,i3-is3,ez)
32  a3Dotu1=a31c*u(i1-is1,i2-is2,i3-is3,ex)+a32c*u(i1-is1,i2-is2,i3-is3,ey)+a33c*u(i1-is1,i2-is2,i3-is3,ez)
33  ! u(i1-is1,i2-is2,i3-is3,k) = b1[k]*x1 +g1[k]
34  b11 =-a11c*(-a1a2*a3a3+a1a3*a2a3)/(-a1a1*a3a3+a1a3**2)+a21c-a31c*(-a1a1*a2a3+a1a3*a1a2)/(-a1a1*a3a3+a1a3**2)
35  b12 =-a12c*(-a1a2*a3a3+a1a3*a2a3)/(-a1a1*a3a3+a1a3**2)+a22c-a32c*(-a1a1*a2a3+a1a3*a1a2)/(-a1a1*a3a3+a1a3**2)
36  b13 =-a13c*(-a1a2*a3a3+a1a3*a2a3)/(-a1a1*a3a3+a1a3**2)+a23c-a33c*(-a1a1*a2a3+a1a3*a1a2)/(-a1a1*a3a3+a1a3**2)
37  g11 =-(-a11c*a1a3*a3Dotu1+a11c*a1Dotu1*a3a3-a31c*a1a3*a1Dotu1+a31c*a1a1*a3Dotu1)/(-a1a1*a3a3+a1a3**2)
38  g12 =-(-a12c*a1a3*a3Dotu1+a12c*a1Dotu1*a3a3-a32c*a1a3*a1Dotu1+a32c*a1a1*a3Dotu1)/(-a1a1*a3a3+a1a3**2)
39  g13 =(a13c*a1a3*a3Dotu1-a13c*a1Dotu1*a3a3+a33c*a1a3*a1Dotu1-a33c*a1a1*a3Dotu1)/(-a1a1*a3a3+a1a3**2)
40 
41 ! ** decompose point u(i1-2*is1,i2-2*is2,i3-2*is3) into components along a1,a2,a3 **
42  a11c=A11D3(i1-2*is1,i2-2*is2,i3-2*is3)
43  a12c=A12D3(i1-2*is1,i2-2*is2,i3-2*is3)
44  a13c=A13D3(i1-2*is1,i2-2*is2,i3-2*is3)
45  a21c=A21D3(i1-2*is1,i2-2*is2,i3-2*is3)
46  a22c=A22D3(i1-2*is1,i2-2*is2,i3-2*is3)
47  a23c=A23D3(i1-2*is1,i2-2*is2,i3-2*is3)
48  a31c=A31D3(i1-2*is1,i2-2*is2,i3-2*is3)
49  a32c=A32D3(i1-2*is1,i2-2*is2,i3-2*is3)
50  a33c=A33D3(i1-2*is1,i2-2*is2,i3-2*is3)
51 
52  a1a1=a11c*a11c+a12c*a12c+a13c*a13c
53  a1a2=a11c*a21c+a12c*a22c+a13c*a23c
54  a1a3=a11c*a31c+a12c*a32c+a13c*a33c
55  a2a3=a21c*a31c+a22c*a32c+a23c*a33c
56  a3a3=a31c*a31c+a32c*a32c+a33c*a33c
57 ! The tangential component is assumed valid:
58  a1Dotu2=a11c*u(i1-2*is1,i2-2*is2,i3-2*is3,ex)+a12c*u(i1-2*is1,i2-2*is2,i3-2*is3,ey)+a13c*u(i1-2*is1,i2-2*is2,i3-2*is3,ez)
59  a3Dotu2=a31c*u(i1-2*is1,i2-2*is2,i3-2*is3,ex)+a32c*u(i1-2*is1,i2-2*is2,i3-2*is3,ey)+a33c*u(i1-2*is1,i2-2*is2,i3-2*is3,ez)
60  ! u(i1-2*is1,i2-2*is2,i3-2*is3,k) = b2[k]*x2 +g2[k]
61  b21 =-a11c*(-a1a2*a3a3+a1a3*a2a3)/(-a1a1*a3a3+a1a3**2)+a21c-a31c*(-a1a1*a2a3+a1a3*a1a2)/(-a1a1*a3a3+a1a3**2)
62  b22 =-a12c*(-a1a2*a3a3+a1a3*a2a3)/(-a1a1*a3a3+a1a3**2)+a22c-a32c*(-a1a1*a2a3+a1a3*a1a2)/(-a1a1*a3a3+a1a3**2)
63  b23 =-a13c*(-a1a2*a3a3+a1a3*a2a3)/(-a1a1*a3a3+a1a3**2)+a23c-a33c*(-a1a1*a2a3+a1a3*a1a2)/(-a1a1*a3a3+a1a3**2)
64  g21 =(a11c*a1a3*a3Dotu2-a11c*a1Dotu2*a3a3+a31c*a1a3*a1Dotu2-a31c*a1a1*a3Dotu2)/(-a1a1*a3a3+a1a3**2)
65  g22 =(a12c*a1a3*a3Dotu2-a12c*a1Dotu2*a3a3+a32c*a1a3*a1Dotu2-a32c*a1a1*a3Dotu2)/(-a1a1*a3a3+a1a3**2)
66  g23 =(a13c*a1a3*a3Dotu2-a13c*a1Dotu2*a3a3+a33c*a1a3*a1Dotu2-a33c*a1a1*a3Dotu2)/(-a1a1*a3a3+a1a3**2)
67 
68 ! ** decompose point u(i1-js1,i2-js2,i3-js3) into components along a1,a2,a3 **
69  a11c=A11D3(i1-js1,i2-js2,i3-js3)
70  a12c=A12D3(i1-js1,i2-js2,i3-js3)
71  a13c=A13D3(i1-js1,i2-js2,i3-js3)
72  a21c=A21D3(i1-js1,i2-js2,i3-js3)
73  a22c=A22D3(i1-js1,i2-js2,i3-js3)
74  a23c=A23D3(i1-js1,i2-js2,i3-js3)
75  a31c=A31D3(i1-js1,i2-js2,i3-js3)
76  a32c=A32D3(i1-js1,i2-js2,i3-js3)
77  a33c=A33D3(i1-js1,i2-js2,i3-js3)
78 
79  a2a2=a21c*a21c+a22c*a22c+a23c*a23c
80  a1a2=a11c*a21c+a12c*a22c+a13c*a23c
81  a1a3=a11c*a31c+a12c*a32c+a13c*a33c
82  a2a3=a21c*a31c+a22c*a32c+a23c*a33c
83  a3a3=a31c*a31c+a32c*a32c+a33c*a33c
84 ! The tangential component is assumed valid:
85  a2Dotu3=a21c*u(i1-js1,i2-js2,i3-js3,ex)+a22c*u(i1-js1,i2-js2,i3-js3,ey)+a23c*u(i1-js1,i2-js2,i3-js3,ez)
86  a3Dotu3=a31c*u(i1-js1,i2-js2,i3-js3,ex)+a32c*u(i1-js1,i2-js2,i3-js3,ey)+a33c*u(i1-js1,i2-js2,i3-js3,ez)
87  ! u(i1-js1,i2-js2,i3-js3,k) = b3[k]*x3 +g3[k]
88  b31 =a11c-a21c/(a2a3**2-a3a3*a2a2)*(-a1a2*a3a3+a1a3*a2a3)+a31c*(-a1a2*a2a3+a2a2*a1a3)/(a2a3**2-a3a3*a2a2)
89  b32 =a12c-a22c/(a2a3**2-a3a3*a2a2)*(-a1a2*a3a3+a1a3*a2a3)+a32c*(-a1a2*a2a3+a2a2*a1a3)/(a2a3**2-a3a3*a2a2)
90  b33 =a13c-a23c/(a2a3**2-a3a3*a2a2)*(-a1a2*a3a3+a1a3*a2a3)+a33c*(-a1a2*a2a3+a2a2*a1a3)/(a2a3**2-a3a3*a2a2)
91  g31 =(-a21c*a3a3*a2Dotu3+a21c*a2a3*a3Dotu3-a31c*a2a2*a3Dotu3+a31c*a2Dotu3*a2a3)/(a2a3**2-a3a3*a2a2)
92  g32 =(-a22c*a3a3*a2Dotu3+a22c*a2a3*a3Dotu3-a32c*a2a2*a3Dotu3+a32c*a2Dotu3*a2a3)/(a2a3**2-a3a3*a2a2)
93  g33 =(-a23c*a3a3*a2Dotu3+a23c*a2a3*a3Dotu3-a33c*a2a2*a3Dotu3+a33c*a2Dotu3*a2a3)/(a2a3**2-a3a3*a2a2)
94 
95 ! ** decompose point u(i1-2*js1,i2-2*js2,i3-2*js3) into components along a1,a2,a3 **
96  a11c=A11D3(i1-2*js1,i2-2*js2,i3-2*js3)
97  a12c=A12D3(i1-2*js1,i2-2*js2,i3-2*js3)
98  a13c=A13D3(i1-2*js1,i2-2*js2,i3-2*js3)
99  a21c=A21D3(i1-2*js1,i2-2*js2,i3-2*js3)
100  a22c=A22D3(i1-2*js1,i2-2*js2,i3-2*js3)
101  a23c=A23D3(i1-2*js1,i2-2*js2,i3-2*js3)
102  a31c=A31D3(i1-2*js1,i2-2*js2,i3-2*js3)
103  a32c=A32D3(i1-2*js1,i2-2*js2,i3-2*js3)
104  a33c=A33D3(i1-2*js1,i2-2*js2,i3-2*js3)
105 
106  a2a2=a21c*a21c+a22c*a22c+a23c*a23c
107  a1a2=a11c*a21c+a12c*a22c+a13c*a23c
108  a1a3=a11c*a31c+a12c*a32c+a13c*a33c
109  a2a3=a21c*a31c+a22c*a32c+a23c*a33c
110  a3a3=a31c*a31c+a32c*a32c+a33c*a33c
111 ! The tangential component is assumed valid:
112  a2Dotu4=a21c*u(i1-2*js1,i2-2*js2,i3-2*js3,ex)+a22c*u(i1-2*js1,i2-2*js2,i3-2*js3,ey)+a23c*u(i1-2*js1,i2-2*js2,i3-2*js3,ez)
113  a3Dotu4=a31c*u(i1-2*js1,i2-2*js2,i3-2*js3,ex)+a32c*u(i1-2*js1,i2-2*js2,i3-2*js3,ey)+a33c*u(i1-2*js1,i2-2*js2,i3-2*js3,ez)
114  ! u(i1-2*js1,i2-2*js2,i3-2*js3,k) = b4[k]*x4 +g4[k]
115  b41 =a11c-a21c/(a2a3**2-a3a3*a2a2)*(-a1a2*a3a3+a1a3*a2a3)+a31c*(-a1a2*a2a3+a2a2*a1a3)/(a2a3**2-a3a3*a2a2)
116  b42 =a12c-a22c/(a2a3**2-a3a3*a2a2)*(-a1a2*a3a3+a1a3*a2a3)+a32c*(-a1a2*a2a3+a2a2*a1a3)/(a2a3**2-a3a3*a2a2)
117  b43 =a13c-a23c/(a2a3**2-a3a3*a2a2)*(-a1a2*a3a3+a1a3*a2a3)+a33c*(-a1a2*a2a3+a2a2*a1a3)/(a2a3**2-a3a3*a2a2)
118  g41 =-(a21c*a3a3*a2Dotu4-a21c*a3Dotu4*a2a3+a31c*a2a2*a3Dotu4-a31c*a2a3*a2Dotu4)/(a2a3**2-a3a3*a2a2)
119  g42 =-(a22c*a3a3*a2Dotu4-a22c*a3Dotu4*a2a3+a32c*a2a2*a3Dotu4-a32c*a2a3*a2Dotu4)/(a2a3**2-a3a3*a2a2)
120  g43 =(-a23c*a3a3*a2Dotu4+a23c*a3Dotu4*a2a3-a33c*a2a2*a3Dotu4+a33c*a2a3*a2Dotu4)/(a2a3**2-a3a3*a2a2)
121 
122  ! Evaluate a1, a2 and a3 at the corner
123  a11=A11D3(i1,i2,i3)
124  a12=A12D3(i1,i2,i3)
125  a13=A13D3(i1,i2,i3)
126  a21=A21D3(i1,i2,i3)
127  a22=A22D3(i1,i2,i3)
128  a23=A23D3(i1,i2,i3)
129  a31=A31D3(i1,i2,i3)
130  a32=A32D3(i1,i2,i3)
131  a33=A33D3(i1,i2,i3)
132 
133  a1DotLu=a11*DeltaU+a12*DeltaV+a13*DeltaW
134  a2DotLu=a21*DeltaU+a22*DeltaV+a23*DeltaW
135 
136 ! a1.Lu = 0
137 ! e1 := cc11a*u(i1-2,i2,i3)+cc12a*v(i1-2,i2,i3)+cc13a*w(i1-2,i2,i3)
138 ! + cc14a*u(i1-1,i2,i3)+cc15a*v(i1-1,i2,i3)+cc16a*w(i1-1,i2,i3)
139 ! + cc11b*u(i1,i2-2,i3)+cc12b*v(i1,i2-2,i3)+cc13b*w(i1,i2-2,i3)
140 ! + cc14b*u(i1,i2-1,i3)+cc15b*v(i1,i2-1,i3)+cc16b*w(i1,i2-1,i3) - f1:
141 ! a2.Lu = 0 :
142 ! e2 := cc21a*u(i1-2,i2,i3)+cc22a*v(i1-2,i2,i3)+cc23a*w(i1-2,i2,i3)
143 ! + cc24a*u(i1-1,i2,i3)+cc25a*v(i1-1,i2,i3)+cc26a*w(i1-1,i2,i3)
144 ! + cc21b*u(i1,i2-2,i3)+cc22b*v(i1,i2-2,i3)+cc23b*w(i1,i2-2,i3)
145 ! + cc24b*u(i1,i2-1,i3)+cc25b*v(i1,i2-1,i3)+cc26b*w(i1,i2-1,i3) - f2:
146  cc11a=(-1/12.*c11/dra**2+1/12.*c1/dra)*a11
147  cc12a=(-1/12.*c11/dra**2+1/12.*c1/dra)*a12
148  cc13a=(-1/12.*c11/dra**2+1/12.*c1/dra)*a13
149  cc14a=(4/3.*c11/dra**2-2/3.*c1/dra)*a11
150  cc15a=(4/3.*c11/dra**2-2/3.*c1/dra)*a12
151  cc16a=(4/3.*c11/dra**2-2/3.*c1/dra)*a13
152  cc11b=(-1/12.*c22/dsa**2+1/12.*c2/dsa)*a11
153  cc12b=(-1/12.*c22/dsa**2+1/12.*c2/dsa)*a12
154  cc13b=(-1/12.*c22/dsa**2+1/12.*c2/dsa)*a13
155  cc14b=(4/3.*c22/dsa**2-2/3.*c2/dsa)*a11
156  cc15b=(4/3.*c22/dsa**2-2/3.*c2/dsa)*a12
157  cc16b=(4/3.*c22/dsa**2-2/3.*c2/dsa)*a13
158  cc21a=(-1/12.*c11/dra**2+1/12.*c1/dra)*a21
159  cc22a=(-1/12.*c11/dra**2+1/12.*c1/dra)*a22
160  cc23a=(-1/12.*c11/dra**2+1/12.*c1/dra)*a23
161  cc24a=(4/3.*c11/dra**2-2/3.*c1/dra)*a21
162  cc25a=(4/3.*c11/dra**2-2/3.*c1/dra)*a22
163  cc26a=(4/3.*c11/dra**2-2/3.*c1/dra)*a23
164  cc21b=(-1/12.*c22/dsa**2+1/12.*c2/dsa)*a21
165  cc22b=(-1/12.*c22/dsa**2+1/12.*c2/dsa)*a22
166  cc23b=(-1/12.*c22/dsa**2+1/12.*c2/dsa)*a23
167  cc24b=(4/3.*c22/dsa**2-2/3.*c2/dsa)*a21
168  cc25b=(4/3.*c22/dsa**2-2/3.*c2/dsa)*a22
169  cc26b=(4/3.*c22/dsa**2-2/3.*c2/dsa)*a23
170 
171  f1=a1DotLu-cc11a*u(i1-2*is1,i2-2*is2,i3-2*is3,ex)-cc12a*u(i1-2*is1,i2-2*is2,i3-2*is3,ey)-cc13a*u(i1-2*is1,i2-2*is2,i3-2*is3,ez)-cc14a*u(i1-is1,i2-is2,i3-is3,ex)-cc15a*u(i1-is1,i2-is2,i3-is3,ey)-cc16a*u(i1-is1,i2-is2,i3-is3,ez)-cc11b*u(i1-2*js1,i2-2*js2,i3-2*js3,ex)-cc12b*u(i1-2*js1,i2-2*js2,i3-2*js3,ey)-cc13b*u(i1-2*js1,i2-2*js2,i3-2*js3,ez)-cc14b*u(i1-js1,i2-js2,i3-js3,ex)-cc15b*u(i1-js1,i2-js2,i3-js3,ey)-cc16b*u(i1-js1,i2-js2,i3-js3,ez)
172  f2=a2DotLu-cc21a*u(i1-2*is1,i2-2*is2,i3-2*is3,ex)-cc22a*u(i1-2*is1,i2-2*is2,i3-2*is3,ey)-cc23a*u(i1-2*is1,i2-2*is2,i3-2*is3,ez)-cc24a*u(i1-is1,i2-is2,i3-is3,ex)-cc25a*u(i1-is1,i2-is2,i3-is3,ey)-cc26a*u(i1-is1,i2-is2,i3-is3,ez)-cc21b*u(i1-2*js1,i2-2*js2,i3-2*js3,ex)-cc22b*u(i1-2*js1,i2-2*js2,i3-2*js3,ey)-cc23b*u(i1-2*js1,i2-2*js2,i3-2*js3,ez)-cc24b*u(i1-js1,i2-js2,i3-js3,ex)-cc25b*u(i1-js1,i2-js2,i3-js3,ey)-cc26b*u(i1-js1,i2-js2,i3-js3,ez)
173  f3=6*a21*u(i1,i2,i3,ex)+6*a22*u(i1,i2,i3,ey)+6*a23*u(i1,i2,i3,ez)-4*a21*u(i1+is1,i2+is2,i3+is3,ex)-4*a22*u(i1+is1,i2+is2,i3+is3,ey)-4*a23*u(i1+is1,i2+is2,i3+is3,ez)+a21*u(i1+2*is1,i2+2*is2,i3+2*is3,ex)+a22*u(i1+2*is1,i2+2*is2,i3+2*is3,ey)+a23*u(i1+2*is1,i2+2*is2,i3+2*is3,ez)-g2f
174  f4=6*a11*u(i1,i2,i3,ex)+6*a12*u(i1,i2,i3,ey)+6*a13*u(i1,i2,i3,ez)-4*a11*u(i1+js1,i2+js2,i3+js3,ex)-4*a12*u(i1+js1,i2+js2,i3+js3,ey)-4*a13*u(i1+js1,i2+js2,i3+js3,ez)+a11*u(i1+2*js1,i2+2*js2,i3+2*js3,ex)+a12*u(i1+2*js1,i2+2*js2,i3+2*js3,ey)+a13*u(i1+2*js1,i2+2*js2,i3+2*js3,ez)-g1f
175 
176  ! Simplfied forms for the 4 equations a1.Lu, a2.Lu, a2.D+r4 u = g2f a1.D+s4 u = g1f
177  ! e1x := dd11*x1+dd12*x2+dd13*x3+dd14*x4+ f1x
178  ! e2x := dd21*x1+dd22*x2+dd23*x3+dd24*x4+ f2x
179  ! e3x := dd31*x1+dd32*x2+dd33*x3+dd34*x4+ f3x
180  ! e4x := dd41*x1+dd42*x2+dd43*x3+dd44*x4+ f4x
181  dd11=cc14a*b11+cc15a*b12+cc16a*b13
182  dd12=cc11a*b21+cc12a*b22+cc13a*b23
183  dd13=cc14b*b31+cc15b*b32+cc16b*b33
184  dd14=cc11b*b41+cc12b*b42+cc13b*b43
185  dd21=cc24a*b11+cc25a*b12+cc26a*b13
186  dd22=cc21a*b21+cc22a*b22+cc23a*b23
187  dd23=cc24b*b31+cc25b*b32+cc26b*b33
188  dd24=cc21b*b41+cc22b*b42+cc23b*b43
189  dd31=-4*a21*b11-4*a22*b12-4*a23*b13
190  dd32=a21*b21+a22*b22+a23*b23
191  dd33=0
192  dd34=0
193  dd41=0
194  dd42=0
195  dd43=-4*a11*b31-4*a12*b32-4*a13*b33
196  dd44=a11*b41+a12*b42+a13*b43
197 
198  f1x=cc11a*g21+cc12a*g22+cc13a*g23+cc14a*g11+cc15a*g12+cc16a*g13+cc11b*g41+cc12b*g42+cc13b*g43+cc14b*g31+cc15b*g32+cc16b*g33+f1
199  f2x=cc21a*g21+cc22a*g22+cc23a*g23+cc24a*g11+cc25a*g12+cc26a*g13+cc21b*g41+cc22b*g42+cc23b*g43+cc24b*g31+cc25b*g32+cc26b*g33+f2
200  f3x=a21*g21+a22*g22+a23*g23-4*a21*g11-4*a22*g12-4*a23*g13+f3
201  f4x=a11*g41+a12*g42+a13*g43-4*a11*g31-4*a12*g32-4*a13*g33+f4
202 
203 ! solution x1,x2,x3,x4:
204  det=-dd32*dd43*dd14*dd21-dd32*dd11*dd23*dd44+dd32*dd11*dd43*dd24+dd43*dd14*dd22*dd31-dd13*dd44*dd22*dd31+dd12*dd31*dd23*dd44+dd32*dd13*dd44*dd21-dd12*dd31*dd43*dd24
205  x1=(-dd32*f2x*dd13*dd44-dd32*dd43*dd24*f1x+dd32*dd23*dd44*f1x+dd32*dd43*f2x*dd14-dd32*dd23*f4x*dd14+dd32*dd24*dd13*f4x-dd23*dd44*dd12*f3x+dd22*f3x*dd13*dd44-dd43*dd22*f3x*dd14+dd43*dd24*dd12*f3x)/det
206  x2=(dd31*f2x*dd13*dd44+dd31*dd43*dd24*f1x-dd31*dd23*dd44*f1x-dd31*dd43*f2x*dd14+dd31*dd23*f4x*dd14-dd31*dd24*dd13*f4x+f3x*dd43*dd14*dd21+f3x*dd11*dd23*dd44-f3x*dd11*dd43*dd24-f3x*dd13*dd44*dd21)/det
207  x3=(dd44*dd32*dd11*f2x-dd44*dd12*dd31*f2x+dd44*dd12*f3x*dd21-dd44*dd32*f1x*dd21-dd44*dd11*dd22*f3x+dd44*f1x*dd22*dd31+f4x*dd32*dd14*dd21-f4x*dd32*dd11*dd24-f4x*dd14*dd22*dd31+f4x*dd12*dd31*dd24)/det
208  x4=(-dd32*dd13*f4x*dd21-dd32*dd11*dd43*f2x+dd12*dd31*dd43*f2x-dd12*dd31*dd23*f4x-dd43*dd12*f3x*dd21+dd32*dd43*f1x*dd21+dd32*dd11*dd23*f4x+dd13*f4x*dd22*dd31+dd11*dd43*dd22*f3x-dd43*f1x*dd22*dd31)/det
209 
210 ! **** Now assign the extended boundary points ****
211  u(i1- is1,i2- is2,i3- is3,ex) = b11*x1+ g11
212  u(i1- is1,i2- is2,i3- is3,ey) = b12*x1+ g12
213  u(i1- is1,i2- is2,i3- is3,ez) = b13*x1+ g13
214  u(i1-2*is1,i2-2*is2,i3-2*is3,ex) = b21*x2+ g21
215  u(i1-2*is1,i2-2*is2,i3-2*is3,ey) = b22*x2+ g22
216  u(i1-2*is1,i2-2*is2,i3-2*is3,ez) = b23*x2+ g23
217  u(i1- js1,i2- js2,i3- js3,ex) = b31*x3+ g31
218  u(i1- js1,i2- js2,i3- js3,ey) = b32*x3+ g32
219  u(i1- js1,i2- js2,i3- js3,ez) = b33*x3+ g33
220  u(i1-2*js1,i2-2*js2,i3-2*js3,ex) = b41*x4+ g41
221  u(i1-2*js1,i2-2*js2,i3-2*js3,ey) = b42*x4+ g42
222  u(i1-2*js1,i2-2*js2,i3-2*js3,ez) = b43*x4+ g43