Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3dCurve.m @ 395:359861563866 feature/beams

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