Mercurial > repos > public > sbplib
comparison +scheme/Hypsyst3dCurve.m @ 350:5d5652fe826a feature/hypsyst
A commit before I try resolving the performance issues
author | Ylva Rydin <ylva.rydin@telia.com> |
---|---|
date | Wed, 02 Nov 2016 00:02:01 +0100 |
parents | |
children | 7cc3d5bd3692 |
comparison
equal
deleted
inserted
replaced
349:cd6a29ab3746 | 350:5d5652fe826a |
---|---|
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 | |
13 | |
14 X_xi, X_eta, X_zeta,Y_xi,Y_eta,Y_zeta,Z_xi,Z_eta,Z_zeta | |
15 Aev | |
16 | |
17 metric_terms | |
18 | |
19 order % Order accuracy for the approximation | |
20 | |
21 D % non-stabalized scheme operator | |
22 Aevaluated, Bevaluated, Cevaluated, Eevaluated | |
23 Ahat, Bhat, Chat, E | |
24 A,B,C | |
25 | |
26 J, Ji | |
27 | |
28 H % Discrete norm | |
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. | |
31 I_xi,I_eta,I_zeta, I_N,onesN | |
32 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 params %parameters for the coeficient matrice | |
35 end | |
36 | |
37 | |
38 methods | |
39 function obj = Hypsyst3dCurve(m, order, A, B,C, E, params,ti) | |
40 xilim ={0 1}; | |
41 etalim = {0 1}; | |
42 zetalim = {0 1}; | |
43 | |
44 if length(m) == 1 | |
45 m = [m m m]; | |
46 end | |
47 m_xi = m(1); | |
48 m_eta = m(2); | |
49 m_zeta=m(3); | |
50 m_tot=m_xi*m_eta*m_zeta; | |
51 obj.params = params; | |
52 obj.n = length(A(obj,0,0,0)); | |
53 | |
54 obj.m=m; | |
55 | |
56 obj.order=order; | |
57 obj.onesN=ones(obj.n); | |
58 ops_xi = sbp.D2Standard(m_xi,xilim,order); | |
59 ops_eta = sbp.D2Standard(m_eta,etalim,order); | |
60 ops_zeta = sbp.D2Standard(m_zeta,zetalim,order); | |
61 | |
62 obj.xi = ops_xi.x; | |
63 obj.eta = ops_eta.x; | |
64 obj.zeta = ops_zeta.x; | |
65 | |
66 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)); | |
68 obj.Zeta = kr(ones(m_xi,1),ones(m_eta,1),obj.zeta); | |
69 | |
70 obj.Eta_xi=kr(obj.eta,ones(m_xi,1)); | |
71 obj.Zeta_xi=kr(ones(m_eta,1),obj.zeta); | |
72 | |
73 obj.Xi_eta=kr(obj.xi,ones(m_zeta,1)); | |
74 obj.Zeta_eta=kr(ones(m_xi,1),obj.zeta); | |
75 | |
76 obj.Xi_zeta=kr(obj.xi,ones(m_eta,1)); | |
77 obj.Eta_zeta=kr(ones(m_zeta,1),obj.eta); | |
78 | |
79 [X,Y,Z] = ti.map(obj.Xi,obj.Eta,obj.Zeta); | |
80 obj.X=X; | |
81 obj.Y=Y; | |
82 obj.Z=Z; | |
83 | |
84 I_n = eye(obj.n); | |
85 I_xi = speye(m_xi); | |
86 obj.I_xi = I_xi; | |
87 I_eta = speye(m_eta); | |
88 obj.I_eta = I_eta; | |
89 I_zeta = speye(m_zeta); | |
90 obj.I_zeta = I_zeta; | |
91 | |
92 | |
93 O_xi=ones(m_xi,1); | |
94 O_eta=ones(m_eta,1); | |
95 O_zeta=ones(m_zeta,1); | |
96 | |
97 D1_xi = kr(ops_xi.D1, I_eta,I_zeta); | |
98 obj.Hxii = kr(I_n, ops_xi.HI, I_eta,I_zeta); | |
99 D1_eta = kr(I_xi, ops_eta.D1,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); | |
103 obj.h=[ops_xi.h ops_eta.h ops_zeta.h]; | |
104 | |
105 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); | |
107 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); | |
109 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); | |
111 | |
112 obj.A=A; | |
113 obj.B=B; | |
114 obj.C=C; | |
115 | |
116 obj.X_xi=D1_xi*X; | |
117 obj.X_eta=D1_eta*X; | |
118 obj.X_zeta=D1_zeta*X; | |
119 obj.Y_xi=D1_xi*Y; | |
120 obj.Y_eta=D1_eta*Y; | |
121 obj.Y_zeta=D1_zeta*Y; | |
122 obj.Z_xi=D1_xi*Z; | |
123 obj.Z_eta=D1_eta*Z; | |
124 obj.Z_zeta=D1_zeta*Z; | |
125 | |
126 D1_xi=kr(I_n,D1_xi); | |
127 D1_eta=kr(I_n,D1_eta); | |
128 D1_zeta=kr(I_n,D1_zeta); | |
129 | |
130 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); | |
132 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); | |
134 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); | |
136 | |
137 | |
138 obj.Ahat=@transform_coefficient_matrix; | |
139 obj.Bhat=@transform_coefficient_matrix; | |
140 obj.Chat=@transform_coefficient_matrix; | |
141 obj.E=@(obj,x,y,z,~,~,~,~,~,~)E(obj,x,y,z); | |
142 | |
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); | |
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); | |
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); | |
146 obj.Eevaluated = obj.evaluateCoefficientMatrix(obj.E, obj.X, obj.Y,obj.Z,[],[],[],[],[],[]); | |
147 | |
148 obj.J=obj.X_xi.*obj.Y_eta.*obj.Z_zeta... | |
149 +obj.X_zeta.*obj.Y_xi.*obj.Z_eta... | |
150 +obj.X_eta.*obj.Y_zeta.*obj.Z_xi... | |
151 -obj.X_xi.*obj.Y_zeta.*obj.Z_eta... | |
152 -obj.X_eta.*obj.Y_xi.*obj.Z_zeta... | |
153 -obj.X_zeta.*obj.Y_eta.*obj.Z_xi; | |
154 | |
155 obj.Ji =kr(I_n,spdiags(1./obj.J,0,m_tot,m_tot)); | |
156 | |
157 obj.D=obj.Ji*(-obj.Aevaluated*D1_xi-obj.Bevaluated*D1_eta -obj.Cevaluated*D1_zeta)-obj.Eevaluated; | |
158 end | |
159 | |
160 function [ret]=transform_coefficient_matrix(obj,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) | |
161 ret=obj.A(obj,x,y,z).*(y_1.*z_2-z_1.*y_2); | |
162 ret=ret+obj.B(obj,x,y,z).*(x_2.*z_1-x_1.*z_2); | |
163 ret=ret+obj.C(obj,x,y,z).*(x_1.*y_2-x_2.*y_1); | |
164 end | |
165 | |
166 | |
167 % 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. | |
169 % 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. | |
171 % data is a function returning the data that should be applied at the boundary. | |
172 function [closure, penalty] = boundary_condition(obj,boundary,type,L) | |
173 default_arg('type','char'); | |
174 BM=boundary_matrices(obj,boundary); | |
175 | |
176 switch type | |
177 case{'c','char'} | |
178 [closure,penalty]=boundary_condition_char(obj,BM); | |
179 case{'general'} | |
180 [closure,penalty]=boundary_condition_general(obj,BM,boundary,L); | |
181 otherwise | |
182 error('No such boundary condition') | |
183 end | |
184 end | |
185 | |
186 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) | |
187 error('An interface function does not exist yet'); | |
188 end | |
189 | |
190 function N = size(obj) | |
191 N = obj.m; | |
192 end | |
193 | |
194 function [ret] = evaluateCoefficientMatrix(obj,mat, X, Y, Z , x_1 , x_2 , y_1 , y_2 , z_1 , z_2) | |
195 params=obj.params; | |
196 side=max(length(X),length(Y)); | |
197 if isa(mat,'function_handle') | |
198 [rows,cols]=size(mat(obj,0,0,0,0,0,0,0,0,0)); | |
199 x_1=kr(obj.onesN,x_1); | |
200 x_2=kr(obj.onesN,x_2); | |
201 y_1=kr(obj.onesN,y_1); | |
202 y_2=kr(obj.onesN,y_2); | |
203 z_1=kr(obj.onesN,z_1); | |
204 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'); | |
206 matVec=sparse(matVec); | |
207 else | |
208 matVec=mat; | |
209 [rows,cols]=size(matVec); | |
210 side=max(length(X),length(Y)); | |
211 cols=cols/side; | |
212 end | |
213 ret=kron(ones(rows,cols),speye(side)); | |
214 | |
215 for ii=1:rows | |
216 for jj=1:cols | |
217 ret((ii-1)*side+1:ii*side,(jj-1)*side+1:jj*side)=diag(matVec(ii,(jj-1)*side+1:jj*side)); | |
218 end | |
219 end | |
220 end | |
221 | |
222 | |
223 function [BM]=boundary_matrices(obj,boundary) | |
224 params=obj.params; | |
225 BM.boundary=boundary; | |
226 switch boundary | |
227 case {'w','W','west'} | |
228 BM.e_=obj.e_w; | |
229 mat=obj.Ahat; | |
230 BM.boundpos='l'; | |
231 BM.Hi=obj.Hxii; | |
232 BM.index=obj.index_w; | |
233 BM.x_1=obj.X_eta(BM.index); | |
234 BM.x_2=obj.X_zeta(BM.index); | |
235 BM.y_1=obj.Y_eta(BM.index); | |
236 BM.y_2=obj.Y_zeta(BM.index); | |
237 BM.z_1=obj.Z_eta(BM.index); | |
238 BM.z_2=obj.Z_zeta(BM.index); | |
239 case {'e','E','east'} | |
240 BM.e_=obj.e_e; | |
241 mat=obj.Ahat; | |
242 BM.boundpos='r'; | |
243 BM.Hi=obj.Hxii; | |
244 BM.index=obj.index_e; | |
245 BM.x_1=obj.X_eta(BM.index); | |
246 BM.x_2=obj.X_zeta(BM.index); | |
247 BM.y_1=obj.Y_eta(BM.index); | |
248 BM.y_2=obj.Y_zeta(BM.index); | |
249 BM.z_1=obj.Z_eta(BM.index); | |
250 BM.z_2=obj.Z_zeta(BM.index); | |
251 case {'s','S','south'} | |
252 BM.e_=obj.e_s; | |
253 mat=obj.Bhat; | |
254 BM.boundpos='l'; | |
255 BM.Hi=obj.Hetai; | |
256 BM.index=obj.index_s; | |
257 BM.x_1=obj.X_zeta(BM.index); | |
258 BM.x_2=obj.X_xi(BM.index); | |
259 BM.y_1=obj.Y_zeta(BM.index); | |
260 BM.y_2=obj.Y_xi(BM.index); | |
261 BM.z_1=obj.Z_zeta(BM.index); | |
262 BM.z_2=obj.Z_xi(BM.index); | |
263 case {'n','N','north'} | |
264 BM.e_=obj.e_n; | |
265 mat=obj.Bhat; | |
266 BM.boundpos='r'; | |
267 BM.Hi=obj.Hetai; | |
268 BM.index=obj.index_n; | |
269 BM.x_1=obj.X_zeta(BM.index); | |
270 BM.x_2=obj.X_xi(BM.index); | |
271 BM.y_1=obj.Y_zeta(BM.index); | |
272 BM.y_2=obj.Y_xi(BM.index); | |
273 BM.z_1=obj.Z_zeta(BM.index); | |
274 BM.z_2=obj.Z_xi(BM.index); | |
275 case{'b','B','Bottom'} | |
276 BM.e_=obj.e_b; | |
277 mat=obj.Chat; | |
278 BM.boundpos='l'; | |
279 BM.Hi=obj.Hzetai; | |
280 BM.index=obj.index_b; | |
281 BM.x_1=obj.X_xi(BM.index); | |
282 BM.x_2=obj.X_eta(BM.index); | |
283 BM.y_1=obj.Y_xi(BM.index); | |
284 BM.y_2=obj.Y_eta(BM.index); | |
285 BM.z_1=obj.Z_xi(BM.index); | |
286 BM.z_2=obj.Z_eta(BM.index); | |
287 case{'t','T','Top'} | |
288 BM.e_=obj.e_t; | |
289 mat=obj.Chat; | |
290 BM.boundpos='r'; | |
291 BM.Hi=obj.Hzetai; | |
292 BM.index=obj.index_t; | |
293 BM.x_1=obj.X_xi(BM.index); | |
294 BM.x_2=obj.X_eta(BM.index); | |
295 BM.y_1=obj.Y_xi(BM.index); | |
296 BM.y_2=obj.Y_eta(BM.index); | |
297 BM.z_1=obj.Z_xi(BM.index); | |
298 BM.z_2=obj.Z_eta(BM.index); | |
299 end | |
300 [BM.V,BM.Vi,BM.D,signVec]=obj.matrixDiag(mat,obj.X(BM.index),obj.Y(BM.index),obj.Z(BM.index),... | |
301 BM.x_1,BM.x_2,BM.y_1,BM.y_2,BM.z_1,BM.z_2); | |
302 BM.side=sum(BM.index); | |
303 BM.pos=signVec(1); BM.zeroval=signVec(2); BM.neg=signVec(3); | |
304 end | |
305 | |
306 | |
307 function [closure, penalty]=boundary_condition_char(obj,BM) | |
308 side = BM.side; | |
309 pos = BM.pos; | |
310 neg = BM.neg; | |
311 zeroval=BM.zeroval; | |
312 V = BM.V; | |
313 Vi = BM.Vi; | |
314 Hi=BM.Hi; | |
315 D=BM.D; | |
316 e_=BM.e_; | |
317 | |
318 switch BM.boundpos | |
319 case {'l'} | |
320 tau=sparse(obj.n*side,pos); | |
321 Vi_plus=Vi(1:pos,:); | |
322 tau(1:pos,:)=-abs(D(1:pos,1:pos)); | |
323 closure=Hi*e_*V*tau*Vi_plus*e_'; | |
324 penalty=-Hi*e_*V*tau*Vi_plus; | |
325 case {'r'} | |
326 tau=sparse(obj.n*side,neg); | |
327 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); | |
328 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); | |
329 closure=Hi*e_*V*tau*Vi_minus*e_'; | |
330 penalty=-Hi*e_*V*tau*Vi_minus; | |
331 end | |
332 end | |
333 | |
334 | |
335 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) | |
336 side = BM.side; | |
337 pos = BM.pos; | |
338 neg = BM.neg; | |
339 zeroval=BM.zeroval; | |
340 V = BM.V; | |
341 Vi = BM.Vi; | |
342 Hi=BM.Hi; | |
343 D=BM.D; | |
344 e_=BM.e_; | |
345 index=BM.index; | |
346 | |
347 switch BM.boundary | |
348 case{'b','B','bottom'} | |
349 Ji_vec=diag(obj.Ji); | |
350 Ji=diag(Ji_vec(index)); | |
351 Zeta_x=Ji*(obj.Y_xi(index).*obj.Z_eta(index)-obj.Z_xi(index).*obj.Y_eta(index)); | |
352 Zeta_y=Ji*(obj.X_eta(index).*obj.Z_xi(index)-obj.X_xi(index).*obj.Z_eta(index)); | |
353 Zeta_z=Ji*(obj.X_xi(index).*obj.Y_eta(index)-obj.Y_xi(index).*obj.X_eta(index)); | |
354 | |
355 L=obj.evaluateCoefficientMatrix(L,Zeta_x,Zeta_y,Zeta_z,[],[],[],[],[],[]); | |
356 end | |
357 | |
358 switch BM.boundpos | |
359 case {'l'} | |
360 tau=sparse(obj.n*side,pos); | |
361 Vi_plus=Vi(1:pos,:); | |
362 Vi_minus=Vi(pos+zeroval+1:obj.n*side,:); | |
363 V_plus=V(:,1:pos); | |
364 V_minus=V(:,(pos+zeroval)+1:obj.n*side); | |
365 | |
366 tau(1:pos,:)=-abs(D(1:pos,1:pos)); | |
367 R=-inv(L*V_plus)*(L*V_minus); | |
368 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; | |
369 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; | |
370 case {'r'} | |
371 tau=sparse(obj.n*side,neg); | |
372 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); | |
373 Vi_plus=Vi(1:pos,:); | |
374 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); | |
375 | |
376 V_plus=V(:,1:pos); | |
377 V_minus=V(:,(pos+zeroval)+1:obj.n*side); | |
378 R=-inv(L*V_minus)*(L*V_plus); | |
379 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; | |
380 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; | |
381 end | |
382 end | |
383 | |
384 | |
385 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2) | |
386 params=obj.params; | |
387 eps=10^(-10); | |
388 if(sum(abs(x_1))>eps) | |
389 syms x_1s | |
390 else | |
391 x_1s=0; | |
392 end | |
393 | |
394 if(sum(abs(x_2))>eps) | |
395 syms x_2s; | |
396 else | |
397 x_2s=0; | |
398 end | |
399 | |
400 | |
401 if(sum(abs(y_1))>eps) | |
402 syms y_1s | |
403 else | |
404 y_1s=0; | |
405 end | |
406 | |
407 if(sum(abs(y_2))>eps) | |
408 syms y_2s; | |
409 else | |
410 y_2s=0; | |
411 end | |
412 | |
413 | |
414 if(sum(abs(z_1))>eps) | |
415 syms z_1s | |
416 else | |
417 z_1s=0; | |
418 end | |
419 | |
420 if(sum(abs(z_2))>eps) | |
421 syms z_2s; | |
422 else | |
423 z_2s=0; | |
424 end | |
425 | |
426 syms xs ys zs | |
427 [V, D]=eig(mat(obj,xs,ys,zs,x_1s,x_2s,y_1s,y_2s,z_1s,z_2s)); | |
428 Vi=inv(V); | |
429 | |
430 syms x_1s x_2s y_1s y_2s z_1s z_2s | |
431 % V= matlabFunction(V); | |
432 % D= matlabFunction(D); | |
433 % Vi= matlabFunction(Vi); | |
434 % | |
435 % xs=x; | |
436 % ys=y; | |
437 % zs=z; | |
438 % x_1s=x_1; | |
439 % x_2s=x_2; | |
440 % y_1s=y_1; | |
441 % y_2s=y_2; | |
442 % z_1s=z_1; | |
443 % z_2s=z_2; | |
444 | |
445 side=max(length(x),length(y)); | |
446 Dret=zeros(obj.n,side*obj.n); | |
447 Vret=zeros(obj.n,side*obj.n); | |
448 Viret=zeros(obj.n,side*obj.n); | |
449 | |
450 for ii=1:obj.n | |
451 for jj=1:obj.n | |
452 Dpart=matlabFunction(D(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); | |
453 Vpart=matlabFunction(V(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); | |
454 Vipart=matlabFunction(V(jj,ii),'Vars',[xs ys zs x_1s x_2s y_1s y_2s z_1s z_2s]); | |
455 Dret(jj,(ii-1)*side+1:side*ii)=sparse(Dpart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)); | |
456 Vret(jj,(ii-1)*side+1:side*ii)=sparse(Vpart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)); | |
457 Viret(jj,(ii-1)*side+1:side*ii)=sparse(Vipart(x,y,z,x_1,x_2,y_1,y_2,z_1,z_2)); | |
458 end | |
459 end | |
460 | |
461 D=sparse(Dret); | |
462 V=sparse(Vret); | |
463 Vi=sparse(Viret); | |
464 V=obj.evaluateCoefficientMatrix(V,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); | |
465 D=obj.evaluateCoefficientMatrix(D,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); | |
466 Vi=obj.evaluateCoefficientMatrix(Vi,x,y,z,x_1,x_2,y_1,y_2,z_1,z_2); | |
467 DD=diag(D); | |
468 | |
469 poseig=(DD>0); | |
470 zeroeig=(DD==0); | |
471 negeig=(DD<0); | |
472 | |
473 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); | |
474 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; | |
475 %Vi=inv(V); | |
476 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; | |
477 end | |
478 end | |
479 end |