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