Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 861:607c631f175e feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Small changes to Elastic2dVariable to facilitate adjoing gradient computation.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 24 Oct 2018 16:17:32 -0700
parents 5751262b323b
children 14fee299ada2
comparison
equal deleted inserted replaced
860:b758d1cf4c8e 861:607c631f175e
41 H_boundary % Boundary inner products 41 H_boundary % Boundary inner products
42 42
43 % Kroneckered norms and coefficients 43 % Kroneckered norms and coefficients
44 RHOi_kron 44 RHOi_kron
45 Hi_kron 45 Hi_kron
46
47 % Structures used for adjoint optimization
48 B
46 end 49 end
47 50
48 methods 51 methods
49 52
50 function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) 53 % The coefficients can either be function handles or grid functions
54 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
51 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); 55 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
52 default_arg('lambda_fun', @(x,y) 0*x+1); 56 default_arg('lambda', @(x,y) 0*x+1);
53 default_arg('mu_fun', @(x,y) 0*x+1); 57 default_arg('mu', @(x,y) 0*x+1);
54 default_arg('rho_fun', @(x,y) 0*x+1); 58 default_arg('rho', @(x,y) 0*x+1);
55 dim = 2; 59 dim = 2;
56 60
57 assert(isa(g, 'grid.Cartesian')) 61 assert(isa(g, 'grid.Cartesian'))
58 62
59 lambda = grid.evalOn(g, lambda_fun); 63 if isa(lambda, 'function_handle')
60 mu = grid.evalOn(g, mu_fun); 64 lambda = grid.evalOn(g, lambda);
61 rho = grid.evalOn(g, rho_fun); 65 end
66 if isa(mu, 'function_handle')
67 mu = grid.evalOn(g, mu);
68 end
69 if isa(rho, 'function_handle')
70 rho = grid.evalOn(g, rho);
71 end
72
62 m = g.size(); 73 m = g.size();
63 m_tot = g.N(); 74 m_tot = g.N();
64 75
65 h = g.scaling(); 76 h = g.scaling();
66 lim = g.lim; 77 lim = g.lim;
200 D2_mu{j}*E{i}' ... 211 D2_mu{j}*E{i}' ...
201 ); 212 );
202 end 213 end
203 end 214 end
204 obj.D = D; 215 obj.D = D;
205 %=========================================% 216 %=========================================%'
206 217
207 % Numerical traction operators for BC. 218 % Numerical traction operators for BC.
208 % Because d1 =/= e0^T*D1, the numerical tractions are different 219 % Because d1 =/= e0^T*D1, the numerical tractions are different
209 % at every boundary. 220 % at every boundary.
210 T_l = cell(dim,1); 221 T_l = cell(dim,1);
261 obj.m = m; 272 obj.m = m;
262 obj.h = h; 273 obj.h = h;
263 obj.order = order; 274 obj.order = order;
264 obj.grid = g; 275 obj.grid = g;
265 obj.dim = dim; 276 obj.dim = dim;
277
278 % Used for adjoint optimization
279 obj.B = cell(1,dim);
280 for i = 1:dim
281 obj.B{i} = zeros(m(i),m(i),m(i));
282 for k = 1:m(i)
283 c = sparse(m(i),1);
284 c(k) = 1;
285 [~, obj.B{i}(:,:,k)] = ops{i}.D2(c);
286 end
287 end
266 288
267 end 289 end
268 290
269 291
270 % Closure functions return the operators applied to the own domain to close the boundary 292 % Closure functions return the operators applied to the own domain to close the boundary
494 case {'w','W','west','West','s','S','south','South'} 516 case {'w','W','west','West','s','S','south','South'}
495 varargout{i} = obj.tau_l{j}; 517 varargout{i} = obj.tau_l{j};
496 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 518 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
497 varargout{i} = obj.tau_r{j}; 519 varargout{i} = obj.tau_r{j};
498 end 520 end
521 case 'alpha'
522 % alpha = alpha(i,j) is the penalty strength for displacement BC.
523 tuning = 1.2;
524 LAMBDA = obj.LAMBDA;
525 MU = obj.MU;
526
527 phi = obj.phi{j};
528 h = obj.h(j);
529 h11 = obj.H11{j}*h;
530 gamma = obj.gamma{j};
531 dim = obj.dim;
532
533 a_lambda = dim/h11 + 1/(h11*phi);
534 a_mu_i = 2/(gamma*h);
535 a_mu_ij = 2/h11 + 1/(h11*phi);
536
537 d = @kroneckerDelta; % Kronecker delta
538 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
539 alpha = @(i,k) d(i,k)*tuning*( d(i,j)* a_lambda*LAMBDA ...
540 + d(i,j)* a_mu_i*MU ...
541 + db(i,j)*a_mu_ij*MU );
542 varargout{i} = alpha;
499 otherwise 543 otherwise
500 error(['No such operator: operator = ' op{i}]); 544 error(['No such operator: operator = ' op{i}]);
501 end 545 end
502 end 546 end
503 547