Mercurial > repos > public > sbplib

comparison +scheme/Hypsyst2d.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
184 % General boundary condition in the form Lu=g(x) 184 % General boundary condition in the form Lu=g(x)
185 function [closure,penalty] = boundary_condition_general(obj,boundary,L) 185 function [closure,penalty] = boundary_condition_general(obj,boundary,L)
186 params = obj.params; 186 params = obj.params;
187 x = obj.x; 187 x = obj.x;
188 y = obj.y; 188 y = obj.y;
189 e_ = obj.getBoundaryOperator('e', boundary);
189 190
190 switch boundary 191 switch boundary
191 case {'w','W','west'} 192 case {'w','W','west'}
192 e_ = obj.e_w;
193 mat = obj.A; 193 mat = obj.A;
194 boundPos = 'l'; 194 boundPos = 'l';
195 Hi = obj.Hxi; 195 Hi = obj.Hxi;
196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y); 196 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(1),y);
197 L = obj.evaluateCoefficientMatrix(L,x(1),y); 197 L = obj.evaluateCoefficientMatrix(L,x(1),y);
198 side = max(length(y)); 198 side = max(length(y));
199 case {'e','E','east'} 199 case {'e','E','east'}
200 e_ = obj.e_e;
201 mat = obj.A; 200 mat = obj.A;
202 boundPos = 'r'; 201 boundPos = 'r';
203 Hi = obj.Hxi; 202 Hi = obj.Hxi;
204 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y); 203 [V,Vi,D,signVec] = obj.matrixDiag(mat,x(end),y);
205 L = obj.evaluateCoefficientMatrix(L,x(end),y); 204 L = obj.evaluateCoefficientMatrix(L,x(end),y);
206 side = max(length(y)); 205 side = max(length(y));
207 case {'s','S','south'} 206 case {'s','S','south'}
208 e_ = obj.e_s;
209 mat = obj.B; 207 mat = obj.B;
210 boundPos = 'l'; 208 boundPos = 'l';
211 Hi = obj.Hyi; 209 Hi = obj.Hyi;
212 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1)); 210 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(1));
213 L = obj.evaluateCoefficientMatrix(L,x,y(1)); 211 L = obj.evaluateCoefficientMatrix(L,x,y(1));
214 side = max(length(x)); 212 side = max(length(x));
215 case {'n','N','north'} 213 case {'n','N','north'}
216 e_ = obj.e_n;
217 mat = obj.B; 214 mat = obj.B;
218 boundPos = 'r'; 215 boundPos = 'r';
219 Hi = obj.Hyi; 216 Hi = obj.Hyi;
220 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end)); 217 [V,Vi,D,signVec] = obj.matrixDiag(mat,x,y(end));
221 L = obj.evaluateCoefficientMatrix(L,x,y(end)); 218 L = obj.evaluateCoefficientMatrix(L,x,y(end));
295 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)]; 292 V = [V(:,poseig) V(:,zeroeig) V(:,negeig)];
296 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; 293 Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
297 signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; 294 signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
298 end 295 end
299 296
297 % Returns the boundary operator op for the boundary specified by the string boundary.
298 % op -- string or a cell array of strings
299 % boundary -- string
300 function varargout = getBoundaryOperator(obj, op, boundary)
301
302 if ~iscell(op)
303 op = {op};
304 end
305
306 for i = 1:numel(op)
307 switch op{i}
308 case 'e'
309 switch boundary
310 case 'w'
311 e = obj.e_w;
312 case 'e'
313 e = obj.e_e;
314 case 's'
315 e = obj.e_s;
316 case 'n'
317 e = obj.e_n;
318 otherwise
319 error('No such boundary: boundary = %s',boundary);
320 end
321 varargout{i} = e;
322 end
323 end
324 end
325
326 % Returns square boundary quadrature matrix, of dimension
327 % corresponding to the number of boundary points
328 %
329 % boundary -- string
330 function H_b = getBoundaryQuadrature(obj, boundary)
331
332 e = obj.getBoundaryOperator('e', boundary);
333
334 switch boundary
335 case 'w'
336 H_b = inv(e'*obj.Hyi*e);
337 case 'e'
338 H_b = inv(e'*obj.Hyi*e);
339 case 's'
340 H_b = inv(e'*obj.Hxi*e);
341 case 'n'
342 H_b = inv(e'*obj.Hxi*e);
343 otherwise
344 error('No such boundary: boundary = %s',boundary);
345 end
346 end
347
300 end 348 end
301 end 349 end