Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 729:aa8cf3851de8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update multiblock.DiffOp to work for systems.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 20 Apr 2018 16:56:49 -0700
parents 6d5953fc090e
children 9f28cf266f86
comparison
equal deleted inserted replaced
728:0aff87f6fb2c 729:aa8cf3851de8
62 m = g.size(); 62 m = g.size();
63 m_tot = g.N(); 63 m_tot = g.N();
64 64
65 h = g.scaling(); 65 h = g.scaling();
66 lim = g.lim; 66 lim = g.lim;
67 if isempty(lim)
68 x = g.x;
69 lim = cell(length(x),1);
70 for i = 1:length(x)
71 lim{i} = {min(x{i}), max(x{i})};
72 end
73 end
67 74
68 % 1D operators 75 % 1D operators
69 ops = cell(dim,1); 76 ops = cell(dim,1);
70 for i = 1:dim 77 for i = 1:dim
71 ops{i} = opSet{i}(m(i), lim{i}, order); 78 ops{i} = opSet{i}(m(i), lim{i}, order);
267 % neighbour_scheme is an instance of Scheme that should be interfaced to. 274 % neighbour_scheme is an instance of Scheme that should be interfaced to.
268 % neighbour_boundary is a string specifying which boundary to interface to. 275 % neighbour_boundary is a string specifying which boundary to interface to.
269 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 276 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
270 default_arg('type',{'free','free'}); 277 default_arg('type',{'free','free'});
271 default_arg('parameter', []); 278 default_arg('parameter', []);
279
280 if ~iscell(type)
281 type = {type, type};
282 end
272 283
273 % j is the coordinate direction of the boundary 284 % j is the coordinate direction of the boundary
274 j = obj.get_boundary_number(boundary); 285 j = obj.get_boundary_number(boundary);
275 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 286 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
276 287
405 closure = sparse(dim*m_tot_u, dim*m_tot_u); 416 closure = sparse(dim*m_tot_u, dim*m_tot_u);
406 penalty = sparse(dim*m_tot_u, dim*m_tot_v); 417 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
407 418
408 % Loop over components that penalties end up on 419 % Loop over components that penalties end up on
409 for i = 1:dim 420 for i = 1:dim
410 closure = closure - sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e'*E{i}'; 421 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
411 penalty = penalty + sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e_v'*E_v{i}'; 422 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
412 423
413 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}'; 424 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
414 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}'; 425 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
415 426
416 % Loop over components that we have interface conditions on 427 % Loop over components that we have interface conditions on
417 for k = 1:dim 428 for k = 1:dim
418 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}'; 429 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
419 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}'; 430 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';