Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2dCurve.m @ 997:78db023a7fe3 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 12 Jan 2019 11:57:50 -0800
parents 706d1c2b4199
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
982:a0b3161e44f3 997:78db023a7fe3
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
380 if ~iscell(op)
381 op = {op};
382 end
383
384 for i = 1:numel(op)
385 switch op{i}
386 case 'e'
387 switch boundary
388 case 'w'
389 e = obj.e_w;
390 case 'e'
391 e = obj.e_e;
392 case 's'
393 e = obj.e_s;
394 case 'n'
395 e = obj.e_n;
396 otherwise
397 error('No such boundary: boundary = %s',boundary);
398 end
399 varargout{i} = e;
400 end
401 end
402 end
403
404 % Returns square boundary quadrature matrix, of dimension
405 % corresponding to the number of boundary points
406 %
407 % boundary -- string
408 function H_b = getBoundaryQuadrature(obj, boundary)
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 otherwise
422 error('No such boundary: boundary = %s',boundary);
423 end
424 end
425
426
377 end 427 end
378 end 428 end