Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 882:14fee299ada2 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
In Elastic2dVariable: Improve notation for borrowing constants. Update interface constant to match corrected derivation. Tests ok.
author Martin Almquist <malmquist@stanford.edu>
date Fri, 26 Oct 2018 16:35:23 -0700
parents 607c631f175e
children e30aaa4a3e09 386ef449df51
comparison
equal deleted inserted replaced
881:67228a10dfad 882:14fee299ada2
29 % Traction operators used for BC 29 % Traction operators used for BC
30 T_l, T_r 30 T_l, T_r
31 tau_l, tau_r 31 tau_l, tau_r
32 32
33 H, Hi, H_1D % Inner products 33 H, Hi, H_1D % Inner products
34 phi % Borrowing constant for (d1 - e^T*D1) from R
35 gamma % Borrowing constant for d1 from M
36 H11 % First element of H
37 e_l, e_r 34 e_l, e_r
38 d1_l, d1_r % Normal derivatives at the boundary 35 d1_l, d1_r % Normal derivatives at the boundary
39 E % E{i}^T picks out component i 36 E % E{i}^T picks out component i
40 37
41 H_boundary % Boundary inner products 38 H_boundary % Boundary inner products
42 39
43 % Kroneckered norms and coefficients 40 % Kroneckered norms and coefficients
44 RHOi_kron 41 RHOi_kron
45 Hi_kron 42 Hi_kron
43
44 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
45 theta_R % Borrowing (d1- D1)^2 from R
46 theta_H % First entry in norm matrix
47 theta_M % Borrowing d1^2 from M.
46 48
47 % Structures used for adjoint optimization 49 % Structures used for adjoint optimization
48 B 50 B
49 end 51 end
50 52
89 ops{i} = opSet{i}(m(i), lim{i}, order); 91 ops{i} = opSet{i}(m(i), lim{i}, order);
90 end 92 end
91 93
92 % Borrowing constants 94 % Borrowing constants
93 for i = 1:dim 95 for i = 1:dim
94 beta = ops{i}.borrowing.R.delta_D; 96 obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D;
95 obj.H11{i} = ops{i}.borrowing.H11; 97 obj.theta_H{i} = h(i)*ops{i}.borrowing.H11;
96 obj.phi{i} = beta/obj.H11{i}; 98 obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1;
97 obj.gamma{i} = ops{i}.borrowing.M.d1;
98 end 99 end
99 100
100 I = cell(dim,1); 101 I = cell(dim,1);
101 D1 = cell(dim,1); 102 D1 = cell(dim,1);
102 D2 = cell(dim,1); 103 D2 = cell(dim,1);
325 switch type 326 switch type
326 327
327 % Dirichlet boundary condition 328 % Dirichlet boundary condition
328 case {'D','d','dirichlet','Dirichlet'} 329 case {'D','d','dirichlet','Dirichlet'}
329 330
330 phi = obj.phi{j}; 331 theta_R = obj.theta_R{j};
331 h = obj.h(j); 332 theta_H = obj.theta_H{j};
332 h11 = obj.H11{j}*h; 333 theta_M = obj.theta_M{j};
333 gamma = obj.gamma{j}; 334
334 335 a_lambda = dim/theta_H + 1/theta_R;
335 a_lambda = dim/h11 + 1/(h11*phi); 336 a_mu_i = 2/theta_M;
336 a_mu_i = 2/(gamma*h); 337 a_mu_ij = 2/theta_H + 1/theta_R;
337 a_mu_ij = 2/h11 + 1/(h11*phi);
338 338
339 d = @kroneckerDelta; % Kronecker delta 339 d = @kroneckerDelta; % Kronecker delta
340 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 340 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
341 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 341 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
342 + d(i,j)* a_mu_i*MU ... 342 + d(i,j)* a_mu_i*MU ...
396 lambda_v = e_v'*LAMBDA_v*e_v; 396 lambda_v = e_v'*LAMBDA_v*e_v;
397 mu_v = e_v'*MU_v*e_v; 397 mu_v = e_v'*MU_v*e_v;
398 %------------------------- 398 %-------------------------
399 399
400 % Borrowing constants 400 % Borrowing constants
401 phi_u = obj.phi{j}; 401 theta_R_u = obj.theta_R{j};
402 h_u = obj.h(j); 402 theta_H_u = obj.theta_H{j};
403 h11_u = obj.H11{j}*h_u; 403 theta_M_u = obj.theta_M{j};
404 gamma_u = obj.gamma{j}; 404
405 405 theta_R_v = neighbour_scheme.theta_R{j_v};
406 phi_v = neighbour_scheme.phi{j_v}; 406 theta_H_v = neighbour_scheme.theta_H{j_v};
407 h_v = neighbour_scheme.h(j_v); 407 theta_M_v = neighbour_scheme.theta_M{j_v};
408 h11_v = neighbour_scheme.H11{j_v}*h_v; 408
409 gamma_v = neighbour_scheme.gamma{j_v}; 409 function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu)
410 410 alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M);
411 % E > sum_i 1/(2*alpha_ij)*(tau_i)^2 411 alpha_ij = mu/(2*th_H) + mu/(4*th_R);
412 function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) 412 end
413 th1 = h11/(2*dim); 413
414 th2 = h11*phi/2; 414 [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u);
415 th3 = h*gamma; 415 [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v);
416 a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3); 416 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v);
417 a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3); 417 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v);
418 alpha_ii = a1 + sqrt(a2 + a1^2);
419
420 alpha_ij = mu*(2/h11 + 1/(phi*h11));
421 end
422
423 [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);
424 [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);
425 sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;
426 sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;
427 418
428 d = @kroneckerDelta; % Kronecker delta 419 d = @kroneckerDelta; % Kronecker delta
429 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 420 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
430 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij); 421 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
431 422
517 varargout{i} = obj.tau_l{j}; 508 varargout{i} = obj.tau_l{j};
518 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 509 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
519 varargout{i} = obj.tau_r{j}; 510 varargout{i} = obj.tau_r{j};
520 end 511 end
521 case 'alpha' 512 case 'alpha'
522 % alpha = alpha(i,j) is the penalty strength for displacement BC. 513 % alpha = alpha(i,j) is the penalty strength for displacement BC.
523 tuning = 1.2; 514 tuning = 1.2;
524 LAMBDA = obj.LAMBDA; 515 LAMBDA = obj.LAMBDA;
525 MU = obj.MU; 516 MU = obj.MU;
526 517
527 phi = obj.phi{j}; 518 phi = obj.phi{j};