Mercurial > repos > public > sbplib

comparison +scheme/Heat2dVariable.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 21394c78c72e
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
982:a0b3161e44f3 997:78db023a7fe3
1 classdef Heat2dVariable < scheme.Scheme 1 classdef Heat2dVariable < scheme.Scheme
2 2
3 % Discretizes the Laplacian with variable coefficent, 3 % Discretizes the Laplacian with variable coefficent,
4 % In the Heat equation way (i.e., the discretization matrix is not necessarily 4 % In the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric) 5 % symmetric)
6 % u_t = div * (kappa * grad u ) 6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This 7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction. 8 % is useful if we have periodic BC in one direction.
9 9
10 properties 10 properties
11 m % Number of points in each direction, possibly a vector 11 m % Number of points in each direction, possibly a vector
27 27
28 H, Hi % Inner products 28 H, Hi % Inner products
29 e_l, e_r 29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary 30 d1_l, d1_r % Normal derivatives at the boundary
31 alpha % Vector of borrowing constants 31 alpha % Vector of borrowing constants
32 32
33 H_boundary % Boundary inner products 33 H_boundary % Boundary inner products
34 34
35 end 35 end
36 36
37 methods 37 methods
160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning) 160 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
161 default_arg('type','Neumann'); 161 default_arg('type','Neumann');
162 default_arg('symmetric', false); 162 default_arg('symmetric', false);
163 default_arg('tuning',1.2); 163 default_arg('tuning',1.2);
164 164
165 % j is the coordinate direction of the boundary 165 % nj: outward unit normal component.
166 % nj: outward unit normal component.
167 % nj = -1 for west, south, bottom boundaries 166 % nj = -1 for west, south, bottom boundaries
168 % nj = 1 for east, north, top boundaries 167 % nj = 1 for east, north, top boundaries
169 [j, nj] = obj.get_boundary_number(boundary); 168 nj = obj.getBoundarySign(boundary);
170 switch nj
171 case 1
172 e = obj.e_r;
173 d = obj.d1_r;
174 case -1
175 e = obj.e_l;
176 d = obj.d1_l;
177 end
178 169
179 Hi = obj.Hi; 170 Hi = obj.Hi;
180 H_gamma = obj.H_boundary{j}; 171 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary);
172 H_gamma = obj.getBoundaryQuadrature(boundary);
173 alpha = obj.getBoundaryBorrowing(boundary);
174
181 KAPPA = obj.KAPPA; 175 KAPPA = obj.KAPPA;
182 kappa_gamma = e{j}'*KAPPA*e{j}; 176 kappa_gamma = e'*KAPPA*e;
183 h = obj.h(j);
184 alpha = h*obj.alpha(j);
185 177
186 switch type 178 switch type
187 179
188 % Dirichlet boundary condition 180 % Dirichlet boundary condition
189 case {'D','d','dirichlet','Dirichlet'} 181 case {'D','d','dirichlet','Dirichlet'}
190 182
191 if ~symmetric 183 if ~symmetric
192 closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); 184 closure = -nj*Hi*d*kappa_gamma*H_gamma*(e' );
193 penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; 185 penalty = nj*Hi*d*kappa_gamma*H_gamma;
194 else 186 else
195 closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )... 187 closure = nj*Hi*d*kappa_gamma*H_gamma*(e' )...
196 -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ; 188 -tuning*2/alpha*Hi*e*kappa_gamma*H_gamma*(e' ) ;
197 penalty = -nj*Hi*d{j}*kappa_gamma*H_gamma ... 189 penalty = -nj*Hi*d*kappa_gamma*H_gamma ...
198 +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma; 190 +tuning*2/alpha*Hi*e*kappa_gamma*H_gamma;
199 end 191 end
200 192
201 % Free boundary condition 193 % Free boundary condition
202 case {'N','n','neumann','Neumann'} 194 case {'N','n','neumann','Neumann'}
203 closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); 195 closure = -nj*Hi*e*kappa_gamma*H_gamma*(d' );
204 penalty = Hi*e{j}*kappa_gamma*H_gamma; 196 penalty = Hi*e*kappa_gamma*H_gamma;
205 % penalty is for normal derivative and not for derivative, hence the sign. 197 % penalty is for normal derivative and not for derivative, hence the sign.
206 198
207 % Unknown boundary condition 199 % Unknown boundary condition
208 otherwise 200 otherwise
209 error('No such boundary condition: type = %s',type); 201 error('No such boundary condition: type = %s',type);
214 % u denotes the solution in the own domain 206 % u denotes the solution in the own domain
215 % v denotes the solution in the neighbour domain 207 % v denotes the solution in the neighbour domain
216 error('Interface not implemented'); 208 error('Interface not implemented');
217 end 209 end
218 210
219 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 211 % Returns the boundary operator op for the boundary specified by the string boundary.
220 function [j, nj] = get_boundary_number(obj, boundary) 212 % op -- string or a cell array of strings
213 % boundary -- string
214 function varargout = getBoundaryOperator(obj, op, boundary)
215
216 if ~iscell(op)
217 op = {op};
218 end
219
220 for i = 1:numel(op)
221 switch op{i}
222 case 'e'
223 switch boundary
224 case 'w'
225 e = obj.e_l{1};
226 case 'e'
227 e = obj.e_r{1};
228 case 's'
229 e = obj.e_l{2};
230 case 'n'
231 e = obj.e_r{2};
232 otherwise
233 error('No such boundary: boundary = %s',boundary);
234 end
235 varargout{i} = e;
236
237 case 'd'
238 switch boundary
239 case 'w'
240 d = obj.d1_l{1};
241 case 'e'
242 d = obj.d1_r{1};
243 case 's'
244 d = obj.d1_l{2};
245 case 'n'
246 d = obj.d1_r{2};
247 otherwise
248 error('No such boundary: boundary = %s',boundary);
249 end
250 varargout{i} = d;
251 end
252 end
253 end
254
255 % Returns square boundary quadrature matrix, of dimension
256 % corresponding to the number of boundary points
257 %
258 % boundary -- string
259 function H_b = getBoundaryQuadrature(obj, boundary)
221 260
222 switch boundary 261 switch boundary
223 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 262 case 'w'
224 j = 1; 263 H_b = obj.H_boundary{1};
225 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 264 case 'e'
226 j = 2; 265 H_b = obj.H_boundary{1};
266 case 's'
267 H_b = obj.H_boundary{2};
268 case 'n'
269 H_b = obj.H_boundary{2};
227 otherwise 270 otherwise
228 error('No such boundary: boundary = %s',boundary); 271 error('No such boundary: boundary = %s',boundary);
229 end 272 end
230 273 end
274
275 % Returns the boundary sign. The right boundary is considered the positive boundary
276 % boundary -- string
277 function s = getBoundarySign(obj, boundary)
231 switch boundary 278 switch boundary
232 case {'w','W','west','West','s','S','south','South'} 279 case {'e','n'}
233 nj = -1; 280 s = 1;
234 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 281 case {'w','s'}
235 nj = 1; 282 s = -1;
236 end
237 end
238
239 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
240 function [return_op] = get_boundary_operator(obj, op, boundary)
241
242 switch boundary
243 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
244 j = 1;
245 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
246 j = 2;
247 otherwise 283 otherwise
248 error('No such boundary: boundary = %s',boundary); 284 error('No such boundary: boundary = %s',boundary);
249 end 285 end
250 286 end
251 switch op 287
252 case 'e' 288 % Returns borrowing constant gamma*h
253 switch boundary 289 % boundary -- string
254 case {'w','W','west','West','s','S','south','South'} 290 function gamm = getBoundaryBorrowing(obj, boundary)
255 return_op = obj.e_l{j}; 291 switch boundary
256 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 292 case {'w','e'}
257 return_op = obj.e_r{j}; 293 gamm = obj.h(1)*obj.alpha(1);
258 end 294 case {'s','n'}
259 case 'd' 295 gamm = obj.h(2)*obj.alpha(2);
260 switch boundary
261 case {'w','W','west','West','s','S','south','South'}
262 return_op = obj.d1_l{j};
263 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
264 return_op = obj.d1_r{j};
265 end
266 otherwise 296 otherwise
267 error(['No such operator: operatr = ' op]); 297 error('No such boundary: boundary = %s',boundary);
268 end 298 end
269
270 end 299 end
271 300
272 function N = size(obj) 301 function N = size(obj)
273 N = prod(obj.m); 302 N = prod(obj.m);
274 end 303 end