Mercurial > repos > public > sbplib
changeset 1059:d8ab528f10f6 feature/poroelastic
Massive boundary operator cleanup in Elastic2dVariable
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sat, 26 Jan 2019 15:38:58 -0800 |
| parents | 84933722ec0e |
| children | e40899094f20 |
| files | +scheme/Elastic2dVariable.m |
| diffstat | 1 files changed, 218 insertions(+), 183 deletions(-) [+] |
line wrap: on
line diff
diff -r 84933722ec0e -r d8ab528f10f6 +scheme/Elastic2dVariable.m --- a/+scheme/Elastic2dVariable.m Fri Jan 25 16:35:53 2019 -0800 +++ b/+scheme/Elastic2dVariable.m Sat Jan 26 15:38:58 2019 -0800 @@ -15,9 +15,9 @@ order % Order of accuracy for the approximation % Diagonal matrices for varible coefficients - LAMBDA % Variable coefficient, related to dilation - MU % Shear modulus, variable coefficient - RHO, RHOi % Density, variable + LAMBDA % Lame's first parameter, related to dilation + MU % Shear modulus + RHO, RHOi, RHOi_kron % Density D % Total operator D1 % First derivatives @@ -26,22 +26,28 @@ D2_lambda D2_mu - % Traction operators used for BC - T_l, T_r - tau_l, tau_r + % Boundary operators in cell format, used for BC + T_w, T_e, T_s, T_n - H, Hi, H_1D % Inner products - e_l, e_r + % Traction operators + tau_w, tau_e, tau_s, tau_n % Return vector field + tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field + tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field + % Inner products + H, Hi, Hi_kron, H_1D - d1_l, d1_r % Normal derivatives at the boundary - E % E{i}^T picks out component i + % Boundary inner products (for scalar field) + H_w, H_e, H_s, H_n - H_boundary % Boundary inner products + % Boundary restriction operators + e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary + e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary + e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary + e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field - % Kroneckered norms and coefficients - RHOi_kron - Hi_kron + % E{i}^T picks out component i + E % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. theta_R % Borrowing (d1- D1)^2 from R @@ -105,10 +111,10 @@ D2 = cell(dim,1); H = cell(dim,1); Hi = cell(dim,1); - e_l = cell(dim,1); - e_r = cell(dim,1); - d1_l = cell(dim,1); - d1_r = cell(dim,1); + e_0 = cell(dim,1); + e_m = cell(dim,1); + d1_0 = cell(dim,1); + d1_m = cell(dim,1); for i = 1:dim I{i} = speye(m(i)); @@ -116,10 +122,10 @@ D2{i} = ops{i}.D2; H{i} = ops{i}.H; Hi{i} = ops{i}.HI; - e_l{i} = ops{i}.e_l; - e_r{i} = ops{i}.e_r; - d1_l{i} = ops{i}.d1_l; - d1_r{i} = ops{i}.d1_r; + e_0{i} = ops{i}.e_l; + e_m{i} = ops{i}.e_r; + d1_0{i} = ops{i}.d1_l; + d1_m{i} = ops{i}.d1_r; end %====== Assemble full operators ======== @@ -134,25 +140,59 @@ obj.D1 = cell(dim,1); obj.D2_lambda = cell(dim,1); obj.D2_mu = cell(dim,1); - obj.e_l = cell(dim,1); - obj.e_r = cell(dim,1); - obj.d1_l = cell(dim,1); - obj.d1_r = cell(dim,1); % D1 obj.D1{1} = kron(D1{1},I{2}); obj.D1{2} = kron(I{1},D1{2}); - % Boundary operators - obj.e_l{1} = kron(e_l{1},I{2}); - obj.e_l{2} = kron(I{1},e_l{2}); - obj.e_r{1} = kron(e_r{1},I{2}); - obj.e_r{2} = kron(I{1},e_r{2}); + % Boundary restriction operators + e_l = cell(dim,1); + e_r = cell(dim,1); + e_l{1} = kron(e_0{1}, I{2}); + e_l{2} = kron(I{1}, e_0{2}); + e_r{1} = kron(e_m{1}, I{2}); + e_r{2} = kron(I{1}, e_m{2}); + + e_scalar_w = e_l{1}; + e_scalar_e = e_r{1}; + e_scalar_s = e_l{2}; + e_scalar_n = e_r{2}; + + I_dim = speye(dim, dim); + e_w = kron(e_scalar_w, I_dim); + e_e = kron(e_scalar_e, I_dim); + e_s = kron(e_scalar_s, I_dim); + e_n = kron(e_scalar_n, I_dim); - obj.d1_l{1} = kron(d1_l{1},I{2}); - obj.d1_l{2} = kron(I{1},d1_l{2}); - obj.d1_r{1} = kron(d1_r{1},I{2}); - obj.d1_r{2} = kron(I{1},d1_r{2}); + % Boundary derivatives + d1_l = cell(dim,1); + d1_r = cell(dim,1); + d1_l{1} = kron(d1_0{1}, I{2}); + d1_l{2} = kron(I{1}, d1_0{2}); + d1_r{1} = kron(d1_m{1}, I{2}); + d1_r{2} = kron(I{1}, d1_m{2}); + + + % E{i}^T picks out component i. + E = cell(dim,1); + I = speye(m_tot,m_tot); + for i = 1:dim + e = sparse(dim,1); + e(i) = 1; + E{i} = kron(I,e); + end + obj.E = E; + + e1_w = (e_scalar_w'*E{1}')'; + e1_e = (e_scalar_e'*E{1}')'; + e1_s = (e_scalar_s'*E{1}')'; + e1_n = (e_scalar_n'*E{1}')'; + + e2_w = (e_scalar_w'*E{2}')'; + e2_e = (e_scalar_e'*E{2}')'; + e2_s = (e_scalar_s'*E{2}')'; + e2_n = (e_scalar_n'*E{2}')'; + % D2 for i = 1:dim @@ -182,21 +222,12 @@ % Quadratures obj.H = kron(H{1},H{2}); obj.Hi = inv(obj.H); - obj.H_boundary = cell(dim,1); - obj.H_boundary{1} = H{2}; - obj.H_boundary{2} = H{1}; + obj.H_w = H{2}; + obj.H_e = H{2}; + obj.H_s = H{1}; + obj.H_n = H{1}; obj.H_1D = {H{1}, H{2}}; - % E{i}^T picks out component i. - E = cell(dim,1); - I = speye(m_tot,m_tot); - for i = 1:dim - e = sparse(dim,1); - e(i) = 1; - E{i} = kron(I,e); - end - obj.E = E; - % Differentiation matrix D (without SAT) D2_lambda = obj.D2_lambda; D2_mu = obj.D2_mu; @@ -221,16 +252,14 @@ % Numerical traction operators for BC. % Because d1 =/= e0^T*D1, the numerical tractions are different % at every boundary. + % + % Formula at boundary j: % tau^{j}_i = sum_k T^{j}_{ik} u_k + % T_l = cell(dim,1); T_r = cell(dim,1); tau_l = cell(dim,1); tau_r = cell(dim,1); - % tau^{j}_i = sum_k T^{j}_{ik} u_k - d1_l = obj.d1_l; - d1_r = obj.d1_r; - e_l = obj.e_l; - e_r = obj.e_r; D1 = obj.D1; % Loop over boundaries @@ -250,45 +279,72 @@ % Loop over components for i = 1:dim - tau_l{j}{i} = sparse(n_l, dim*m_tot); - tau_r{j}{i} = sparse(n_r, dim*m_tot); + tau_l{j}{i} = sparse(dim*m_tot, n_l); + tau_r{j}{i} = sparse(dim*m_tot, n_r); for k = 1:dim T_l{j}{i,k} = ... - -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})... - -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})... - -d(i,k)*MU_l*d1_l{j}'; + (-d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})... + -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})... + -d(i,k)*MU_l*d1_l{j}')'; T_r{j}{i,k} = ... - d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})... + (d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})... +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})... - +d(i,k)*MU_r*d1_r{j}'; + +d(i,k)*MU_r*d1_r{j}')'; - tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; - tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; + tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,k}'*E{k}')'; + tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,k}'*E{k}')'; end end end - % Transpose T and tau to match boundary operator convention - for i = 1:dim - for j = 1:dim - tau_l{i}{j} = transpose(tau_l{i}{j}); - tau_r{i}{j} = transpose(tau_r{i}{j}); - for k = 1:dim - T_l{i}{j,k} = transpose(T_l{i}{j,k}); - T_r{i}{j,k} = transpose(T_r{i}{j,k}); - end - end - end + % Traction tensors, T_ij + obj.T_w = T_l{1}; + obj.T_e = T_r{1}; + obj.T_s = T_l{2}; + obj.T_n = T_r{2}; + + % Restriction operators + obj.e_w = e_w; + obj.e_e = e_e; + obj.e_s = e_s; + obj.e_n = e_n; + + obj.e1_w = e1_w; + obj.e1_e = e1_e; + obj.e1_s = e1_s; + obj.e1_n = e1_n; + + obj.e2_w = e2_w; + obj.e2_e = e2_e; + obj.e2_s = e2_s; + obj.e2_n = e2_n; - obj.T_l = T_l; - obj.T_r = T_r; - obj.tau_l = tau_l; - obj.tau_r = tau_r; + obj.e_scalar_w = e_scalar_w; + obj.e_scalar_e = e_scalar_e; + obj.e_scalar_s = e_scalar_s; + obj.e_scalar_n = e_scalar_n; + + % First component of traction + obj.tau1_w = tau_l{1}{1}; + obj.tau1_e = tau_r{1}{1}; + obj.tau1_s = tau_l{2}{1}; + obj.tau1_n = tau_r{2}{1}; + + % Second component of traction + obj.tau2_w = tau_l{1}{2}; + obj.tau2_e = tau_r{1}{2}; + obj.tau2_s = tau_l{2}{2}; + obj.tau2_n = tau_r{2}{2}; + + % Traction vectors + obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')'; + obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')'; + obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')'; + obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')'; % Kroneckered norms and coefficients - I_dim = speye(dim); obj.RHOi_kron = kron(obj.RHOi, I_dim); obj.Hi_kron = kron(obj.Hi, I_dim); @@ -360,10 +416,11 @@ comp = obj.getComponent(comp, boundary); end - e = obj.getBoundaryOperator('e', boundary); - T = obj.getBoundaryOperator('T', boundary); - tau = obj.getBoundaryOperator('tau', boundary); - H_gamma = obj.getBoundaryOperator('H', boundary); + e = obj.getBoundaryOperatorForScalarField('e', boundary); + tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary); + T = obj.getBoundaryTractionOperator(boundary); + alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary); + H_gamma = obj.getBoundaryQuadratureForScalarField(boundary); E = obj.E; Hi = obj.Hi; @@ -375,7 +432,7 @@ m_tot = obj.grid.N(); % Preallocate - [~, col] = size(tau{comp}); + [~, col] = size(tau); closure = sparse(dim*m_tot, dim*m_tot); penalty = sparse(dim*m_tot, col); @@ -385,8 +442,6 @@ % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - alpha = obj.getBoundaryOperator('alpha', boundary); - % Loop over components that Dirichlet penalties end up on for i = 1:dim C = transpose(T{k,i}); @@ -398,7 +453,7 @@ % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}'; + closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau'; penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; % Unknown boundary condition @@ -435,11 +490,11 @@ % Operators without subscripts are from the own domain. % Get boundary operators - e = obj.getBoundaryOperator('e_tot', boundary); - tau = obj.getBoundaryOperator('tau_tot', boundary); + e = obj.getBoundaryOperator('e', boundary); + tau = obj.getBoundaryOperator('tau', boundary); - e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary); - tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary); + e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); + tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary); H_gamma = obj.getBoundaryQuadrature(boundary); @@ -448,8 +503,8 @@ RHOi = obj.RHOi_kron; % Penalty strength operators - alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary); - alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary); + alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary); + alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary); closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e'); penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v'); @@ -492,83 +547,62 @@ % Returns the boundary operator op for the boundary specified by the string boundary. % op -- string - % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() function o = getBoundaryOperator(obj, op, boundary) assertIsMember(boundary, {'w', 'e', 's', 'n'}) - assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'}) - - switch boundary - case {'w', 'e'} - j = 1; - case {'s', 'n'} - j = 2; - end + assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha'}) switch op - case 'e' - switch boundary - case {'w', 's'} - o = obj.e_l{j}; - case {'e', 'n'} - o = obj.e_r{j}; - end - case 'e_tot' - e = obj.getBoundaryOperator('e', boundary); - I_dim = speye(obj.dim, obj.dim); - o = kron(e, I_dim); - - case 'd' - switch boundary - case {'w', 's'} - o = obj.d1_l{j}; - case {'e', 'n'} - o = obj.d1_r{j}; - end - - case 'T' - switch boundary - case {'w', 's'} - o = obj.T_l{j}; - case {'e', 'n'} - o = obj.T_r{j}; - end - - case 'tau' - switch boundary - case {'w', 's'} - o = obj.tau_l{j}; - case {'e', 'n'} - o = obj.tau_r{j}; - end - - case 'tau_tot' - e = obj.getBoundaryOperator('e', boundary); - tau = obj.getBoundaryOperator('tau', boundary); - - I_dim = speye(obj.dim, obj.dim); - e_tot = kron(e, I_dim); - E = obj.E; - tau_tot = (e_tot'*E{1}*e*tau{1}')'; - for i = 2:obj.dim - tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')'; - end - o = tau_tot; - - case 'H' - o = obj.H_boundary{j}; + case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'} + o = obj.([op, '_', boundary]); case 'alpha' % alpha = alpha(i,j) is the penalty strength for displacement BC. - e = obj.getBoundaryOperator('e', boundary); + e = obj.getBoundaryOperatorForScalarField('e', boundary); + e_tot = obj.getBoundaryOperator('e', boundary); + alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary); + E = obj.E; + [m, n] = size(alpha{1,1}); + alpha_tot = sparse(m*obj.dim, n*obj.dim); + for i = 1:obj.dim + for l = 1:obj.dim + alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; + end + end + o = alpha_tot; + end + + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + function o = getBoundaryOperatorForScalarField(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e', 'alpha'}) + + switch op + + case 'e' + o = obj.(['e_scalar', '_', boundary]); + + case 'alpha' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + e = obj.getBoundaryOperatorForScalarField('e', boundary); LAMBDA = obj.LAMBDA; MU = obj.MU; + dim = obj.dim; - dim = obj.dim; - theta_R = obj.theta_R{j}; - theta_H = obj.theta_H{j}; - theta_M = obj.theta_M{j}; + switch boundary + case {'w', 'e'} + k = 1; + case {'s', 'n'} + k = 2; + end + + theta_R = obj.theta_R{k}; + theta_H = obj.theta_H{k}; + theta_M = obj.theta_M{k}; a_lambda = dim/theta_H + 1/theta_R; a_mu_i = 2/theta_M; @@ -582,45 +616,46 @@ + d(i,j)* a_mu_i*MU ... + db(i,j)*a_mu_ij*MU; for i = 1:obj.dim - for l = 1:obj.dim - alpha{i,l} = d(i,l)*alpha_func(i,j)*e; + for j = 1:obj.dim + alpha{i,j} = d(i,j)*alpha_func(i,k)*e; end end o = alpha; - - case 'alpha_tot' - % alpha = alpha(i,j) is the penalty strength for displacement BC. - e = obj.getBoundaryOperator('e', boundary); - e_tot = obj.getBoundaryOperator('e_tot', boundary); - alpha = obj.getBoundaryOperator('alpha', boundary); - E = obj.E; - [m, n] = size(alpha{1,1}); - alpha_tot = sparse(m*obj.dim, n*obj.dim); - for i = 1:obj.dim - for l = 1:obj.dim - alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; - end - end - o = alpha_tot; end end + % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary. + % op -- string + function T = getBoundaryTractionOperator(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + T = obj.(['T', '_', boundary]); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary unknowns + % + % boundary -- string function H = getBoundaryQuadrature(obj, boundary) - switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} - j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} - j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); - end - H = obj.H_boundary{j}; + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H = obj.getBoundaryQuadratureForScalarField(boundary); I_dim = speye(obj.dim, obj.dim); H = kron(H, I_dim); end + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary grid points + % + % boundary -- string + function H_b = getBoundaryQuadratureForScalarField(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H_b = obj.(['H_', boundary]); + end + function N = size(obj) N = obj.dim*prod(obj.m); end