Mercurial > repos > public > sbplib
changeset 1346:464e6f65c2c6 feature/poroelastic
Refactor computation of borrowing from C. Update to correct number of boundary points for borrowing from R.
| author | Martin Almquist <martin.almquist@it.uu.se> |
|---|---|
| date | Tue, 30 Apr 2024 14:18:47 +0200 |
| parents | 14f44e81e1e3 |
| children | ac54767ae1fb |
| files | +scheme/Elastic2dVariableAnisotropicIncompatible.m |
| diffstat | 1 files changed, 67 insertions(+), 56 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Elastic2dVariableAnisotropicIncompatible.m Tue Apr 30 13:33:39 2024 +0200 +++ b/+scheme/Elastic2dVariableAnisotropicIncompatible.m Tue Apr 30 14:18:47 2024 +0200 @@ -50,7 +50,7 @@ h11 % First entry in norm matrix theta_R % Borrowing from R - nBP % Number of boundary points, related to borrowing. + nBP_R % Number of boundary points in R, related to borrowing from R. e_shifted_scalar_w, e_shifted_scalar_e, e_shifted_scalar_s, e_shifted_scalar_n; % Act on scalar field, return scalar field end @@ -124,12 +124,15 @@ case 2 width = 3; nBP = 2; + nBP_R = 2; case 4 width = 5; nBP = 6; + nBP_R = 4; case 6 width = 7; nBP = 9; + nBP_R = 7; end I = cell(dim,1); @@ -187,11 +190,11 @@ e_n = kron(e_scalar_n, I_dim); % Boundary restriction, shifted up to nBP-1 points - e_shifted_scalar_w = cell(nBP, 1); - e_shifted_scalar_e = cell(nBP, 1); - e_shifted_scalar_s = cell(nBP, 1); - e_shifted_scalar_n = cell(nBP, 1); - for i = 1:nBP-1 + e_shifted_scalar_w = cell(nBP_R, 1); + e_shifted_scalar_e = cell(nBP_R, 1); + e_shifted_scalar_s = cell(nBP_R, 1); + e_shifted_scalar_n = cell(nBP_R, 1); + for i = 1:nBP_R-1 e_0_nBP = {0*e_0{1}, 0*e_0{2}}; e_m_nBP = {0*e_m{1}, 0*e_m{2}}; e_0_nBP{1}(i) = 1; @@ -407,7 +410,7 @@ obj.order = order; obj.grid = g; obj.dim = dim; - obj.nBP = nBP; + obj.nBP_R = nBP_R; end @@ -429,7 +432,7 @@ % dComps = vector of components with displacement BC. Default: 1:dim. % In this way, we can specify one BC at a time even though the SATs depend on all BC. function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) - default_arg('tuning', 1.0); + default_arg('tuning', 1.2); assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); comp = bc{1}; @@ -439,7 +442,6 @@ end e = obj.getBoundaryOperatorForScalarField('e', boundary); - e_sh = obj.getBoundaryOperatorForScalarField('e_shifted', boundary); tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary); T = obj.getBoundaryTractionOperator(boundary); [h11, th_R] = obj.getBorrowing(boundary); @@ -487,54 +489,11 @@ % Symmetric part only required on components with displacement BC. % (Otherwise it's not symmetric.) - - % 1. Compute beta such that C_bJbL_boundary - \beta*C_bJbL_i >= 0, i = 0, 1, ..., nBP. - [~, N] = size(e); - beta_vec = ones(N, 1); - b = obj.getComponent('normal', boundary); - C_boundary = cell(dim, dim, dim, dim); - for i = 1:dim - for j = 1:dim - for k = 1:dim - for l = 1:dim - C_boundary{i,j,k,l} = e'*C{i,j,k,l}*e ; - end - end - end - end - - % Loop over shift inward from boundary - for m = 1:obj.nBP-1 - C_shifted = cell(dim, dim, dim, dim); - for i = 1:dim - for j = 1:dim - for k = 1:dim - for l = 1:dim - C_shifted{i,j,k,l} = e_sh{m}'*C{i,j,k,l}*e_sh{m}; - end - end - end - end - - % Loop along boundary - for i = 1:N - C_boundary_mat = [C_boundary{b,1,b,1}(i,i), C_boundary{b,1,b,2}(i,i); ... - C_boundary{b,2,b,1}(i,i), C_boundary{b,2,b,2}(i,i)]; - - C_shifted_mat = [C_shifted{b,1,b,1}(i,i), C_shifted{b,1,b,2}(i,i); ... - C_shifted{b,2,b,1}(i,i), C_shifted{b,2,b,2}(i,i)]; - - beta = obj.computeBorrowingFromC(C_boundary_mat, C_shifted_mat); - beta_vec(i) = min(beta_vec(i), beta); - end - end - if min(beta_vec) < 1e-10 - disp('WARNING! No borrowing from C seems to be possible.') - disp(['min(beta) = ' num2str(min(beta_vec))]) - end + beta_vec = obj.computeBorrowingFromC(boundary); beta_inv = e*spdiag(1./beta_vec)*e'; % Compensate at corners + [~, N] = size(e); corner_weight = ones(N, 1); corner_weight(1) = 2; corner_weight(end) = 2; @@ -550,7 +509,7 @@ Z = Z + nu(l)*C{l,i,k,j}*nu(k); end end - Z = -1.2*tuning*(corner_weight/h11 + beta_inv/th_R)*Z; + Z = -tuning*(corner_weight/h11 + beta_inv/th_R)*Z; X = Z; closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; @@ -705,7 +664,7 @@ end % Computes the largest beta such that C_b - beta*C_s >= 0, using bisection. - function beta = computeBorrowingFromC(obj, C_b, C_s) + function beta = borrowingBisection(obj, C_b, C_s) tol = 1e-8; beta_min = 0; beta_max = 1; @@ -725,6 +684,58 @@ end end + % Computes largest beta such that C_bJbL_boundary - \beta*C_bJbL_i >= 0, i = 0, 1, ..., nBP. + function beta = computeBorrowingFromC(obj, boundary) + + e = obj.getBoundaryOperatorForScalarField('e', boundary); + e_sh = obj.getBoundaryOperatorForScalarField('e_shifted', boundary); + dim = obj.dim; + + [~, N] = size(e); + beta = ones(N, 1); + b = obj.getComponent('normal', boundary); + C_boundary = cell(dim, dim, dim, dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + C_boundary{i,j,k,l} = e'*obj.C{i,j,k,l}*e ; + end + end + end + end + + % Loop over shift inward from boundary + for m = 1:obj.nBP_R-1 + C_shifted = cell(dim, dim, dim, dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + C_shifted{i,j,k,l} = e_sh{m}'*obj.C{i,j,k,l}*e_sh{m}; + end + end + end + end + + % Loop along boundary + for i = 1:N + C_boundary_mat = [C_boundary{b,1,b,1}(i,i), C_boundary{b,1,b,2}(i,i); ... + C_boundary{b,2,b,1}(i,i), C_boundary{b,2,b,2}(i,i)]; + + C_shifted_mat = [C_shifted{b,1,b,1}(i,i), C_shifted{b,1,b,2}(i,i); ... + C_shifted{b,2,b,1}(i,i), C_shifted{b,2,b,2}(i,i)]; + + beta_i = obj.borrowingBisection(C_boundary_mat, C_shifted_mat); + beta(i) = min(beta(i), beta_i); + end + end + if min(beta) < 1e-10 + disp('WARNING! No borrowing from C seems to be possible.') + disp(['min(beta) = ' num2str(min(beta))]) + end + end + % Returns the component number that is the tangential/normal component % at the specified boundary function comp = getComponent(obj, comp_str, boundary)