view +sbp/+implementations/d2_blocknorm_10.m @ 1347:ac54767ae1fb feature/poroelastic tip

Add interface, not fully compatible.
author Martin Almquist <martin.almquist@it.uu.se>
date Tue, 30 Apr 2024 14:58:35 +0200
parents f7ac3cd6eeaa
children
line wrap: on
line source

function [H, HI, D1, D2, e_1, e_m, M, Q, S_1, S_m] = d2_blocknorm_10(m,h)

    BP = 9;
    if(m<2*BP)
        error(['Operator requires at least ' num2str(2*BP) ' grid points']);
    end

    H_U=[0.428081020217e12 / 0.2633637888000e13 0.779032713983e12 / 0.2633637888000e13 -0.1187642619571e13 / 0.2633637888000e13 0.1642279196603e13 / 0.2633637888000e13 -0.339289243121e12 / 0.526727577600e12 0.1261055176253e13 / 0.2633637888000e13 -0.658216413073e12 / 0.2633637888000e13 0.33968779823e11 / 0.376233984000e12 -0.764998e6 / 0.40186125e8; 0.779032713983e12 / 0.2633637888000e13 0.317907052061e12 / 0.164602368000e12 -0.1082918052397e13 / 0.658409472000e12 0.176473501369e12 / 0.82301184000e11 -0.521625191587e12 / 0.263363788800e12 0.200523313337e12 / 0.164602368000e12 -0.279496000009e12 / 0.658409472000e12 -0.832e3 / 0.40186125e8 0.129794760887e12 / 0.2633637888000e13; -0.1187642619571e13 / 0.2633637888000e13 -0.1082918052397e13 / 0.658409472000e12 0.1304108863849e13 / 0.329204736000e12 -0.533093695961e12 / 0.131681894400e12 0.213949854133e12 / 0.52672757760e11 -0.1883976009151e13 / 0.658409472000e12 0.51084128e8 / 0.40186125e8 -0.92763684343e11 / 0.658409472000e12 -0.59058717923e11 / 0.526727577600e12; 0.1642279196603e13 / 0.2633637888000e13 0.176473501369e12 / 0.82301184000e11 -0.533093695961e12 / 0.131681894400e12 0.217677310051e12 / 0.32920473600e11 -0.1509120465127e13 / 0.263363788800e12 0.170839232e9 / 0.40186125e8 -0.1404096707137e13 / 0.658409472000e12 0.11789520859e11 / 0.32920473600e11 0.85652315431e11 / 0.526727577600e12; -0.339289243121e12 / 0.526727577600e12 -0.521625191587e12 / 0.263363788800e12 0.213949854133e12 / 0.52672757760e11 -0.1509120465127e13 / 0.263363788800e12 0.21849109e8 / 0.3214890e7 -0.1134422468377e13 / 0.263363788800e12 0.602448430967e12 / 0.263363788800e12 -0.4910542309e10 / 0.10534551552e11 -0.83039945231e11 / 0.526727577600e12; 0.1261055176253e13 / 0.2633637888000e13 0.200523313337e12 / 0.164602368000e12 -0.1883976009151e13 / 0.658409472000e12 0.170839232e9 / 0.40186125e8 -0.1134422468377e13 / 0.263363788800e12 0.681437038097e12 / 0.164602368000e12 -0.1108257453763e13 / 0.658409472000e12 0.31631872327e11 / 0.82301184000e11 0.37820115539e11 / 0.376233984000e12; -0.658216413073e12 / 0.2633637888000e13 -0.279496000009e12 / 0.658409472000e12 0.51084128e8 / 0.40186125e8 -0.1404096707137e13 / 0.658409472000e12 0.602448430967e12 / 0.263363788800e12 -0.1108257453763e13 / 0.658409472000e12 0.623491124887e12 / 0.329204736000e12 -0.146643738067e12 / 0.658409472000e12 -0.98874149197e11 / 0.2633637888000e13; 0.33968779823e11 / 0.376233984000e12 -0.832e3 / 0.40186125e8 -0.92763684343e11 / 0.658409472000e12 0.11789520859e11 / 0.32920473600e11 -0.4910542309e10 / 0.10534551552e11 0.31631872327e11 / 0.82301184000e11 -0.146643738067e12 / 0.658409472000e12 0.174599973347e12 / 0.164602368000e12 0.4625165773e10 / 0.526727577600e12; -0.764998e6 / 0.40186125e8 0.129794760887e12 / 0.2633637888000e13 -0.59058717923e11 / 0.526727577600e12 0.85652315431e11 / 0.526727577600e12 -0.83039945231e11 / 0.526727577600e12 0.37820115539e11 / 0.376233984000e12 -0.98874149197e11 / 0.2633637888000e13 0.4625165773e10 / 0.526727577600e12 0.525286231387e12 / 0.526727577600e12;];


    H=speye(m);
    H(1:9,1:9)=H_U;
    H(m-8:m,m-8:m)=rot90( H_U(1:9,1:9) ,2 );
    H=H*h;
    HI=inv(H);

%     Q=1/1260*diag(ones(m-5,1),5)-5/504*diag(ones(m-4,1),4)+5/84*diag(ones(m-3,1),3)-5/21*diag(ones(m-2,1),2)+5/6*diag(ones(m-1,1),1)-5/6*diag(ones(m-1,1),-1)+5/21*diag(ones(m-2,1),-2)-5/84*diag(ones(m-3,1),-3)+5/504*diag(ones(m-4,1),-4)-1/1260*diag(ones(m-5,1),-5) ;
    diags   = -5:5;
    stencil = [-1/1260,5/504,-5/84,5/21,-5/6,0,5/6,-5/21,5/84,-5/504,1/1260];
    Q = stripeMatrix(stencil, diags, m);

    Q_U = [-0.1e1 / 0.2e1 0.78249683e8 / 0.71442000e8 -0.28290472447e11 / 0.18289152000e11 0.4285528063e10 / 0.2032128000e10 -0.3924872557e10 / 0.1828915200e10 0.2856344621e10 / 0.1828915200e10 -0.1598284927e10 / 0.2032128000e10 0.4652402687e10 / 0.18289152000e11 -0.731623e6 / 0.17860500e8; -0.78249683e8 / 0.71442000e8 0 0.5351945471e10 / 0.2032128000e10 -0.4078428731e10 / 0.1143072000e10 0.4879509877e10 / 0.1219276800e10 -0.1273709579e10 / 0.381024000e9 0.7433128649e10 / 0.3657830400e10 -0.1239008e7 / 0.1488375e7 0.21251483e8 / 0.124416000e9; 0.28290472447e11 / 0.18289152000e11 -0.5351945471e10 / 0.2032128000e10 0 0.6701592799e10 / 0.2612736000e10 -0.57535927e8 / 0.17418240e8 0.3058732543e10 / 0.870912000e9 -0.28821953e8 / 0.10206000e8 0.1858901437e10 / 0.1219276800e10 -0.2363118211e10 / 0.6096384000e10; -0.4285528063e10 / 0.2032128000e10 0.4078428731e10 / 0.1143072000e10 -0.6701592799e10 / 0.2612736000e10 0 0.52016695e8 / 0.20901888e8 -0.3777923e7 / 0.1275750e7 0.123333949e9 / 0.41472000e8 -0.2256321727e10 / 0.1143072000e10 0.1058261459e10 / 0.1828915200e10; 0.3924872557e10 / 0.1828915200e10 -0.4879509877e10 / 0.1219276800e10 0.57535927e8 / 0.17418240e8 -0.52016695e8 / 0.20901888e8 0 0.125390297e9 / 0.58060800e8 -0.609569351e9 / 0.261273600e9 0.438488399e9 / 0.243855360e9 -0.14244569e8 / 0.24385536e8; -0.2856344621e10 / 0.1828915200e10 0.1273709579e10 / 0.381024000e9 -0.3058732543e10 / 0.870912000e9 0.3777923e7 / 0.1275750e7 -0.125390297e9 / 0.58060800e8 0 0.4624381729e10 / 0.2612736000e10 -0.484096919e9 / 0.381024000e9 0.2676438019e10 / 0.6096384000e10; 0.1598284927e10 / 0.2032128000e10 -0.7433128649e10 / 0.3657830400e10 0.28821953e8 / 0.10206000e8 -0.123333949e9 / 0.41472000e8 0.609569351e9 / 0.261273600e9 -0.4624381729e10 / 0.2612736000e10 0 0.21500967689e11 / 0.18289152000e11 -0.7199454721e10 / 0.18289152000e11; -0.4652402687e10 / 0.18289152000e11 0.1239008e7 / 0.1488375e7 -0.1858901437e10 / 0.1219276800e10 0.2256321727e10 / 0.1143072000e10 -0.438488399e9 / 0.243855360e9 0.484096919e9 / 0.381024000e9 -0.21500967689e11 / 0.18289152000e11 0 0.761653e6 / 0.882000e6; 0.731623e6 / 0.17860500e8 -0.21251483e8 / 0.124416000e9 0.2363118211e10 / 0.6096384000e10 -0.1058261459e10 / 0.1828915200e10 0.14244569e8 / 0.24385536e8 -0.2676438019e10 / 0.6096384000e10 0.7199454721e10 / 0.18289152000e11 -0.761653e6 / 0.882000e6 0;];

    Q(1:9,1:9)=Q_U;
    Q(m-8:m,m-8:m)=rot90( -Q_U(1:9,1:9) ,2 );

    D1=H\Q;

    M_U =[0.3812926003e10 / 0.2438553600e10 -0.5433856529e10 / 0.1741824000e10 0.4187462879e10 / 0.1045094400e10 -0.65635105447e11 / 0.12192768000e11 0.1457682577e10 / 0.270950400e9 -0.27884016067e11 / 0.7315660800e10 0.22304839493e11 / 0.12192768000e11 -0.188132543e9 / 0.348364800e9 0.42711619e8 / 0.571536000e9; -0.5433885329e10 / 0.1741824000e10 0.23985229969e11 / 0.2286144000e10 -0.10208460799e11 / 0.571536000e9 0.8828370001e10 / 0.381024000e9 -0.12306735263e11 / 0.522547200e9 0.39313626089e11 / 0.2286144000e10 -0.4380200287e10 / 0.508032000e9 0.192498023e9 / 0.71442000e8 -0.14565232681e11 / 0.36578304000e11; 0.29313778073e11 / 0.7315660800e10 -0.10209211399e11 / 0.571536000e9 0.24157533391e11 / 0.653184000e9 -0.6561725111e10 / 0.130636800e9 0.26490755639e11 / 0.522547200e9 -0.352801289e9 / 0.9331200e7 0.401374423e9 / 0.20412000e8 -0.29541854057e11 / 0.4572288000e10 0.7396364989e10 / 0.7315660800e10; -0.65653163047e11 / 0.12192768000e11 0.8832999601e10 / 0.381024000e9 -0.6566554871e10 / 0.130636800e9 0.1575758731e10 / 0.21772800e8 -0.2571648133e10 / 0.34836480e8 0.556969019e9 / 0.10206000e8 -0.3139265911e10 / 0.108864000e9 0.165424529e9 / 0.16934400e8 -0.11623567549e11 / 0.7315660800e10; 0.1458632977e10 / 0.270950400e9 -0.86260417241e11 / 0.3657830400e10 0.3792749777e10 / 0.74649600e8 -0.2576699653e10 / 0.34836480e8 0.157840723e9 / 0.2041200e7 -0.29838889141e11 / 0.522547200e9 0.5142742211e10 / 0.174182400e9 -0.36792522023e11 / 0.3657830400e10 0.2438136689e10 / 0.1463132160e10; -0.27909676867e11 / 0.7315660800e10 0.39389053289e11 / 0.2286144000e10 -0.2478149663e10 / 0.65318400e8 0.559214069e9 / 0.10206000e8 -0.29914661941e11 / 0.522547200e9 0.2843819551e10 / 0.65318400e8 -0.14816015149e11 / 0.653184000e9 0.1677001673e10 / 0.228614400e9 -0.44441740171e11 / 0.36578304000e11; 0.3188985299e10 / 0.1741824000e10 -0.626864041e9 / 0.72576000e8 0.57539389e8 / 0.2916000e7 -0.3153500311e10 / 0.108864000e9 0.5162239811e10 / 0.174182400e9 -0.14840163949e11 / 0.653184000e9 0.419025709e9 / 0.31104000e8 -0.138945749e9 / 0.27216000e8 0.591880819e9 / 0.746496000e9; -0.1317440441e10 / 0.2438553600e10 0.192696473e9 / 0.71442000e8 -0.4230355151e10 / 0.653184000e9 0.165983249e9 / 0.16934400e8 -0.36905792423e11 / 0.3657830400e10 0.1679779433e10 / 0.228614400e9 -0.972870443e9 / 0.190512000e9 0.9129544111e10 / 0.2286144000e10 -0.13387742111e11 / 0.7315660800e10; 0.42721069e8 / 0.571536000e9 -0.14572922281e11 / 0.36578304000e11 0.7407199549e10 / 0.7315660800e10 -0.11649228349e11 / 0.7315660800e10 0.349038407e9 / 0.209018880e9 -0.44495912971e11 / 0.36578304000e11 0.29009849731e11 / 0.36578304000e11 -0.13387863071e11 / 0.7315660800e10 0.21585797479e11 / 0.7315660800e10;];


%     T=[-0.1e1 / 0.3150e4 0.5e1 / 0.1008e4 -0.5e1 / 0.126e3 0.5e1 / 0.21e2 -0.5e1 / 0.3e1 0.5269e4 / 0.1800e4 -0.5e1 / 0.3e1 0.5e1 / 0.21e2 -0.5e1 / 0.126e3 0.5e1 / 0.1008e4 -0.1e1 / 0.3150e4;];
%     M=(T(1)*diag(ones(m-5,1),5)+T(2)*diag(ones(m-4,1),4)+T(3)*diag(ones(m-3,1),3)+T(4)*diag(ones(m-2,1),2)+T(5)*diag(ones(m-1,1),1)+T(7)*diag(ones(m-1,1),-1)+T(8)*diag(ones(m-2,1),-2)+T(9)*diag(ones(m-3,1),-3)+T(10)*diag(ones(m-4,1),-4)+T(11)*diag(ones(m-5,1),-5)+T(6)*diag(ones(m,1),0));

    diags   = -5:5;
    left_stencil = [-1/3150,5/1008,-5/126,5/21,-5/3];
    stencil = [left_stencil,5269/1800,fliplr(left_stencil)];
    M = stripeMatrix(stencil, diags, m);

    M(1:9,1:9)=M_U;

    M(m-8:m,m-8:m)=rot90(  M_U ,2 );
    M=M/h;

    DS_U=[0.761e3 / 0.280e3 -8 14 -0.56e2 / 0.3e1 0.35e2 / 0.2e1 -0.56e2 / 0.5e1 0.14e2 / 0.3e1 -0.8e1 / 0.7e1 0.1e1 / 0.8e1;];
    DS=sparse(m,m);
    DS(1,1:9)=DS_U;
    DS(m,m-8:m)=fliplr(DS_U);
    DS=DS/h;

    D2=H\(-M+DS);

    % Try adding AD to boundary D1=HI*(Q-D9'*D9)
%     DD_9=sparse(m);
%     d9=[-1 9 -36 84 -126 126 -84 36 -9 1];t9=sum(abs(d9));%d9=d9/t9;
%     DD_9(1:1,1:10)=[d9];
%     DD_9(m:m,m-9:m)=[d9];
% 
% 
%     ADD=30*h/(t9)*DD_9'*DD_9;

    e_1 = sparse(m,1);
    e_1(1)= 1;
    e_m = sparse(m,1);
    e_m(end)= 1;
    S_1 = -DS(1,:)';
    S_m =  DS(end,:)';

    Q = H*D1-(-(e_1*e_1') + (e_m*e_m'));
    M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
end