Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst3d.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 9b3d7fc61a36
children 69b078cf8072
comparison
equal deleted inserted replaced
353:fd4f1c80755d 354:dbac99d2c318
12 A, B, C, E 12 A, B, C, E
13 Aevaluated,Bevaluated,Cevaluated, Eevaluated 13 Aevaluated,Bevaluated,Cevaluated, Eevaluated
14 14
15 H % Discrete norm 15 H % Discrete norm
16 % Norms in the x, y and z directions 16 % Norms in the x, y and z directions
17 Hx, Hy, Hz
17 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 18 Hxi,Hyi, Hzi % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
18 I_x,I_y, I_z, I_N 19 I_x,I_y, I_z, I_N
19 e_w, e_e, e_s, e_n, e_b, e_t 20 e_w, e_e, e_s, e_n, e_b, e_t
20 params %parameters for the coeficient matrice 21 params %parameters for the coeficient matrice
21 end 22 end
22 23
23 24
24 methods 25 methods
25 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params) 26 function obj = Hypsyst3d(m, lim, order, A, B,C, E, params,operator)
26 default_arg('E', []) 27 default_arg('E', [])
28 default_arg('operatpr',[])
27 xlim = lim{1}; 29 xlim = lim{1};
28 ylim = lim{2}; 30 ylim = lim{2};
29 zlim = lim{3}; 31 zlim = lim{3};
30 32
31 if length(m) == 1 33 if length(m) == 1
32 m = [m m m]; 34 m = [m m m];
33 end 35 end
34 36
35 obj.A=A; 37 obj.A = A;
36 obj.B=B; 38 obj.B = B;
37 obj.C=C; 39 obj.C = C;
38 obj.E=E; 40 obj.E = E;
39 m_x = m(1); 41 m_x = m(1);
40 m_y = m(2); 42 m_y = m(2);
41 m_z=m(3); 43 m_z=m(3);
42 obj.params = params; 44 obj.params = params;
43 45
44 ops_x = sbp.D2Standard(m_x,xlim,order); 46 switch operator
45 ops_y = sbp.D2Standard(m_y,ylim,order); 47 case 'upwind'
46 ops_z = sbp.D2Standard(m_z,zlim,order); 48 ops_x = sbp.D1Upwind(m_x,xlim,order);
49 ops_y = sbp.D1Upwind(m_y,ylim,order);
50 ops_z = sbp.D1Upwind(m_z,zlim,order);
51 otherwise
52 ops_x = sbp.D2Standard(m_x,xlim,order);
53 ops_y = sbp.D2Standard(m_y,ylim,order);
54 ops_z = sbp.D2Standard(m_z,zlim,order);
55 end
47 56
48 obj.x = ops_x.x; 57 obj.x = ops_x.x;
49 obj.y = ops_y.x; 58 obj.y = ops_y.x;
50 obj.z = ops_z.x; 59 obj.z = ops_z.x;
51 60
52 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));%% Que pasa? 61 obj.X = kr(obj.x,ones(m_y,1),ones(m_z,1));%% Que pasa?
53 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1)); 62 obj.Y = kr(ones(m_x,1),obj.y,ones(m_z,1));
54 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z); 63 obj.Z = kr(ones(m_x,1),ones(m_y,1),obj.z);
55 64
56 obj.Yx=kr(obj.y,ones(m_z,1)); 65 obj.Yx = kr(obj.y,ones(m_z,1));
57 obj.Zx=kr(ones(m_y,1),obj.z); 66 obj.Zx = kr(ones(m_y,1),obj.z);
58 67
59 obj.Xy=kr(obj.x,ones(m_z,1)); 68 obj.Xy = kr(obj.x,ones(m_z,1));
60 obj.Zy=kr(ones(m_x,1),obj.z); 69 obj.Zy = kr(ones(m_x,1),obj.z);
61 70
62 obj.Xz=kr(obj.x,ones(m_y,1)); 71 obj.Xz = kr(obj.x,ones(m_y,1));
63 obj.Yz=kr(ones(m_z,1),obj.y); 72 obj.Yz = kr(ones(m_z,1),obj.y);
64 73
65 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z); 74 obj.Aevaluated = obj.evaluateCoefficientMatrix(A, obj.X, obj.Y,obj.Z);
66 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z); 75 obj.Bevaluated = obj.evaluateCoefficientMatrix(B, obj.X, obj.Y,obj.Z);
67 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z); 76 obj.Cevaluated = obj.evaluateCoefficientMatrix(C, obj.X, obj.Y,obj.Z);
68 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z); 77 obj.Eevaluated = obj.evaluateCoefficientMatrix(E, obj.X, obj.Y,obj.Z);
74 obj.I_x = I_x; 83 obj.I_x = I_x;
75 I_y = speye(m_y); 84 I_y = speye(m_y);
76 obj.I_y = I_y; 85 obj.I_y = I_y;
77 I_z = speye(m_z); 86 I_z = speye(m_z);
78 obj.I_z = I_z; 87 obj.I_z = I_z;
79 88 I_N=kr(I_n,I_x,I_y,I_z);
80 89
81 D1_x = kr(I_n, ops_x.D1, I_y,I_z);
82 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z); 90 obj.Hxi = kr(I_n, ops_x.HI, I_y,I_z);
83 D1_y = kr(I_n, I_x, ops_y.D1,I_z); 91 obj.Hx = ops_x.H;
84 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z); 92 obj.Hyi = kr(I_n, I_x, ops_y.HI,I_z);
85 D1_z = kr(I_n, I_x, I_y,ops_z.D1); 93 obj.Hy = ops_y.H
86 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI); 94 obj.Hzi = kr(I_n, I_x,I_y, ops_z.HI);
95 obj.Hz = ops_z.H;
87 96
88 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z); 97 obj.e_w = kr(I_n, ops_x.e_l, I_y,I_z);
89 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z); 98 obj.e_e = kr(I_n, ops_x.e_r, I_y,I_z);
90 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z); 99 obj.e_s = kr(I_n, I_x, ops_y.e_l,I_z);
91 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z); 100 obj.e_n = kr(I_n, I_x, ops_y.e_r,I_z);
92 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l); 101 obj.e_b = kr(I_n, I_x, I_y, ops_z.e_l);
93 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r); 102 obj.e_t = kr(I_n, I_x, I_y, ops_z.e_r);
94 103
95 obj.m=m; 104 obj.m = m;
96 obj.h=[ops_x.h ops_y.h ops_x.h]; 105 obj.h = [ops_x.h ops_y.h ops_x.h];
97 obj.order=order; 106 obj.order = order;
98 107
99 obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated; 108 switch operator
109 case 'upwind'
110 alphaA = max(eig(A(params,obj.x(end),obj.y(end),obj.z(end))));
111 alphaB = max(eig(B(params,obj.x(end),obj.y(end),obj.z(end))));
112 alphaC = max(eig(C(params,obj.x(end),obj.y(end),obj.z(end))));
113
114 Ap = (obj.Aevaluated+alphaA*I_N)/2;
115 Am = (obj.Aevaluated-alphaA*I_N)/2;
116 Bp = (obj.Bevaluated+alphaB*I_N)/2;
117 Bm = (obj.Bevaluated-alphaB*I_N)/2;
118 Cp = (obj.Cevaluated+alphaC*I_N)/2;
119 Cm = (obj.Cevaluated-alphaC*I_N)/2;
120
121 Dpx = kr(I_n, ops_x.Dp, I_y,I_z);
122 Dmx = kr(I_n, ops_x.Dm, I_y,I_z);
123 Dpy = kr(I_n, I_x, ops_y.Dp,I_z);
124 Dmy = kr(I_n, I_x, ops_y.Dm,I_z);
125 Dpz = kr(I_n, I_x, I_y,ops_z.Dp);
126 Dmz = kr(I_n, I_x, I_y,ops_z.Dm);
127
128 obj.D=-Am*Dpx-Ap*Dmx-Bm*Dpy-Bp*Dmy-Cm*Dpz-Cp*Dmz-obj.Eevaluated;
129 otherwise
130 D1_x = kr(I_n, ops_x.D1, I_y,I_z);
131 D1_y = kr(I_n, I_x, ops_y.D1,I_z);
132 D1_z = kr(I_n, I_x, I_y,ops_z.D1);
133 obj.D=-obj.Aevaluated*D1_x-obj.Bevaluated*D1_y-obj.Cevaluated*D1_z-obj.Eevaluated;
134 end
100 end 135 end
101 136
102 % Closure functions return the opertors applied to the own doamin to close the boundary 137 % Closure functions return the opertors applied to the own doamin to close the boundary
103 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. 138 % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin.
104 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 139 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
108 default_arg('type','char'); 143 default_arg('type','char');
109 BM=boundary_matrices(obj,boundary); 144 BM=boundary_matrices(obj,boundary);
110 145
111 switch type 146 switch type
112 case{'c','char'} 147 case{'c','char'}
113 [closure,penalty]=boundary_condition_char(obj,BM); 148 [closure,penalty] = boundary_condition_char(obj,BM);
114 case{'general'} 149 case{'general'}
115 [closure,penalty]=boundary_condition_general(obj,BM,boundary,L); 150 [closure,penalty] = boundary_condition_general(obj,BM,boundary,L);
116 otherwise 151 otherwise
117 error('No such boundary condition') 152 error('No such boundary condition')
118 end 153 end
119 end 154 end
120 155
125 function N = size(obj) 160 function N = size(obj)
126 N = obj.m; 161 N = obj.m;
127 end 162 end
128 163
129 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z) 164 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y, Z)
130 params=obj.params; 165 params = obj.params;
131 side=max(length(X),length(Y)); 166 side = max(length(X),length(Y));
132 if isa(mat,'function_handle') 167 if isa(mat,'function_handle')
133 [rows,cols]=size(mat(params,0,0,0)); 168 [rows,cols] = size(mat(params,0,0,0));
134 matVec=mat(params,X',Y',Z'); 169 matVec = mat(params,X',Y',Z');
135 matVec=sparse(matVec); 170 matVec=sparse(matVec);
136 else 171 else
137 matVec=mat; 172 matVec = mat;
138 [rows,cols]=size(matVec); 173 [rows,cols]=size(matVec);
139 side=max(length(X),length(Y)); 174 side=max(length(X),length(Y));
140 cols=cols/side; 175 cols=cols/side;
141 end 176 end
142 ret=cell(rows,cols); 177 ret=cell(rows,cols);
213 D=BM.D; 248 D=BM.D;
214 e_=BM.e_; 249 e_=BM.e_;
215 250
216 switch BM.boundpos 251 switch BM.boundpos
217 case {'l'} 252 case {'l'}
218 tau=sparse(obj.n*side,pos); 253 tau = sparse(obj.n*side,pos);
219 Vi_plus=Vi(1:pos,:); 254 Vi_plus = Vi(1:pos,:);
220 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 255 tau(1:pos,:) = -abs(D(1:pos,1:pos));
221 closure=Hi*e_*V*tau*Vi_plus*e_'; 256 closure = Hi*e_*V*tau*Vi_plus*e_';
222 penalty=-Hi*e_*V*tau*Vi_plus; 257 penalty = -Hi*e_*V*tau*Vi_plus;
223 case {'r'} 258 case {'r'}
224 tau=sparse(obj.n*side,neg); 259 tau = sparse(obj.n*side,neg);
225 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 260 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
226 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 261 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
227 closure=Hi*e_*V*tau*Vi_minus*e_'; 262 closure = Hi*e_*V*tau*Vi_minus*e_';
228 penalty=-Hi*e_*V*tau*Vi_minus; 263 penalty = -Hi*e_*V*tau*Vi_minus;
229 end 264 end
230 end 265 end
231 266
232 267
233 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L) 268 function [closure,penalty]=boundary_condition_general(obj,BM,boundary,L)
239 Vi = BM.Vi; 274 Vi = BM.Vi;
240 Hi=BM.Hi; 275 Hi=BM.Hi;
241 D=BM.D; 276 D=BM.D;
242 e_=BM.e_; 277 e_=BM.e_;
243 switch boundary 278 switch boundary
244 case {'w','W','west'} 279 case {'w','W','west'}
245 L=obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx); 280 L = obj.evaluateCoefficientMatrix(L,obj.x(1),obj.Yx,obj.Zx);
246 case {'e','E','east'} 281 case {'e','E','east'}
247 L=obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx); 282 L = obj.evaluateCoefficientMatrix(L,obj.x(end),obj.Yx,obj.Zx);
248 case {'s','S','south'} 283 case {'s','S','south'}
249 L=obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy); 284 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(1),obj.Zy);
250 case {'n','N','north'} 285 case {'n','N','north'}
251 L=obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy); 286 L = obj.evaluateCoefficientMatrix(L,obj.Xy,obj.y(end),obj.Zy);
252 case {'b','B','bottom'} 287 case {'b','B','bottom'}
253 L=obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1)); 288 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(1));
254 case {'t','T','top'} 289 case {'t','T','top'}
255 L=obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end)); 290 L = obj.evaluateCoefficientMatrix(L,obj.Xz,obj.Yz,obj.z(end));
256 end 291 end
257 292
258 switch BM.boundpos 293 switch BM.boundpos
259 case {'l'} 294 case {'l'}
260 tau=sparse(obj.n*side,pos); 295 tau = sparse(obj.n*side,pos);
261 Vi_plus=Vi(1:pos,:); 296 Vi_plus = Vi(1:pos,:);
262 Vi_minus=Vi(pos+zeroval+1:obj.n*side,:); 297 Vi_minus = Vi(pos+zeroval+1:obj.n*side,:);
263 V_plus=V(:,1:pos); 298 V_plus = V(:,1:pos);
264 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 299 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
265 300
266 tau(1:pos,:)=-abs(D(1:pos,1:pos)); 301 tau(1:pos,:) = -abs(D(1:pos,1:pos));
267 R=-inv(L*V_plus)*(L*V_minus); 302 R = -inv(L*V_plus)*(L*V_minus);
268 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 303 closure = Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
269 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 304 penalty = -Hi*e_*V*tau*inv(L*V_plus)*L;
270 case {'r'} 305 case {'r'}
271 tau=sparse(obj.n*side,neg); 306 tau = sparse(obj.n*side,neg);
272 tau((pos+zeroval)+1:obj.n*side,:)=-abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side)); 307 tau((pos+zeroval)+1:obj.n*side,:) = -abs(D((pos+zeroval)+1:obj.n*side,(pos+zeroval)+1:obj.n*side));
273 Vi_plus=Vi(1:pos,:); 308 Vi_plus = Vi(1:pos,:);
274 Vi_minus=Vi((pos+zeroval)+1:obj.n*side,:); 309 Vi_minus = Vi((pos+zeroval)+1:obj.n*side,:);
275 310
276 V_plus=V(:,1:pos); 311 V_plus = V(:,1:pos);
277 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 312 V_minus = V(:,(pos+zeroval)+1:obj.n*side);
278 R=-inv(L*V_minus)*(L*V_plus); 313 R = -inv(L*V_minus)*(L*V_plus);
279 closure=Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_'; 314 closure = Hi*e_*V*tau*(Vi_minus-R*Vi_plus)*e_';
280 penalty=-Hi*e_*V*tau*inv(L*V_minus)*L; 315 penalty = -Hi*e_*V*tau*inv(L*V_minus)*L;
281 end 316 end
282 end 317 end
283 318
284 319
285 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z) 320 function [V,Vi, D,signVec]=matrixDiag(obj,mat,x,y,z)
286 params=obj.params; 321 params = obj.params;
287 syms xs ys zs 322 syms xs ys zs
288 [V, D]=eig(mat(params,xs,ys,zs)); 323 [V, D] = eig(mat(params,xs,ys,zs));
289 xs=x; 324 Vi=inv(V);
290 ys=y; 325 xs = x;
291 zs=z; 326 ys = y;
292 327 zs = z;
293 328
294 side=max(length(x),length(y)); 329
295 Dret=zeros(obj.n,side*obj.n); 330 side = max(length(x),length(y));
296 Vret=zeros(obj.n,side*obj.n); 331 Dret = zeros(obj.n,side*obj.n);
332 Vret = zeros(obj.n,side*obj.n);
333 Viret= zeros(obj.n,side*obj.n);
297 for ii=1:obj.n 334 for ii=1:obj.n
298 for jj=1:obj.n 335 for jj=1:obj.n
299 Dret(jj,(ii-1)*side+1:side*ii)=eval(D(jj,ii)); 336 Dret(jj,(ii-1)*side+1:side*ii) = eval(D(jj,ii));
300 Vret(jj,(ii-1)*side+1:side*ii)=eval(V(jj,ii)); 337 Vret(jj,(ii-1)*side+1:side*ii) = eval(V(jj,ii));
338 Viret(jj,(ii-1)*side+1:side*ii) = eval(Vi(jj,ii));
301 end 339 end
302 end 340 end
303 341
304 D=sparse(Dret); 342 D = sparse(Dret);
305 V=sparse(Vret); 343 V = sparse(Vret);
306 V=obj.evaluateCoefficientMatrix(V,x,y,z); 344 Vi = sparse(Viret);
307 D=obj.evaluateCoefficientMatrix(D,x,y,z); 345 V = obj.evaluateCoefficientMatrix(V,x,y,z);
308 DD=diag(D); 346 Vi= obj.evaluateCoefficientMatrix(Vi,x,y,z);
309 347 D = obj.evaluateCoefficientMatrix(D,x,y,z);
310 poseig=(DD>0); 348 DD = diag(D);
311 zeroeig=(DD==0); 349
312 negeig=(DD<0); 350 poseig = (DD>0);
313 351 zeroeig = (DD==0);
314 D=diag([DD(poseig); DD(zeroeig); DD(negeig)]); 352 negeig = (DD<0);
315 V=[V(:,poseig) V(:,zeroeig) V(:,negeig)]; 353
316 Vi=inv(V); 354 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
317 signVec=[sum(poseig),sum(zeroeig),sum(negeig)]; 355 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
356 Vi= [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
357 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
318 end 358 end
319 end 359 end
320 end 360 end