Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3dCurve.m @ 354:dbac99d2c318 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Removed inv(Vi) to save time
author Ylva Rydin <ylva.rydin@telia.com>
date Mon, 28 Nov 2016 08:46:28 +0100
parents 7cc3d5bd3692
children 69b078cf8072
comparison
equal deleted inserted replaced
353:fd4f1c80755d 354:dbac99d2c318
26 J, Ji 26 J, Ji
27 27
28 H % Discrete norm 28 H % Discrete norm
29 % Norms in the x, y and z directions 29 % Norms in the x, y and z directions
30 Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 30 Hxii,Hetai,Hzetai, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
31 Hxi,Heta,Hzeta
31 I_xi,I_eta,I_zeta, I_N,onesN 32 I_xi,I_eta,I_zeta, I_N,onesN
32 e_w, e_e, e_s, e_n, e_b, e_t 33 e_w, e_e, e_s, e_n, e_b, e_t
33 index_w, index_e,index_s,index_n, index_b, index_t 34 index_w, index_e,index_s,index_n, index_b, index_t
34 params %parameters for the coeficient matrice 35 params %parameters for the coeficient matrice
35 end 36 end
36 37
37 38
38 methods 39 methods
39 function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti) 40 function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti,operator)
40 xilim ={0 1}; 41 xilim ={0 1};
41 etalim = {0 1}; 42 etalim = {0 1};
42 zetalim = {0 1}; 43 zetalim = {0 1};
43 44
44 if length(m) == 1 45 if length(m) == 1
45 m = [m m m]; 46 m = [m m m];
46 end 47 end
47 m_xi = m(1); 48 m_xi = m(1);
48 m_eta = m(2); 49 m_eta = m(2);
49 m_zeta=m(3); 50 m_zeta = m(3);
50 m_tot=m_xi*m_eta*m_zeta; 51 m_tot = m_xi*m_eta*m_zeta;
51 obj.params = params; 52 obj.params = params;
52 obj.n = length(A(obj,0,0,0)); 53 obj.n = length(A(obj,0,0,0));
53 54
54 obj.m=m; 55 obj.m = m;
55 56
56 obj.order=order; 57 obj.order = order;
57 obj.onesN=ones(obj.n); 58 obj.onesN = ones(obj.n);
58 ops_xi = sbp.D2Standard(m_xi,xilim,order); 59
59 ops_eta = sbp.D2Standard(m_eta,etalim,order); 60 switch operator
60 ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); 61 case 'upwind'
62 ops_xi = sbp.D1Upwind(m_xi,xilim,order);
63 ops_eta = sbp.D1Upwind(m_eta,etalim,order);
64 ops_zeta = sbp.D1Upwind(m_zeta,zetalim,order);
65 case 'standard'
66 ops_xi = sbp.D2Standard(m_xi,xilim,order);
67 ops_eta = sbp.D2Standard(m_eta,etalim,order);
68 ops_zeta = sbp.D2Standard(m_zeta,zetalim,order);
69 otherwise
70 error('Operator not available')
71 end
61 72
62 obj.xi = ops_xi.x; 73 obj.xi = ops_xi.x;
63 obj.eta = ops_eta.x; 74 obj.eta = ops_eta.x;
64 obj.zeta = ops_zeta.x; 75 obj.zeta = ops_zeta.x;
65 76
66 obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1));%% Que pasa? 77 obj.Xi = kr(obj.xi,ones(m_eta,1),ones(m_zeta,1));%% Que pasa?
67 obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1)); 78 obj.Eta = kr(ones(m_xi,1),obj.eta,ones(m_zeta,1));
68 obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); 79 obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta);
69 80
70 obj.Eta_xi=kr(obj.eta,ones(m_xi,1)); 81 obj.Eta_xi = kr(obj.eta,ones(m_xi,1));
71 obj.Zeta_xi=kr(ones(m_eta,1),obj.zeta); 82 obj.Zeta_xi = kr(ones(m_eta,1),obj.zeta);
72 83
73 obj.Xi_eta=kr(obj.xi,ones(m_zeta,1)); 84 obj.Xi_eta = kr(obj.xi,ones(m_zeta,1));
74 obj.Zeta_eta=kr(ones(m_xi,1),obj.zeta); 85 obj.Zeta_eta = kr(ones(m_xi,1),obj.zeta);
75 86
76 obj.Xi_zeta=kr(obj.xi,ones(m_eta,1)); 87 obj.Xi_zeta = kr(obj.xi,ones(m_eta,1));
77 obj.Eta_zeta=kr(ones(m_zeta,1),obj.eta); 88 obj.Eta_zeta = kr(ones(m_zeta,1),obj.eta);
78 89
79 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); 90 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta);
80 obj.X=X; 91 obj.X = X;
81 obj.Y=Y; 92 obj.Y = Y;
82 obj.Z=Z; 93 obj.Z = Z;
83 94
84 I_n = eye(obj.n); 95 I_n = eye(obj.n);
85 I_xi = speye(m_xi); 96 I_xi = speye(m_xi);
86 obj.I_xi = I_xi; 97 obj.I_xi = I_xi;
87 I_eta = speye(m_eta); 98 I_eta = speye(m_eta);
88 obj.I_eta = I_eta; 99 obj.I_eta = I_eta;
89 I_zeta = speye(m_zeta); 100 I_zeta = speye(m_zeta);
90 obj.I_zeta = I_zeta; 101 obj.I_zeta = I_zeta;
91 102
92 103 I_N=kr(I_n,I_xi,I_eta,I_zeta);
93 O_xi=ones(m_xi,1); 104
94 O_eta=ones(m_eta,1); 105 O_xi = ones(m_xi,1);
95 O_zeta=ones(m_zeta,1); 106 O_eta = ones(m_eta,1);
96 107 O_zeta = ones(m_zeta,1);
97 D1_xi = kr(ops_xi.D1, I_eta,I_zeta); 108
98 obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); 109 obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta);
99 D1_eta = kr(I_xi, ops_eta.D1,I_zeta); 110 obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta);
100 obj.Hetai = kr(I_n, I_xi, ops_eta.HI,I_zeta);
101 D1_zeta = kr(I_xi, I_eta,ops_zeta.D1);
102 obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI); 111 obj.Hzetai = kr(I_n, I_xi,I_eta, ops_zeta.HI);
103 obj.h=[ops_xi.h ops_eta.h ops_zeta.h]; 112 obj.Hxi = ops_xi.H;
113 obj.Heta = ops_eta.H;
114 obj.Hzeta = ops_zeta.H;
115 obj.h = [ops_xi.h ops_eta.h ops_zeta.h];
116
117 switch operator
118 case 'upwind'
119 D1_xi = kr((ops_xi.Dp+ops_xi.Dm)/2, I_eta,I_zeta);
120 D1_eta = kr(I_xi, (ops_eta.Dp+ops_eta.Dm)/2,I_zeta);
121 D1_zeta = kr(I_xi, I_eta,(ops_zeta.Dp+ops_zeta.Dm)/2);
122 otherwise
123 D1_xi = kr(ops_xi.D1, I_eta,I_zeta);
124 D1_eta = kr(I_xi, ops_eta.D1,I_zeta);
125 D1_zeta = kr(I_xi, I_eta,ops_zeta.D1);
126 end
104 127
105 obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta); 128 obj.e_w = kr(I_n, ops_xi.e_l, I_eta,I_zeta);
106 obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta); 129 obj.e_e = kr(I_n, ops_xi.e_r, I_eta,I_zeta);
107 obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta); 130 obj.e_s = kr(I_n, I_xi, ops_eta.e_l,I_zeta);
108 obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta); 131 obj.e_n = kr(I_n, I_xi, ops_eta.e_r,I_zeta);
109 obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l); 132 obj.e_b = kr(I_n, I_xi, I_eta, ops_zeta.e_l);
110 obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r); 133 obj.e_t = kr(I_n, I_xi, I_eta, ops_zeta.e_r);
111 134
112 obj.A=A; 135 obj.A = A;
113 obj.B=B; 136 obj.B = B;
114 obj.C=C; 137 obj.C = C;
115 138
116 obj.X_xi=D1_xi*X; 139 obj.X_xi = D1_xi*X;
117 obj.X_eta=D1_eta*X; 140 obj.X_eta = D1_eta*X;
118 obj.X_zeta=D1_zeta*X; 141 obj.X_zeta = D1_zeta*X;
119 obj.Y_xi=D1_xi*Y; 142 obj.Y_xi = D1_xi*Y;
120 obj.Y_eta=D1_eta*Y; 143 obj.Y_eta = D1_eta*Y;
121 obj.Y_zeta=D1_zeta*Y; 144 obj.Y_zeta = D1_zeta*Y;
122 obj.Z_xi=D1_xi*Z; 145 obj.Z_xi = D1_xi*Z;
123 obj.Z_eta=D1_eta*Z; 146 obj.Z_eta = D1_eta*Z;
124 obj.Z_zeta=D1_zeta*Z; 147 obj.Z_zeta = D1_zeta*Z;
125 148
126 D1_xi=kr(I_n,D1_xi); 149 D1_xi = kr(I_n,D1_xi);
127 D1_eta=kr(I_n,D1_eta); 150 D1_eta = kr(I_n,D1_eta);
128 D1_zeta=kr(I_n,D1_zeta); 151 D1_zeta = kr(I_n,D1_zeta);
129 152
130 obj.index_w=(kr(ops_xi.e_l, O_eta,O_zeta)==1); 153 obj.index_w = (kr(ops_xi.e_l, O_eta,O_zeta)==1);
131 obj.index_e=(kr(ops_xi.e_r, O_eta,O_zeta)==1); 154 obj.index_e = (kr(ops_xi.e_r, O_eta,O_zeta)==1);
132 obj.index_s=(kr(O_xi, ops_eta.e_l,O_zeta)==1); 155 obj.index_s = (kr(O_xi, ops_eta.e_l,O_zeta)==1);
133 obj.index_n=(kr(O_xi, ops_eta.e_r,O_zeta)==1); 156 obj.index_n = (kr(O_xi, ops_eta.e_r,O_zeta)==1);
134 obj.index_b=(kr(O_xi, O_eta, ops_zeta.e_l)==1); 157 obj.index_b = (kr(O_xi, O_eta, ops_zeta.e_l)==1);
135 obj.index_t=(kr(O_xi, O_eta, ops_zeta.e_r)==1); 158 obj.index_t = (kr(O_xi, O_eta, ops_zeta.e_r)==1);
136 159
137 160
138 obj.Ahat=@transform_coefficient_matrix; 161 obj.Ahat = @transform_coefficient_matrix;
139 obj.Bhat=@transform_coefficient_matrix; 162 obj.Bhat = @transform_coefficient_matrix;
140 obj.Chat=@transform_coefficient_matrix; 163 obj.Chat = @transform_coefficient_matrix;
141 obj.E=@(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); 164 obj.E = @(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z);
142 165
143 obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta); 166 obj.Aevaluated = obj.evaluateCoefficientMatrix(obj.Ahat,obj.X, obj.Y,obj.Z, obj.X_eta,obj.X_zeta,obj.Y_eta,obj.Y_zeta,obj.Z_eta,obj.Z_zeta);
144 obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi); 167 obj.Bevaluated = obj.evaluateCoefficientMatrix(obj.Bhat,obj.X, obj.Y,obj.Z, obj.X_zeta,obj.X_xi,obj.Y_zeta,obj.Y_xi,obj.Z_zeta,obj.Z_xi);
145 obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta); 168 obj.Cevaluated = obj.evaluateCoefficientMatrix(obj.Chat,obj.X,obj.Y,obj.Z, obj.X_xi,obj.X_eta,obj.Y_xi,obj.Y_eta,obj.Z_xi,obj.Z_eta);
146 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); 169 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]);
147 170
148 obj.J=obj.X_xi.*obj.Y_eta.*obj.Z_zeta... 171 obj.J = obj.X_xi.*obj.Y_eta.*obj.Z_zeta...
149 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... 172 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta...
150 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... 173 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi...
151 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... 174 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta...
152 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... 175 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta...
153 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; 176 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi;
154 177
155 obj.Ji =kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); 178 obj.Ji = kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot));
156 179
157 obj.D=obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; 180
158 end 181 switch operator
159 182 case 'upwind'
160 function [ret]=transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) 183 alphaA = max(eig(obj.Ahat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_eta(end),obj.X_zeta(end),obj.Y_eta(end),obj.Y_zeta(end),obj.Z_eta(end),obj.Z_zeta(end))));
161 ret=obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); 184 alphaB = max(eig(obj.Bhat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_zeta(end),obj.X_xi(end),obj.Y_zeta(end),obj.Y_xi(end),obj.Z_zeta(end),obj.Z_xi(end))));
162 ret=ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2); 185 alphaC = max(eig(obj.Chat(obj,obj.X(end), obj.Y(end),obj.Z(end), obj.X_xi(end),obj.X_eta(end),obj.Y_xi(end),obj.Y_eta(end),obj.Z_xi(end),obj.Z_eta(end))));
163 ret=ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1); 186
187 Ap = (obj.Aevaluated+alphaA*I_N)/2;
188 Am = (obj.Aevaluated-alphaA*I_N)/2;
189 Bp = (obj.Bevaluated+alphaB*I_N)/2;
190 Bm = (obj.Bevaluated-alphaB*I_N)/2;
191 Cp = (obj.Cevaluated+alphaC*I_N)/2;
192 Cm = (obj.Cevaluated-alphaC*I_N)/2;
193
194 Dpxi = kr(I_n, ops_xi.Dp, I_eta,I_zeta);
195 Dmxi = kr(I_n, ops_xi.Dm, I_eta,I_zeta);
196 Dpeta = kr(I_n, I_xi, ops_eta.Dp,I_zeta);
197 Dmeta = kr(I_n, I_xi, ops_eta.Dm,I_zeta);
198 Dpzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dp);
199 Dmzeta = kr(I_n, I_xi, I_eta,ops_zeta.Dm);
200
201 obj.D = obj.Ji*(-Am*Dpxi-Ap*Dmxi-Bm*Dpeta-Bp*Dmeta-Cm*Dpzeta-Cp*Dmzeta)-obj.Eevaluated;
202 otherwise
203 obj.D = obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated;
204 end
205
206
207 end
208
209 function [ret] = transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)
210 ret = obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2);
211 ret = ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2);
212 ret = ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1);
164 end 213 end
165 214
166 215
167 % Closure functions return the opertors applied to the own doamin to close the boundary 216 % Closure functions return the opertors applied to the own doamin to close the boundary
168 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 217 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
169 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 218 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
170 % type is a string specifying the type of boundary condition if there are several. 219 % type is a string specifying the type of boundary condition if there are several.
171 % data is a function returning the data that should be applied at the boundary. 220 % data is a function returning the data that should be applied at the boundary.
172 function [closure, penalty] = boundary_condition(obj,boundary,type,L) 221 function [closure, penalty] = boundary_condition(obj,boundary,type,L)
173 default_arg('type','char'); 222 default_arg('type','char');
174 BM=boundary_matrices(obj,boundary); 223 BM = boundary_matrices(obj,boundary);
175 224
176 switch type 225 switch type
177 case{'c','char'} 226 case{'c','char'}
178 [closure,penalty]=boundary_condition_char(obj,BM); 227 [closure,penalty] = boundary_condition_char(obj,BM);
179 case{'general'} 228 case{'general'}
180 [closure,penalty]=boundary_condition_general(obj,BM,boundary,L); 229 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L);
181 otherwise 230 otherwise
182 error('No such boundary condition') 231 error('No such boundary condition')
183 end 232 end
184 end 233 end
185 234
190 function N = size(obj) 239 function N = size(obj)
191 N = obj.m; 240 N = obj.m;
192 end 241 end
193 242
194 function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2) 243 function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2)
195 params=obj.params; 244 params = obj.params;
196 side=max(length(X),length(Y)); 245 side = max(length(X),length(Y));
197 if isa(mat,'function_handle') 246 if isa(mat,'function_handle')
198 [rows,cols]=size(mat(obj,0,0,0,0,0,0,0,0,0)); 247 [rows,cols] = size(mat(obj,0,0,0,0,0,0,0,0,0));
199 x_1=kr(obj.onesN,x_1); 248 x_1 = kr(obj.onesN,x_1);
200 x_2=kr(obj.onesN,x_2); 249 x_2 = kr(obj.onesN,x_2);
201 y_1=kr(obj.onesN,y_1); 250 y_1 = kr(obj.onesN,y_1);
202 y_2=kr(obj.onesN,y_2); 251 y_2 = kr(obj.onesN,y_2);
203 z_1=kr(obj.onesN,z_1); 252 z_1 = kr(obj.onesN,z_1);
204 z_2=kr(obj.onesN,z_2); 253 z_2 = kr(obj.onesN,z_2);
205 matVec=mat(obj,X',Y',Z',x_1',x_2',y_1',y_2',z_1',z_2'); 254 matVec = mat(obj,X',Y',Z',x_1',x_2',y_1',y_2',z_1',z_2');
206 matVec=sparse(matVec); 255 matVec = sparse(matVec);
207 else 256 else
208 matVec=mat; 257 matVec = mat;
209 [rows,cols]=size(matVec); 258 [rows,cols] = size(matVec);
210 side=max(length(X),length(Y)); 259 side = max(length(X),length(Y));
211 cols=cols/side; 260 cols = cols/side;
212 end 261 end
213 ret=cell(rows,cols); 262 ret = cell(rows,cols);
214 263
215 264
216 for ii=1:rows 265 for ii=1:rows
217 for jj=1:cols 266 for jj=1:cols
218 ret{ii,jj}=diag(matVec(ii,(jj-1)*side+1:jj*side)); 267 ret{ii,jj} = diag(matVec(ii,(jj-1)*side+1:jj*side));
219 end 268 end
220 end 269 end
221 270
222 ret=cell2mat(ret); 271 ret = cell2mat(ret);
223 end 272 end
224 273
225 274
226 function [BM]=boundary_matrices(obj,boundary) 275 function [BM] = boundary_matrices(obj,boundary)
227 params=obj.params; 276 params = obj.params;
228 BM.boundary=boundary; 277 BM.boundary = boundary;
229 switch boundary 278 switch boundary
230 case {'w','W','west'} 279 case {'w','W','west'}
231 BM.e_=obj.e_w; 280 BM.e_ = obj.e_w;
232 mat=obj.Ahat; 281 mat = obj.Ahat;
233 BM.boundpos='l'; 282 BM.boundpos = 'l';
234 BM.Hi=obj.Hxii; 283 BM.Hi = obj.Hxii;
235 BM.index=obj.index_w; 284 BM.index = obj.index_w;
236 BM.x_1=obj.X_eta(BM.index); 285 BM.x_1 = obj.X_eta(BM.index);
237 BM.x_2=obj.X_zeta(BM.index); 286 BM.x_2 = obj.X_zeta(BM.index);
238 BM.y_1=obj.Y_eta(BM.index); 287 BM.y_1 = obj.Y_eta(BM.index);
239 BM.y_2=obj.Y_zeta(BM.index); 288 BM.y_2 = obj.Y_zeta(BM.index);
240 BM.z_1=obj.Z_eta(BM.index); 289 BM.z_1 = obj.Z_eta(BM.index);
241 BM.z_2=obj.Z_zeta(BM.index); 290 BM.z_2 = obj.Z_zeta(BM.index);
242 case {'e','E','east'} 291 case {'e','E','east'}
243 BM.e_=obj.e_e; 292 BM.e_ = obj.e_e;
244 mat=obj.Ahat; 293 mat = obj.Ahat;
245 BM.boundpos='r'; 294 BM.boundpos = 'r';
246 BM.Hi=obj.Hxii; 295 BM.Hi = obj.Hxii;
247 BM.index=obj.index_e; 296 BM.index = obj.index_e;
248 BM.x_1=obj.X_eta(BM.index); 297 BM.x_1 = obj.X_eta(BM.index);
249 BM.x_2=obj.X_zeta(BM.index); 298 BM.x_2 = obj.X_zeta(BM.index);
250 BM.y_1=obj.Y_eta(BM.index); 299 BM.y_1 = obj.Y_eta(BM.index);
251 BM.y_2=obj.Y_zeta(BM.index); 300 BM.y_2 = obj.Y_zeta(BM.index);
252 BM.z_1=obj.Z_eta(BM.index); 301 BM.z_1 = obj.Z_eta(BM.index);
253 BM.z_2=obj.Z_zeta(BM.index); 302 BM.z_2 = obj.Z_zeta(BM.index);
254 case {'s','S','south'} 303 case {'s','S','south'}
255 BM.e_=obj.e_s; 304 BM.e_ = obj.e_s;
256 mat=obj.Bhat; 305 mat = obj.Bhat;
257 BM.boundpos='l'; 306 BM.boundpos = 'l';
258 BM.Hi=obj.Hetai; 307 BM.Hi = obj.Hetai;
259 BM.index=obj.index_s; 308 BM.index = obj.index_s;
260 BM.x_1=obj.X_zeta(BM.index); 309 BM.x_1 = obj.X_zeta(BM.index);
261 BM.x_2=obj.X_xi(BM.index); 310 BM.x_2 = obj.X_xi(BM.index);
262 BM.y_1=obj.Y_zeta(BM.index); 311 BM.y_1 = obj.Y_zeta(BM.index);
263 BM.y_2=obj.Y_xi(BM.index); 312 BM.y_2 = obj.Y_xi(BM.index);
264 BM.z_1=obj.Z_zeta(BM.index); 313 BM.z_1 = obj.Z_zeta(BM.index);
265 BM.z_2=obj.Z_xi(BM.index); 314 BM.z_2 = obj.Z_xi(BM.index);
266 case {'n','N','north'} 315 case {'n','N','north'}
267 BM.e_=obj.e_n; 316 BM.e_ = obj.e_n;
268 mat=obj.Bhat; 317 mat = obj.Bhat;
269 BM.boundpos='r'; 318 BM.boundpos = 'r';
270 BM.Hi=obj.Hetai; 319 BM.Hi = obj.Hetai;
271 BM.index=obj.index_n; 320 BM.index = obj.index_n;
272 BM.x_1=obj.X_zeta(BM.index); 321 BM.x_1 = obj.X_zeta(BM.index);
273 BM.x_2=obj.X_xi(BM.index); 322 BM.x_2 = obj.X_xi(BM.index);
274 BM.y_1=obj.Y_zeta(BM.index); 323 BM.y_1 = obj.Y_zeta(BM.index);
275 BM.y_2=obj.Y_xi(BM.index); 324 BM.y_2 = obj.Y_xi(BM.index);
276 BM.z_1=obj.Z_zeta(BM.index); 325 BM.z_1 = obj.Z_zeta(BM.index);
277 BM.z_2=obj.Z_xi(BM.index); 326 BM.z_2 = obj.Z_xi(BM.index);
278 case{'b','B','Bottom'} 327 case{'b','B','Bottom'}
279 BM.e_=obj.e_b; 328 BM.e_ = obj.e_b;
280 mat=obj.Chat; 329 mat = obj.Chat;
281 BM.boundpos='l'; 330 BM.boundpos = 'l';
282 BM.Hi=obj.Hzetai; 331 BM.Hi = obj.Hzetai;
283 BM.index=obj.index_b; 332 BM.index = obj.index_b;
284 BM.x_1=obj.X_xi(BM.index); 333 BM.x_1 = obj.X_xi(BM.index);
285 BM.x_2=obj.X_eta(BM.index); 334 BM.x_2 = obj.X_eta(BM.index);
286 BM.y_1=obj.Y_xi(BM.index); 335 BM.y_1 = obj.Y_xi(BM.index);
287 BM.y_2=obj.Y_eta(BM.index); 336 BM.y_2 = obj.Y_eta(BM.index);
288 BM.z_1=obj.Z_xi(BM.index); 337 BM.z_1 = obj.Z_xi(BM.index);
289 BM.z_2=obj.Z_eta(BM.index); 338 BM.z_2 = obj.Z_eta(BM.index);
290 case{'t','T','Top'} 339 case{'t','T','Top'}
291 BM.e_=obj.e_t; 340 BM.e_ = obj.e_t;
292 mat=obj.Chat; 341 mat = obj.Chat;
293 BM.boundpos='r'; 342 BM.boundpos = 'r';
294 BM.Hi=obj.Hzetai; 343 BM.Hi = obj.Hzetai;
295 BM.index=obj.index_t; 344 BM.index = obj.index_t;
296 BM.x_1=obj.X_xi(BM.index); 345 BM.x_1 = obj.X_xi(BM.index);
297 BM.x_2=obj.X_eta(BM.index); 346 BM.x_2 = obj.X_eta(BM.index);
298 BM.y_1=obj.Y_xi(BM.index); 347 BM.y_1 = obj.Y_xi(BM.index);
299 BM.y_2=obj.Y_eta(BM.index); 348 BM.y_2 = obj.Y_eta(BM.index);
300 BM.z_1=obj.Z_xi(BM.index); 349 BM.z_1 = obj.Z_xi(BM.index);
301 BM.z_2=obj.Z_eta(BM.index); 350 BM.z_2 = obj.Z_eta(BM.index);
302 end 351 end
303 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(BM.index),obj.Y(BM.index),obj.Z(BM.index),... 352 [BM.V,BM.Vi,BM.D,signVec] = obj.matrixDiag(mat,obj.X(BM.index),obj.Y(BM.index),obj.Z(BM.index),...
304 BM.x_1,BM.x_2,BM.y_1,BM.y_2,BM.z_1,BM.z_2); 353 BM.x_1,BM.x_2,BM.y_1,BM.y_2,BM.z_1,BM.z_2);
305 BM.side=sum(BM.index); 354 BM.side = sum(BM.index);
306 BM.pos=signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); 355 BM.pos = signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3);
307 end 356 end
308 357
309 358
310 function [closure, penalty]=boundary_condition_char(obj,BM) 359 function [closure, penalty]=boundary_condition_char(obj,BM)
311 side = BM.side; 360 side = BM.side;
312 pos = BM.pos; 361 pos = BM.pos;
313 neg = BM.neg; 362 neg = BM.neg;
314 zeroval=BM.zeroval; 363 zeroval = BM.zeroval;
315 V = BM.V; 364 V = BM.V;
316 Vi = BM.Vi; 365 Vi = BM.Vi;
317 Hi=BM.Hi; 366 Hi = BM.Hi;
318 D=BM.D; 367 D = BM.D;
319 e_=BM.e_; 368 e_ = BM.e_;
320 369
321 switch BM.boundpos 370 switch BM.boundpos
322 case {'l'} 371 case {'l'}
323 tau=sparse(obj.n*side,pos); 372 tau = sparse(obj.n*side,pos);
324 Vi_plus=Vi(1:pos,:); 373 Vi_plus = Vi(1:pos,:);
325 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 374 tau(1:pos,:) = -abs(D(1:pos,1:pos));
326 closure=Hi*e_*V*tau*Vi_plus*e_'; 375 closure = Hi*e_*V*tau*Vi_plus*e_';
327 penalty=-Hi*e_*V*tau*Vi_plus; 376 penalty = -Hi*e_*V*tau*Vi_plus;
328 case {'r'} 377 case {'r'}
329 tau=sparse(obj.n*side,neg); 378 tau = sparse(obj.n*side,neg);
330 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 379 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
331 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 380 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
332 closure=Hi*e_*V*tau*Vi_minus*e_'; 381 closure = Hi*e_*V*tau*Vi_minus*e_';
333 penalty=-Hi*e_*V*tau*Vi_minus; 382 penalty = -Hi*e_*V*tau*Vi_minus;
334 end 383 end
335 end 384 end
336 385
337 386
338 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) 387 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L)
339 side = BM.side; 388 side = BM.side;
340 pos = BM.pos; 389 pos = BM.pos;
341 neg = BM.neg; 390 neg = BM.neg;
342 zeroval=BM.zeroval; 391 zeroval = BM.zeroval;
343 V = BM.V; 392 V = BM.V;
344 Vi = BM.Vi; 393 Vi = BM.Vi;
345 Hi=BM.Hi; 394 Hi = BM.Hi;
346 D=BM.D; 395 D = BM.D;
347 e_=BM.e_; 396 e_ = BM.e_;
348 index=BM.index; 397 index = BM.index;
349 398
350 switch BM.boundary 399 switch BM.boundary
351 case{'b','B','bottom'} 400 case{'b','B','bottom'}
352 Ji_vec=diag(obj.Ji); 401 Ji_vec = diag(obj.Ji);
353 Ji=diag(Ji_vec(index)); 402 Ji = diag(Ji_vec(index));
354 Zeta_x=Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index)); 403 Zeta_x = Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index));
355 Zeta_y=Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index)); 404 Zeta_y = Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index));
356 Zeta_z=Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index)); 405 Zeta_z = Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index));
357 406
358 L=obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]); 407 L = obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]);
359 end 408 end
360 409
361 switch BM.boundpos 410 switch BM.boundpos
362 case {'l'} 411 case {'l'}
363 tau=sparse(obj.n*side,pos); 412 tau = sparse(obj.n*side,pos);
364 Vi_plus=Vi(1:pos,:); 413 Vi_plus = Vi(1:pos,:);
365 Vi_minus=Vi(pos+zeroval+1:obj.n*side,:); 414 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
366 V_plus=V(:,1:pos); 415 V_plus = V(:,1:pos);
367 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 416 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
368 417
369 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 418 tau(1:pos,:) = -abs(D(1:pos,1:pos));
370 R=-inv(L*V_plus)*(L*V_minus); 419 R = -inv(L*V_plus)*(L*V_minus);
371 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 420 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
372 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 421 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
373 case {'r'} 422 case {'r'}
374 tau=sparse(obj.n*side,neg); 423 tau = sparse(obj.n*side,neg);
375 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 424 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
376 Vi_plus=Vi(1:pos,:); 425 Vi_plus = Vi(1:pos,:);
377 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 426 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
378 427
379 V_plus=V(:,1:pos); 428 V_plus = V(:,1:pos);
380 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 429 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
381 R=-inv(L*V_minus)*(L*V_plus); 430 R = -inv(L*V_minus)*(L*V_plus);
382 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 431 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
383 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 432 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
384 end 433 end
385 end 434 end
386 435
387 436
388 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) 437 function [V,Vi, D,signVec] = matrixDiag(obj,mat,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)
389 params=obj.params; 438 params = obj.params;
390 eps=10^(-10); 439 eps = 10^(-10);
391 if(sum(abs(x_1))>eps) 440 if(sum(abs(x_1))>eps)
392 syms x_1s 441 syms x_1s
393 else 442 else
394 x_1s=0; 443 x_1s = 0;
395 end 444 end
396 445
397 if(sum(abs(x_2))>eps) 446 if(sum(abs(x_2))>eps)
398 syms x_2s; 447 syms x_2s;
399 else 448 else
400 x_2s=0; 449 x_2s = 0;
401 end 450 end
402 451
403 452
404 if(sum(abs(y_1))>eps) 453 if(sum(abs(y_1))>eps)
405 syms y_1s 454 syms y_1s
406 else 455 else
407 y_1s=0; 456 y_1s = 0;
408 end 457 end
409 458
410 if(sum(abs(y_2))>eps) 459 if(sum(abs(y_2))>eps)
411 syms y_2s; 460 syms y_2s;
412 else 461 else
413 y_2s=0; 462 y_2s = 0;
414 end 463 end
415 464
416 if(sum(abs(z_1))>eps) 465 if(sum(abs(z_1))>eps)
417 syms z_1s 466 syms z_1s
418 else 467 else
419 z_1s=0; 468 z_1s = 0;
420 end 469 end
421 470
422 if(sum(abs(z_2))>eps) 471 if(sum(abs(z_2))>eps)
423 syms z_2s; 472 syms z_2s;
424 else 473 else
425 z_2s=0; 474 z_2s = 0;
426 end 475 end
427 476
428 syms xs ys zs 477 syms xs ys zs
429 [V, D]=eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); 478 [V, D] = eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s));
430 Vi=inv(V); 479 Vi = inv(V);
431 % syms x_1s x_2s y_1s y_2s z_1s z_2s 480 % syms x_1s x_2s y_1s y_2s z_1s z_2s
432 xs=x; 481 xs = x;
433 ys=y; 482 ys = y;
434 zs=z; 483 zs = z;
435 x_1s=x_1; 484 x_1s = x_1;
436 x_2s=x_2; 485 x_2s = x_2;
437 y_1s=y_1; 486 y_1s = y_1;
438 y_2s=y_2; 487 y_2s = y_2;
439 z_1s=z_1; 488 z_1s = z_1;
440 z_2s=z_2; 489 z_2s = z_2;
441 490
442 side=max(length(x),length(y)); 491 side = max(length(x),length(y));
443 Dret=zeros(obj.n,side*obj.n); 492 Dret = zeros(obj.n,side*obj.n);
444 Vret=zeros(obj.n,side*obj.n); 493 Vret = zeros(obj.n,side*obj.n);
445 Viret=zeros(obj.n,side*obj.n); 494 Viret = zeros(obj.n,side*obj.n);
446 495
447 for ii=1:obj.n 496 for ii=1:obj.n
448 for jj=1:obj.n 497 for jj=1:obj.n
449 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 498 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
450 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 499 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
451 Viret(jj,(ii-1)*side+1:side*ii)=eval(Vi(jj,ii)); 500 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
452 end 501 end
453 end 502 end
454 503
455 D=sparse(Dret); 504 D = sparse(Dret);
456 V=sparse(Vret); 505 V = sparse(Vret);
457 Vi=sparse(Viret); 506 Vi = sparse(Viret);
458 V=obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); 507 V = obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2);
459 D=obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); 508 D = obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2);
460 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); 509 Vi = obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2);
461 DD=diag(D); 510 DD = diag(D);
462 511
463 poseig=(DD>0); 512 poseig = (DD>0);
464 zeroeig=(DD==0); 513 zeroeig = (DD==0);
465 negeig=(DD<0); 514 negeig = (DD<0);
466 515
467 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); 516 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
468 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 517 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
469 Vi=[Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 518 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
470 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; 519 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
471 end 520 end
472 end 521 end
473 end 522 end