Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 1118:07d0caf915e4 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Introduce optFlag so that one can choose not to build the optimization-related cell of matrices called B. It is too computationally costly and should probably be done in a different way.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 05 May 2019 19:05:31 -0700
parents e40899094f20
children d9da4c1cdaa0
comparison
equal deleted inserted replaced
1117:27019aca2f17 1118:07d0caf915e4
59 end 59 end
60 60
61 methods 61 methods
62 62
63 % The coefficients can either be function handles or grid functions 63 % The coefficients can either be function handles or grid functions
64 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet) 64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
65 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet, optFlag)
65 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); 66 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
66 default_arg('lambda', @(x,y) 0*x+1); 67 default_arg('lambda', @(x,y) 0*x+1);
67 default_arg('mu', @(x,y) 0*x+1); 68 default_arg('mu', @(x,y) 0*x+1);
68 default_arg('rho', @(x,y) 0*x+1); 69 default_arg('rho', @(x,y) 0*x+1);
70 default_arg('optFlag', false);
69 dim = 2; 71 dim = 2;
70 72
71 assert(isa(g, 'grid.Cartesian')) 73 assert(isa(g, 'grid.Cartesian'))
72 74
73 if isa(lambda, 'function_handle') 75 if isa(lambda, 'function_handle')
354 obj.order = order; 356 obj.order = order;
355 obj.grid = g; 357 obj.grid = g;
356 obj.dim = dim; 358 obj.dim = dim;
357 359
358 % B, used for adjoint optimization 360 % B, used for adjoint optimization
359 B = cell(dim, 1); 361 B = [];
360 for i = 1:dim 362 if optFlag
361 B{i} = cell(m_tot, 1); 363 B = cell(dim, 1);
362 end 364 for i = 1:dim
363 365 B{i} = cell(m_tot, 1);
364 for i = 1:dim 366 end
365 for j = 1:m_tot 367
366 B{i}{j} = sparse(m_tot, m_tot); 368 B0 = sparse(m_tot, m_tot);
367 end 369 for i = 1:dim
368 end 370 for j = 1:m_tot
369 371 B{i}{j} = B0;
370 ind = grid.funcToMatrix(g, 1:m_tot); 372 end
371 373 end
372 % Direction 1 374
373 for k = 1:m(1) 375 ind = grid.funcToMatrix(g, 1:m_tot);
374 c = sparse(m(1),1); 376
375 c(k) = 1; 377 % Direction 1
376 [~, B_1D] = ops{1}.D2(c); 378 for k = 1:m(1)
377 for l = 1:m(2) 379 c = sparse(m(1),1);
378 p = ind(:,l); 380 c(k) = 1;
379 B{1}{(k-1)*m(2) + l}(p, p) = B_1D; 381 [~, B_1D] = ops{1}.D2(c);
380 end 382 for l = 1:m(2)
381 end 383 p = ind(:,l);
382 384 B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
383 % Direction 2 385 end
384 for k = 1:m(2) 386 end
385 c = sparse(m(2),1); 387
386 c(k) = 1; 388 % Direction 2
387 [~, B_1D] = ops{2}.D2(c); 389 for k = 1:m(2)
388 for l = 1:m(1) 390 c = sparse(m(2),1);
389 p = ind(l,:); 391 c(k) = 1;
390 B{2}{(l-1)*m(2) + k}(p, p) = B_1D; 392 [~, B_1D] = ops{2}.D2(c);
391 end 393 for l = 1:m(1)
392 end 394 p = ind(l,:);
393 395 B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
396 end
397 end
398 end
394 obj.B = B; 399 obj.B = B;
400
395 401
396 end 402 end
397 403
398 404
399 % Closure functions return the operators applied to the own domain to close the boundary 405 % Closure functions return the operators applied to the own domain to close the boundary