Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 1197:433c89bf19e0 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 07 Aug 2019 15:23:42 +0200
parents 8d73fcdb07a5
children
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
167 params = obj.params; 167 params = obj.params;
168 X = obj.X; 168 X = obj.X;
169 Y = obj.Y; 169 Y = obj.Y;
170 xi = obj.xi; 170 xi = obj.xi;
171 eta = obj.eta; 171 eta = obj.eta;
172 e_ = obj.getBoundaryOperator('e', boundary);
172 173
173 switch boundary 174 switch boundary
174 case {'w','W','west'} 175 case {'w','W','west'}
175 e_ = obj.e_w;
176 mat = obj.Ahat; 176 mat = obj.Ahat;
177 boundPos = 'l'; 177 boundPos = 'l';
178 Hi = obj.Hxii; 178 Hi = obj.Hxii;
179 [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)); 179 [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));
180 side = max(length(eta)); 180 side = max(length(eta));
181 case {'e','E','east'} 181 case {'e','E','east'}
182 e_ = obj.e_e;
183 mat = obj.Ahat; 182 mat = obj.Ahat;
184 boundPos = 'r'; 183 boundPos = 'r';
185 Hi = obj.Hxii; 184 Hi = obj.Hxii;
186 [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)); 185 [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));
187 side = max(length(eta)); 186 side = max(length(eta));
188 case {'s','S','south'} 187 case {'s','S','south'}
189 e_ = obj.e_s;
190 mat = obj.Bhat; 188 mat = obj.Bhat;
191 boundPos = 'l'; 189 boundPos = 'l';
192 Hi = obj.Hetai; 190 Hi = obj.Hetai;
193 [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)); 191 [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));
194 side = max(length(xi)); 192 side = max(length(xi));
195 case {'n','N','north'} 193 case {'n','N','north'}
196 e_ = obj.e_n;
197 mat = obj.Bhat; 194 mat = obj.Bhat;
198 boundPos = 'r'; 195 boundPos = 'r';
199 Hi = obj.Hetai; 196 Hi = obj.Hetai;
200 [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)); 197 [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));
201 side = max(length(xi)); 198 side = max(length(xi));
372 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]); 369 D = diag([DD(poseig); DD(zeroeig); DD(negeig)]);
373 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; 370 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
374 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 371 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
375 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; 372 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
376 end 373 end
374
375 % Returns the boundary operator op for the boundary specified by the string boundary.
376 % op -- string or a cell array of strings
377 % boundary -- string
378 function varargout = getBoundaryOperator(obj, op, boundary)
379 assertIsMember(boundary, {'w', 'e', 's', 'n'})
380
381 if ~iscell(op)
382 op = {op};
383 end
384
385 for i = 1:numel(op)
386 switch op{i}
387 case 'e'
388 switch boundary
389 case 'w'
390 e = obj.e_w;
391 case 'e'
392 e = obj.e_e;
393 case 's'
394 e = obj.e_s;
395 case 'n'
396 e = obj.e_n;
397 end
398 varargout{i} = e;
399 end
400 end
401 end
402
403 % Returns square boundary quadrature matrix, of dimension
404 % corresponding to the number of boundary points
405 %
406 % boundary -- string
407 function H_b = getBoundaryQuadrature(obj, boundary)
408 assertIsMember(boundary, {'w', 'e', 's', 'n'})
409
410 e = obj.getBoundaryOperator('e', boundary);
411
412 switch boundary
413 case 'w'
414 H_b = inv(e'*obj.Hetai*e);
415 case 'e'
416 H_b = inv(e'*obj.Hetai*e);
417 case 's'
418 H_b = inv(e'*obj.Hxii*e);
419 case 'n'
420 H_b = inv(e'*obj.Hxii*e);
421 end
422 end
423
424
377 end 425 end
378 end 426 end