Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 301:d9860ebc3148 feature/hypsyst

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
HypsystCurve2D Seems to work (Converges with MMS)
author Ylva Rydin <ylva.rydin@telia.com>
date Wed, 05 Oct 2016 17:36:34 +0200
parents 32f3ee81bc37
children 5d5652fe826a
comparison
equal deleted inserted replaced
300:32f3ee81bc37 301:d9860ebc3148
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 n %size of system 4 n %size of system
5 h % Grid spacing 5 h % Grid spacing
6 X,Y % Values of x and y for each grid point 6 X,Y % Values of x and y for each grid point
7 7
8 J, Ji 8 J, Ji %Jacobaian and inverse Jacobian
9 xi,eta 9 xi,eta
10 Xi,Eta 10 Xi,Eta
11 11
12 A,B 12 A,B
13 X_eta, Y_eta 13 X_eta, Y_eta
16 16
17 D % non-stabalized scheme operator 17 D % non-stabalized scheme operator
18 Ahat, Bhat, E 18 Ahat, Bhat, E
19 19
20 H % Discrete norm 20 H % Discrete norm
21 % Norms in the x and y directions 21 Hxii,Hetai % Kroneckerd norms in xi and eta.
22 Hxii,Hetai % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
23 I_xi,I_eta, I_N, onesN 22 I_xi,I_eta, I_N, onesN
24 e_w, e_e, e_s, e_n 23 e_w, e_e, e_s, e_n
25 index_w, index_e,index_s,index_n 24 index_w, index_e,index_s,index_n
26 params %parameters for the coeficient matrice 25 params %parameters for the coeficient matrice
27 end 26 end
37 m = [m m]; 36 m = [m m];
38 end 37 end
39 obj.params = params; 38 obj.params = params;
40 obj.A=A; 39 obj.A=A;
41 obj.B=B; 40 obj.B=B;
42 41
43
44 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta); 42 obj.Ahat=@(params,x,y,x_eta,y_eta)(A(params,x,y).*y_eta-B(params,x,y).*x_eta);
45 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi); 43 obj.Bhat=@(params,x,y,x_xi,y_xi)(B(params,x,y).*x_xi-A(params,x,y).*y_xi);
46 obj.E=@(params,x,y,~,~)E(params,x,y); 44 obj.E=@(params,x,y,~,~)E(params,x,y);
47 45
48 m_xi = m(1); 46 m_xi = m(1);
59 obj.Eta = kr(ones(m_xi,1),obj.eta); 57 obj.Eta = kr(ones(m_xi,1),obj.eta);
60 58
61 obj.n = length(A(obj.params,0,0)); 59 obj.n = length(A(obj.params,0,0));
62 obj.onesN=ones(obj.n); 60 obj.onesN=ones(obj.n);
63 61
64 obj.index_w=1:m_xi; 62 obj.index_w=1:m_eta;
65 obj.index_e=(m_tot-m_xi+1):m_tot; 63 obj.index_e=(m_tot-m_eta+1):m_tot;
66 obj.index_s=1:m_xi:(m_tot-m_xi+1); 64 obj.index_s=1:m_eta:(m_tot-m_eta+1);
67 obj.index_n=(m_xi):m_xi:m_tot; 65 obj.index_n=(m_eta):m_eta:m_tot;
68 66
69 I_n = eye(obj.n); 67 I_n = eye(obj.n);
70 I_xi = speye(m_xi); 68 I_xi = speye(m_xi);
71 obj.I_xi = I_xi; 69 obj.I_xi = I_xi;
72 I_eta = speye(m_eta); 70 I_eta = speye(m_eta);
136 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_) 134 function [ret] = evaluateCoefficientMatrix(obj, mat, X, Y,x_,y_)
137 params=obj.params; 135 params=obj.params;
138 136
139 if isa(mat,'function_handle') 137 if isa(mat,'function_handle')
140 [rows,cols]=size(mat(params,0,0,0,0)); 138 [rows,cols]=size(mat(params,0,0,0,0));
141 x_=kr(x_,obj.onesN); 139 x_=kr(obj.onesN,x_);
142 y_=kr(y_,obj.onesN); 140 y_=kr(obj.onesN,y_);
143 matVec=mat(params,X',Y',x_',y_'); 141 matVec=mat(params,X',Y',x_',y_');
144 matVec=sparse(matVec); 142 matVec=sparse(matVec);
145 side=max(length(X),length(Y)); 143 side=max(length(X),length(Y));
146 else 144 else
147 matVec=mat; 145 matVec=mat;
161 159
162 function [closure, penalty]=boundary_condition_char(obj,boundary) 160 function [closure, penalty]=boundary_condition_char(obj,boundary)
163 params=obj.params; 161 params=obj.params;
164 X=obj.X; Y=obj.Y; 162 X=obj.X; Y=obj.Y;
165 xi=obj.xi; eta=obj.eta; 163 xi=obj.xi; eta=obj.eta;
166 side=max(length(xi),length(eta)); 164
167 165
168 switch boundary 166 switch boundary
169 case {'w','W','west'} 167 case {'w','W','west'}
170 e_=obj.e_w; 168 e_=obj.e_w;
171 mat=obj.Ahat; 169 mat=obj.Ahat;
172 boundPos='l'; 170 boundPos='l';
173 Hi=obj.Hxii; 171 Hi=obj.Hxii;
174 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); 172 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
173 side=max(length(eta));
175 case {'e','E','east'} 174 case {'e','E','east'}
176 e_=obj.e_e; 175 e_=obj.e_e;
177 mat=obj.Ahat; 176 mat=obj.Ahat;
178 boundPos='r'; 177 boundPos='r';
179 Hi=obj.Hxii; 178 Hi=obj.Hxii;
180 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); 179 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
180 side=max(length(eta));
181 case {'s','S','south'} 181 case {'s','S','south'}
182 e_=obj.e_s; 182 e_=obj.e_s;
183 mat=obj.Bhat; 183 mat=obj.Bhat;
184 boundPos='l'; 184 boundPos='l';
185 Hi=obj.Hetai; 185 Hi=obj.Hetai;
186 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); 186 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
187 side=max(length(xi));
187 case {'n','N','north'} 188 case {'n','N','north'}
188 e_=obj.e_n; 189 e_=obj.e_n;
189 mat=obj.Bhat; 190 mat=obj.Bhat;
190 boundPos='r'; 191 boundPos='r';
191 Hi=obj.Hetai; 192 Hi=obj.Hetai;
192 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 193 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
194 side=max(length(xi));
193 end 195 end
194 196
195 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 197 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
196 198
197 switch boundPos 199 switch boundPos
212 214
213 215
214 function [closure,penalty]=boundary_condition_general(obj,boundary,L) 216 function [closure,penalty]=boundary_condition_general(obj,boundary,L)
215 params=obj.params; 217 params=obj.params;
216 X=obj.X; Y=obj.Y; 218 X=obj.X; Y=obj.Y;
217 side=max(length(xi),length(eta)); 219 xi=obj.xi; eta=obj.eta;
218 220
219 switch boundary 221 switch boundary
220 case {'w','W','west'} 222 case {'w','W','west'}
221 e_=obj.e_w; 223 e_=obj.e_w;
222 mat=obj.Ahat; 224 mat=obj.Ahat;
223 boundPos='l'; 225 boundPos='l';
224 Hi=obj.Hxii; 226 Hi=obj.Hxii;
225 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); 227 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
226 L=obj.evaluateCoefficientMatrix(L,X(obj.index_w),Y(obj.index_w)); 228
229 Ji_vec=diag(obj.Ji);
230 Ji=diag(Ji_vec(obj.index_w));
231 xi_x=Ji*obj.Y_eta(obj.index_w);
232 xi_y=-Ji*obj.X_eta(obj.index_w);
233 L=obj.evaluateCoefficientMatrix(L,xi_x,xi_y,[],[]);
234 side=max(length(eta));
227 case {'e','E','east'} 235 case {'e','E','east'}
228 e_=obj.e_e; 236 e_=obj.e_e;
229 mat=obj.Ahat; 237 mat=obj.Ahat;
230 boundPos='r'; 238 boundPos='r';
231 Hi=obj.Hxii; 239 Hi=obj.Hxii;
232 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); 240 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
233 L=obj.evaluateCoefficientMatrix(L,X(obj.index_e),Y(obj.index_e),[],[]); 241
242 Ji_vec=diag(obj.Ji);
243 Ji=diag(Ji_vec(obj.index_e));
244 xi_x=Ji*obj.Y_eta(obj.index_e);
245 xi_y=-Ji*obj.X_eta(obj.index_e);
246 L=obj.evaluateCoefficientMatrix(L,-xi_x,-xi_y,[],[]);
247 side=max(length(eta));
234 case {'s','S','south'} 248 case {'s','S','south'}
235 e_=obj.e_s; 249 e_=obj.e_s;
236 mat=obj.Bhat; 250 mat=obj.Bhat;
237 boundPos='l'; 251 boundPos='l';
238 Hi=obj.Hetai; 252 Hi=obj.Hetai;
239 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); 253 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
240 L=obj.evaluateCoefficientMatrix(L,X(obj.index_s),Y(obj.index_s),[],[]); 254
255 Ji_vec=diag(obj.Ji);
256 Ji=diag(Ji_vec(obj.index_s));
257 eta_x=Ji*obj.Y_xi(obj.index_s);
258 eta_y=-Ji*obj.X_xi(obj.index_s);
259 L=obj.evaluateCoefficientMatrix(L,eta_x,eta_y,[],[]);
260 side=max(length(xi));
241 case {'n','N','north'} 261 case {'n','N','north'}
242 e_=obj.e_n; 262 e_=obj.e_n;
243 mat=obj.Bhat; 263 mat=obj.Bhat;
244 boundPos='r'; 264 boundPos='r';
245 Hi=obj.Hetai; 265 Hi=obj.Hetai;
246 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n)); 266 [V,Vi,D,signVec]=obj.matrixDiag(mat,X(obj.index_n),Y(obj.index_n),obj.X_xi(obj.index_n),obj.Y_xi(obj.index_n));
247 L=obj.evaluateCoefficientMatrix(L,X(obj.index_n),Y(obj.index_n)); 267
268 Ji_vec=diag(obj.Ji);
269 Ji=diag(Ji_vec(obj.index_n));
270 eta_x=Ji*obj.Y_xi(obj.index_n);
271 eta_y=-Ji*obj.X_xi(obj.index_n);
272 L=obj.evaluateCoefficientMatrix(L,-eta_x,-eta_y,[],[]);
273
274 side=max(length(xi));
275
248 end 276 end
249 277
250 pos=signVec(1); zeroval=signVec(2); neg=signVec(3); 278 pos=signVec(1); zeroval=signVec(2); neg=signVec(3);
251 279
252 switch boundPos 280 switch boundPos
253 case {'l'} 281 case {'l'}
254 tau=sparse(obj.n*side,pos); 282 tau=sparse(obj.n*side,pos);
255 Vi_plus=Vi(1:pos,:); 283 Vi_plus=Vi(1:pos,:);
256 Vi_minus=Vi(pos+1:obj.n*side,:); 284 Vi_minus=Vi(pos+1:obj.n*side,:);
257 V_plus=V(:,1:pos); 285 V_plus=V(:,1:pos);
258 V_minus=V(:,(pos+zeroval)+1:obj.n*side); 286 V_minus=V(:,(pos)+1:obj.n*side);
259 287
260 tau(1:pos*side,:)=-abs(D(1:pos,1:pos)); 288 tau(1:pos,:)=-abs(D(1:pos,1:pos));
261 R=-inv(L*V_plus)*(L*V_minus); 289 R=-inv(L*V_plus)*(L*V_minus);
262 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_'; 290 closure=Hi*e_*V*tau*(Vi_plus-R*Vi_minus)*e_';
263 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L; 291 penalty=-Hi*e_*V*tau*inv(L*V_plus)*L;
264 case {'r'} 292 case {'r'}
265 tau=sparse(obj.n*side,neg); 293 tau=sparse(obj.n*side,neg);