Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 727:6d5953fc090e feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First implementation of elastic interface
author Martin Almquist <malmquist@stanford.edu>
date Fri, 20 Apr 2018 11:32:04 -0700
parents 37d5d69b1a0d
children aa8cf3851de8
comparison
equal deleted inserted replaced
726:37d5d69b1a0d 727:6d5953fc090e
270 default_arg('type',{'free','free'}); 270 default_arg('type',{'free','free'});
271 default_arg('parameter', []); 271 default_arg('parameter', []);
272 272
273 % j is the coordinate direction of the boundary 273 % j is the coordinate direction of the boundary
274 j = obj.get_boundary_number(boundary); 274 j = obj.get_boundary_number(boundary);
275 [e, T, tau] = obj.get_boundary_operator({'e','T','tau'}, boundary); 275 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
276 276
277 E = obj.E; 277 E = obj.E;
278 Hi = obj.Hi; 278 Hi = obj.Hi;
279 H_gamma = obj.H_boundary{j};
280 LAMBDA = obj.LAMBDA; 279 LAMBDA = obj.LAMBDA;
281 MU = obj.MU; 280 MU = obj.MU;
282 RHOi = obj.RHOi; 281 RHOi = obj.RHOi;
283 282
284 dim = obj.dim; 283 dim = obj.dim;
285 m_tot = obj.grid.N(); 284 m_tot = obj.grid.N();
286
287 RHOi_kron = obj.RHOi_kron;
288 Hi_kron = obj.Hi_kron;
289 285
290 % Preallocate 286 % Preallocate
291 closure = sparse(dim*m_tot, dim*m_tot); 287 closure = sparse(dim*m_tot, dim*m_tot);
292 penalty = cell(dim,1); 288 penalty = cell(dim,1);
293 for k = 1:dim 289 for k = 1:dim
339 end 335 end
340 336
341 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 337 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
342 % u denotes the solution in the own domain 338 % u denotes the solution in the own domain
343 % v denotes the solution in the neighbour domain 339 % v denotes the solution in the neighbour domain
340 % Operators without subscripts are from the own domain.
344 tuning = 1.2; 341 tuning = 1.2;
345 % tuning = 20.2; 342
346 error('Interface not implemented'); 343 % j is the coordinate direction of the boundary
344 j = obj.get_boundary_number(boundary);
345 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
346
347 % Get boundary operators
348 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
349 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
350
351 % Operators and quantities that correspond to the own domain only
352 Hi = obj.Hi;
353 RHOi = obj.RHOi;
354 dim = obj.dim;
355
356 %--- Other operators ----
357 m_tot_u = obj.grid.N();
358 E = obj.E;
359 LAMBDA_u = obj.LAMBDA;
360 MU_u = obj.MU;
361 lambda_u = e'*LAMBDA_u*e;
362 mu_u = e'*MU_u*e;
363
364 m_tot_v = neighbour_scheme.grid.N();
365 E_v = neighbour_scheme.E;
366 LAMBDA_v = neighbour_scheme.LAMBDA;
367 MU_v = neighbour_scheme.MU;
368 lambda_v = e_v'*LAMBDA_v*e_v;
369 mu_v = e_v'*MU_v*e_v;
370 %-------------------------
371
372 % Borrowing constants
373 phi_u = obj.phi{j};
374 h_u = obj.h(j);
375 h11_u = obj.H11{j}*h_u;
376 gamma_u = obj.gamma{j};
377
378 phi_v = neighbour_scheme.phi{j_v};
379 h_v = neighbour_scheme.h(j_v);
380 h11_v = neighbour_scheme.H11{j_v}*h_v;
381 gamma_v = neighbour_scheme.gamma{j_v};
382
383 % E > sum_i 1/(2*alpha_ij)*(tau_i)^2
384 function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu)
385 th1 = h11/(2*dim);
386 th2 = h11*phi/2;
387 th3 = h*gamma;
388 a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3);
389 a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3);
390 alpha_ii = a1 + sqrt(a2 + a1^2);
391
392 alpha_ij = 2/h11 + 1/(phi*h11);
393 end
394
395 [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);
396 [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);
397 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;
398 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;
399
400 d = @kroneckerDelta; % Kronecker delta
401 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
402 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
403
404 % Preallocate
405 closure = sparse(dim*m_tot_u, dim*m_tot_u);
406 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
407
408 % Loop over components that penalties end up on
409 for i = 1:dim
410 closure = closure - sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e'*E{i}';
411 penalty = penalty + sigma(i,j)*E{i}*RHOi*Hi*e*H_gamma*e_v'*E_v{i}';
412
413 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}';
415
416 % Loop over components that we have interface conditions on
417 for k = 1:dim
418 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}';
420 end
421 end
347 end 422 end
348 423
349 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 424 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
350 function [j, nj] = get_boundary_number(obj, boundary) 425 function [j, nj] = get_boundary_number(obj, boundary)
351 426