sbplib: +scheme/Elastic2dVariableAnisotropicIncompatible.m comparison

comparison +scheme/Elastic2dVariableAnisotropicIncompatible.m @ 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

comparison

equal deleted inserted replaced

-:14f44e81e1e3
+:464e6f65c2c6
 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
 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
 methods
 switch order
 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);
 D1 = cell(dim,1);
 D2 = cell(dim,1);
 e_e = kron(e_scalar_e, I_dim);
 e_s = kron(e_scalar_s, I_dim);
 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_w = cell(nBP_R, 1);
-e_shifted_scalar_e = cell(nBP, 1);
+e_shifted_scalar_e = cell(nBP_R, 1);
-e_shifted_scalar_s = cell(nBP, 1);
+e_shifted_scalar_s = cell(nBP_R, 1);
-e_shifted_scalar_n = cell(nBP, 1);
+e_shifted_scalar_n = cell(nBP_R, 1);
-for i = 1:nBP-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;
 e_0_nBP{2}(i) = 1;
 e_m_nBP{1}(end+1-i) = 1;
 obj.m = m;
 obj.h = h;
 obj.order = order;
 obj.grid = g;
 obj.dim = dim;
-obj.nBP = nBP;
+obj.nBP_R = nBP_R;
 end
 % Closure functions return the operators applied to the own domain to close the boundary
 % bc = {comp, 'd', dComps},
 % where
 % 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};
 type = bc{2};
 if ischar(comp)
 comp = obj.getComponent(comp, boundary);
 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);
 H_gamma     = obj.getBoundaryQuadratureForScalarField(boundary);
 nu          = obj.getNormal(boundary);
 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
 end
 % Symmetric part only required on components with displacement BC.
 % (Otherwise it's not symmetric.)
+beta_vec = obj.computeBorrowingFromC(boundary);
-% 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_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;
 corner_weight = e*spdiag(corner_weight)*e';
 for l = 1:dim
 for k = 1:dim
 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;
 end
 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 error('Non-conforming interfaces not implemented yet.');
 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;
 beta = 0.5;
 err = (beta_max - beta_min)/2;
 end
 err = err/2;
 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)
 assertIsMember(comp_str, {'normal', 'tangential'});
 assertIsMember(boundary, {'w', 'e', 's', 'n'});

Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropicIncompatible.m @ 1346:464e6f65c2c6 feature/poroelastic