Mercurial > repos > public > sbplib
changeset 882:14fee299ada2 feature/poroelastic
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 | 67228a10dfad |
| children | 76efb6a7b466 |
| files | +scheme/Elastic2dVariable.m |
| diffstat | 1 files changed, 28 insertions(+), 37 deletions(-) [+] |
line wrap: on
line diff
diff -r 67228a10dfad -r 14fee299ada2 +scheme/Elastic2dVariable.m --- a/+scheme/Elastic2dVariable.m Fri Oct 26 15:41:26 2018 -0700 +++ b/+scheme/Elastic2dVariable.m Fri Oct 26 16:35:23 2018 -0700 @@ -31,9 +31,6 @@ tau_l, tau_r H, Hi, H_1D % Inner products - phi % Borrowing constant for (d1 - e^T*D1) from R - gamma % Borrowing constant for d1 from M - H11 % First element of H e_l, e_r d1_l, d1_r % Normal derivatives at the boundary E % E{i}^T picks out component i @@ -44,6 +41,11 @@ RHOi_kron Hi_kron + % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. + theta_R % Borrowing (d1- D1)^2 from R + theta_H % First entry in norm matrix + theta_M % Borrowing d1^2 from M. + % Structures used for adjoint optimization B end @@ -91,10 +93,9 @@ % Borrowing constants for i = 1:dim - beta = ops{i}.borrowing.R.delta_D; - obj.H11{i} = ops{i}.borrowing.H11; - obj.phi{i} = beta/obj.H11{i}; - obj.gamma{i} = ops{i}.borrowing.M.d1; + obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D; + obj.theta_H{i} = h(i)*ops{i}.borrowing.H11; + obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1; end I = cell(dim,1); @@ -327,14 +328,13 @@ % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - phi = obj.phi{j}; - h = obj.h(j); - h11 = obj.H11{j}*h; - gamma = obj.gamma{j}; + theta_R = obj.theta_R{j}; + theta_H = obj.theta_H{j}; + theta_M = obj.theta_M{j}; - a_lambda = dim/h11 + 1/(h11*phi); - a_mu_i = 2/(gamma*h); - a_mu_ij = 2/h11 + 1/(h11*phi); + a_lambda = dim/theta_H + 1/theta_R; + a_mu_i = 2/theta_M; + a_mu_ij = 2/theta_H + 1/theta_R; d = @kroneckerDelta; % Kronecker delta db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta @@ -398,32 +398,23 @@ %------------------------- % Borrowing constants - phi_u = obj.phi{j}; - h_u = obj.h(j); - h11_u = obj.H11{j}*h_u; - gamma_u = obj.gamma{j}; - - phi_v = neighbour_scheme.phi{j_v}; - h_v = neighbour_scheme.h(j_v); - h11_v = neighbour_scheme.H11{j_v}*h_v; - gamma_v = neighbour_scheme.gamma{j_v}; + theta_R_u = obj.theta_R{j}; + theta_H_u = obj.theta_H{j}; + theta_M_u = obj.theta_M{j}; - % E > sum_i 1/(2*alpha_ij)*(tau_i)^2 - function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) - th1 = h11/(2*dim); - th2 = h11*phi/2; - th3 = h*gamma; - a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3); - a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3); - alpha_ii = a1 + sqrt(a2 + a1^2); + theta_R_v = neighbour_scheme.theta_R{j_v}; + theta_H_v = neighbour_scheme.theta_H{j_v}; + theta_M_v = neighbour_scheme.theta_M{j_v}; - alpha_ij = mu*(2/h11 + 1/(phi*h11)); + function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu) + alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M); + alpha_ij = mu/(2*th_H) + mu/(4*th_R); end - [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u); - [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v); - sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4; - sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4; + [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u); + [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v); + sigma_ii = tuning*(alpha_ii_u + alpha_ii_v); + sigma_ij = tuning*(alpha_ij_u + alpha_ij_v); d = @kroneckerDelta; % Kronecker delta db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta @@ -519,7 +510,7 @@ varargout{i} = obj.tau_r{j}; end case 'alpha' - % alpha = alpha(i,j) is the penalty strength for displacement BC. + % alpha = alpha(i,j) is the penalty strength for displacement BC. tuning = 1.2; LAMBDA = obj.LAMBDA; MU = obj.MU;