Mercurial > repos > public > sbplib
view +sbp/+implementations/intOpOP_orders_6to6_ratio2to1.m @ 1322:412b8ceafbc6 feature/poroelastic
Add Zt to output args for Elastic2dCurvilinearAnisotropic.interfaceNormalTangential
author | Martin Almquist <malmquist@stanford.edu> |
---|---|
date | Sat, 24 Oct 2020 20:58:26 -0700 |
parents | 4291731570bb |
children |
line wrap: on
line source
function [IC2F,IF2C,Hc,Hf] = IntOpOP_orders_6to6_ratio2to1(mc,hc,ACC) % ACC is a string. % ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C. % ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F. ratio = 2; mf = ratio*mc-1; hf = hc/ratio; switch ACC case 'F2C' [stencil_F2C,BC_F2C,HcU,HfU] = ... sbp.implementations.intOpOP_orders_6to6_ratio_2to1_accC2F3_accF2C4; case 'C2F' [stencil_F2C,BC_F2C,HcU,HfU] = ... sbp.implementations.intOpOP_orders_6to6_ratio_2to1_accC2F4_accF2C3; end stencil_width = length(stencil_F2C); stencil_hw = (stencil_width-1)/2; [BC_rows,BC_cols] = size(BC_F2C); %%% Norm matrices %%% Hc = speye(mc,mc); HcUm = length(HcU); Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm); Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm); Hc = Hc*hc; Hf = speye(mf,mf); HfUm = length(HfU); Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm); Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm); Hf = Hf*hf; %%%%%%%%%%%%%%%%%%%%%% %%% Create IF2C from stencil and BC IF2C = sparse(mc,mf); for i = BC_rows+1 : mc-BC_rows IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX> end IF2C(1:BC_rows,1:BC_cols) = BC_F2C; IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Create IC2F using symmetry condition %%%% IC2F = Hf\IF2C.'*Hc;