sbplib: +scheme/Elastic2dStaggeredAnisotropic.m comparison

comparison +scheme/Elastic2dStaggeredAnisotropic.m @ 1279:546ee16760d5 feature/poroelastic

Add interfaces in Elastic2dStaggeredAnisotropic

author	Martin Almquist <malmquist@stanford.edu>
date	Sun, 07 Jun 2020 11:30:49 -0700
parents	e61083f178be
children

comparison

equal deleted inserted replaced

-:ce422ba8964e
+:546ee16760d5
 for a = 1:nGrids
 for i = 1:dim
 for j = 1:dim
 for k = 1:dim
 for l = 1:dim
-C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
+if numel(C) == nGrids
+C_mat{a}{i,j,k,l} = spdiag(C{a}{i,j,k,l});
+else
+C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
+end
 end
 end
 end
 end
 end
 obj.e_s_s = e_s_s;
 obj.e_n_s = e_n_s;
 % D1_u2s{a, b}{i} approximates ddi and
-%takes from u grid number b to s grid number a
+% takes from u grid number b to s grid number a
 % Some of D1_x2y{a, b} are 0.
 D1_u2s = cell(nGrids, nGrids);
 D1_s2u = cell(nGrids, nGrids);
 N_u = cell(nGrids, 1);
 m_tot = obj.N;
 % Preallocate
 [~, col] = size(tau{1}{1});
 closure = sparse(m_tot, m_tot);
-penalty = cell(nGrids, 1);
+penalty = cell(1, nGrids);
 for a = 1:nGrids
 [~, col] = size(e_u{a});
 penalty{a} = sparse(m_tot, col);
 end
 % Unknown boundary condition
 otherwise
 error('No such boundary condition: type = %s',type);
 end
+penalty = cell2mat(penalty);
 end
 % type     Struct that specifies the interface coupling.
 %          Fields:
 %          -- tuning:           penalty strength, defaults to 1.0
 u = obj;
 v = neighbour_scheme;
 % Operators, u side
-e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
+eu_u       = u.getBoundaryOperatorForScalarField('e_u', boundary);
+es_u       = u.getBoundaryOperatorForScalarField('e_s', boundary);
 tau_u     = u.getBoundaryOperator('tau', boundary);
-h11_u     = u.getBorrowing(boundary);
 nu_u      = u.getNormal(boundary);
-E_u = u.E;
+G_u = u.G;
+U_u = u.U;
 C_u = u.C;
-m_tot_u = u.grid.N();
+m_tot_u = u.N;
 % Operators, v side
-e_v       = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
+eu_v       = v.getBoundaryOperatorForScalarField('e_u', neighbour_boundary);
+es_v       = v.getBoundaryOperatorForScalarField('e_s', neighbour_boundary);
 tau_v     = v.getBoundaryOperator('tau', neighbour_boundary);
-h11_v     = v.getBorrowing(neighbour_boundary);
 nu_v      = v.getNormal(neighbour_boundary);
-E_v = v.E;
+G_v = v.G;
+U_v = v.U;
 C_v = v.C;
-m_tot_v = v.grid.N();
+m_tot_v = v.N;
-% Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
-flipFlag = false;
-e_v_flip = e_v;
-if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ...
-(strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ...
-(strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ...
-(strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ...
-(strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ...
-(strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ...
-(strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ...
-(strcmp(boundary,'e') && strcmp(neighbour_boundary,'e'))
-flipFlag = true;
-e_v_flip = fliplr(e_v);
-t1 = tau_v(:,1:2:end-1);
-t2 = tau_v(:,2:2:end);
-t1 = fliplr(t1);
-t2 = fliplr(t2);
-tau_v(:,1:2:end-1) = t1;
-tau_v(:,2:2:end) = t2;
-end
 % Operators that are only required for own domain
-Hi      = u.Hi_kron;
+% Hi      = u.Hi_kron;
-RHOi    = u.RHOi_kron;
+% RHOi    = u.RHOi_kron;
-e_kron  = u.getBoundaryOperator('e', boundary);
+% e_kron  = u.getBoundaryOperator('e', boundary);
+H       = u.H_u;
+RHO     = u.RHO;
 T_u     = u.getBoundaryTractionOperator(boundary);
 % Shared operators
 H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
-H_gamma_kron    = u.getBoundaryQuadrature(boundary);
+% H_gamma_kron    = u.getBoundaryQuadrature(boundary);
 dim             = u.dim;
+nGrids          = obj.nGrids;
 % Preallocate
-[~, m_int] = size(H_gamma);
+% [~, m_int] = size(H_gamma);
-closure = sparse(dim*m_tot_u, dim*m_tot_u);
+closure = sparse(m_tot_u, m_tot_u);
-penalty = sparse(dim*m_tot_u, dim*m_tot_v);
+penalty = sparse(m_tot_u, m_tot_v);
-% ---- Continuity of displacement ------
+%---- Grid layout -------
+% gu1 = xp o yp;
-% Y: symmetrizing part of penalty
+% gu2 = xd o yd;
-% Z: symmetric part of penalty
+% gs1 = xd o yp;
-% X = Y + Z.
+% gs2 = xp o yd;
+%------------------------
-% Loop over components to couple across interface
+switch boundary
+case {'w', 'e'}
+switch neighbour_boundary
+case {'w', 'e'}
+gridCombos = {{1,1,1}, {2,2,2}};
+case {'s', 'n'}
+gridCombos = {{1,1,2}, {2,2,1}};
+end
+case {'s', 'n'}
+switch neighbour_boundary
+case {'s', 'n'}
+gridCombos = {{2,1,1}, {1,2,2}};
+case {'w', 'e'}
+gridCombos = {{2,1,2}, {1,2,1}};
+end
+end
+% Symmetrizing part
+for c = 1:nGrids
+for m = 1:numel(gridCombos)
+gc = gridCombos{m};
+a = gc{1};
+b = gc{2};
+for i = 1:dim
+for j = 1:dim
+Y = 1/2*T_u{a,c}{j,i}';
+closure = closure + G_u{c}*U_u{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(eu_u{b}'*U_u{b}{j}'*G_u{b}') ));
+penalty = penalty - G_u{c}*U_u{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(eu_v{b}'*U_v{b}{j}'*G_v{b}') ));
+end
+end
+end
+end
+% Symmetric part
+for m = 1:numel(gridCombos)
+gc = gridCombos{m};
+a = gc{1};
+b = gc{2};
+bv = gc{3};
+h11_u = u.getBorrowing(b, boundary);
+h11_v = v.getBorrowing(bv, neighbour_boundary);
+for i = 1:dim
+for j = 1:dim
+Z_u = 0*es_u{b}'*es_u{b};
+Z_v = 0*es_v{bv}'*es_v{bv};
+for l = 1:dim
+for k = 1:dim
+Z_u = Z_u + es_u{b}'*nu_u(l)*C_u{b}{l,i,k,j}*nu_u(k)*es_u{b};
+Z_v = Z_v + es_v{bv}'*nu_v(l)*C_v{bv}{l,i,k,j}*nu_v(k)*es_v{bv};
+end
+end
+Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
+% X = es_u{b}'*Z*es_u{b};
+% X = Z;
+closure = closure + G_u{a}*U_u{a}{i}*((RHO{a}*H{a})\(eu_u{a}*Z'*H_gamma{b}*(eu_u{a}'*U_u{a}{j}'*G_u{a}' ) ));
+penalty = penalty - G_u{a}*U_u{a}{i}*((RHO{a}*H{a})\(eu_u{a}*Z'*H_gamma{b}*(eu_v{a}'*U_v{a}{j}'*G_v{a}' ) ));
+end
+end
+end
+% Continuity of traction
 for j = 1:dim
+for m = 1:numel(gridCombos)
-% Loop over components that penalties end up on
+gc = gridCombos{m};
-for i = 1:dim
+a = gc{1};
-Y = 1/2*T_u{j,i}';
+b = gc{2};
-Z_u = sparse(m_int, m_int);
+bv = gc{3};
-Z_v = sparse(m_int, m_int);
+closure = closure - 1/2*G_u{a}*U_u{a}{j}*(RHO{a}\(H{a}\(eu_u{a}*H_gamma{b}*tau_u{b}{j}')));
-for l = 1:dim
+penalty = penalty - 1/2*G_u{a}*U_u{a}{j}*(RHO{a}\(H{a}\(eu_u{a}*H_gamma{b}*tau_v{bv}{j}')));
-for k = 1:dim
+end
-Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
+end
-Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
-end
-end
-if flipFlag
-Z_v = rot90(Z_v,2);
-end
-Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
-X = Y + Z*e_u';
-closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
-penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}';
-end
-end
-% ---- Continuity of traction ------
-closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
-penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
-% ---- Multiply by inverse of density x quadraure ----
-closure = RHOi*Hi*closure;
-penalty = RHOi*Hi*penalty;
 end
 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
 error('Non-conforming interfaces not implemented yet.');
 H_b = obj.(['H_', boundary, '_s']);
 end
 function N = size(obj)
-N = obj.dim*prod(obj.m);
+N = length(obj.D);
 end
 end
 end

Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dStaggeredAnisotropic.m @ 1279:546ee16760d5 feature/poroelastic