Description of system_mat_1st

0001 function s_mat= system_mat_1st_order( fwd_model, img)
0002 % SYSTEM_MAT_1ST_ORDER: SS= system_mat_1st_order( fwd_model, img)
0003 % Calc system matrix for Andy Adler's EIT code
0004 % img.fwd_model = forward model
0005 % img       = image background for system matrix calc
0006 % s_mat.E = CC' * SS * conductivites * CC;
0007 % where:
0008 %   SS  = Unconnected system Matrix
0009 %   CC  = Connectivity Matrix
0010 
0011 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: system_mat_1st_order.m 5980 2019-06-23 07:34:48Z aadler $
0013 
0014 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0015 
0016 if nargin == 1
0017    img= fwd_model;
0018 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0019    warning('EIDORS:DeprecatedInterface', ...
0020       ['Calling SYSTEM_MAT_1ST_ORDER with two arguments is deprecated and will cause' ...
0021        ' an error in a future version. First argument ignored.']);
0022 end
0023 fwd_model= img.fwd_model;
0024 
0025 % check parametrizations
0026 if isfield(img,'current_params') && ...
0027      (~ischar(img.current_params) || ~strcmp(img.current_params,'conductivity'))
0028     if ischar(img.current_params)
0029        error('system_mat_1st_order does not work for %s',img.current_params);
0030     else
0031        img.current_params
0032        error('system_mat_1st_order does not work for the above');
0033     end
0034 end
0035 
0036 FC= system_mat_fields( fwd_model);
0037 lFC= size(FC,1);
0038 ED = elem_dim(fwd_model);
0039 lNE= ED*num_elems(fwd_model);
0040 
0041 elem_data = check_elem_data(fwd_model, img);
0042 if size(elem_data,3) == 1
0043 % Scalar conductivity == isotropic
0044    elem_sigma = kron( elem_data, ones(ED,1) );
0045    elem_sigma(end+1:lFC) = 1; % add ones for CEM
0046 
0047    ES= spdiags(elem_sigma,0,lFC,lFC);
0048 else
0049    switch elem_dim(fwd_model)
0050      case 2;
0051        idx = 1:2:lNE;
0052        ES= sparse([idx,idx+1,idx,idx+1]', ...
0053                   [idx,idx,idx+1,idx+1]', elem_data(:), lFC,lFC);
0054        
0055        ES(lNE+1:lFC,lNE+1:lFC) = speye(lFC-lNE);
0056      case 3;
0057        idx = 1:3:lNE;
0058        ES= sparse([idx,idx+1,idx+2,idx,idx+1,idx+2,idx,idx+1,idx+2]', ...
0059                   [idx,idx,idx,idx+1,idx+1,idx+1,idx+2,idx+2,idx+2]', ...
0060                    elem_data(:), lFC,lFC);
0061        
0062        ES(lNE+1:lFC,lNE+1:lFC) = speye(lFC-lNE);
0063      otherwise; 
0064        error('%d D anisotropic elements not implemented', elem_dim(fwd_model));
0065    end
0066 
0067 end
0068 
0069 % System matrix is symmetric (but not hermitian)
0070 %   so use the .'  (not ')
0071 E= FC.' * ES * FC;
0072 % This will force E to by symmetric, so that the
0073 % solver won't mistakenly use LU rather than Chol
0074 s_mat.E= 1/2*(E.' + E);
0075 
0076 function elem_data = check_elem_data(fwd_model, img);
0077    elem_data = img.elem_data; 
0078    sz_elem_data = size(elem_data);
0079    if sz_elem_data(2) ~= 1;
0080       error('system_mat_1st_order: can only solve one image (sz_elem_data=%)', ...
0081             sz_elem_data);
0082    end
0083 
0084    if isfield(fwd_model, 'coarse2fine');
0085      c2f = fwd_model.coarse2fine;
0086      sz_c2f = size(c2f);
0087      switch sz_elem_data(1)
0088        case sz_c2f(1); % Ok
0089        case sz_c2f(2); elem_data = c2f * elem_data;
0090           if isfield(fwd_model, 'background')
0091               elem_data = elem_data + fwd_model.background;
0092           end
0093 
0094        otherwise;
0095            error(['system_mat_1st_order: provided elem_data ' ...
0096                   ' (sz=%d) does not match c2f (sz=[%d %d])'], ...
0097                    sz_elem_data(1), sz_c2f(1), sz_c2f(2));
0098      end
0099    else
0100      if sz_elem_data(1) ~= num_elems(fwd_model)
0101        error(['system_mat_1st_order: provided elem_data (sz=%d) does ' ...
0102           ' not match fwd_model (sz=%d)'], sz_elem_data(1), num_elems(fwd_model));
0103      end
0104    end
0105 
0106 
0107 % Need tests for each error case
0108 % Need tests for background
0109 
0110 function do_unit_test
0111    imdl=  mk_common_model('a2c2',16);
0112    img=  mk_image(imdl);
0113    S1 = system_mat_1st_order(img.fwd_model,img);
0114    unit_test_cmp('sys_mat1', S1.E(1,:), [4,-1,-1,-1,-1,zeros(1,36)],1e-14);
0115 
0116    img.elem_data([1:16]) = 2;
0117    img.calc_colours.clim = 4; show_fem(img,[0,0,3]);
0118 
0119    S2 = system_mat_1st_order(img.fwd_model,img);
0120    unit_test_cmp('sys_mat2', S2.E(1,:), 2*[4,-1,-1,-1,-1,zeros(1,36)],1e-14);
0121 
0122    imdl=  mk_common_model('a2C2',16);
0123    img=  mk_image(imdl);
0124    S1 = system_mat_1st_order(img.fwd_model,img);
0125    unit_test_cmp('sys_mat1', S1.E(1,:), [4,-1,-1,-1,-1,zeros(1,52)],1e-14);
0126 
0127    img.elem_data([1:16]) = 2;
0128    img.calc_colours.clim = 4; show_fem(img,[0,0,3]);
0129 
0130    S2 = system_mat_1st_order(img.fwd_model,img);
0131    unit_test_cmp('sys_mat2', S2.E(1,:), 2*[4,-1,-1,-1,-1,zeros(1,52)],1e-14);
0132    
0133    idx = 41-(0:15);
0134    unit_test_cmp('sys_mat4', S2.E(idx,idx),   S1.E(idx,idx),1e-14);
0135    idx = 1:5;
0136    unit_test_cmp('sys_mat3', S2.E(idx,idx), 2*S1.E(idx,idx),1e-14);
0137 
0138    img.elem_data([1:36]) = 2;
0139    S2 = system_mat_1st_order(img.fwd_model,img);
0140    idx = 1:5;
0141    unit_test_cmp('sys_mat3', S2.E(idx,idx), 2*S1.E(idx,idx),1e-14);
0142 
0143    img.elem_data = ones(63,1);
0144    try
0145      % this should give an error, since elem_data is wrong size
0146       S2 = system_mat_1st_order(img.fwd_model,img);
0147       unit_test_cmp('sys_mat: test for size error', 1,0);
0148    catch
0149       unit_test_cmp('sys_mat: test for size error', 1,1);
0150    end
0151 
0152    imdl=  mk_common_model('a2c2',16);
0153    cmdl = mk_circ_tank(4,[],0);
0154    c2f = mk_coarse_fine_mapping( img.fwd_model, cmdl); % spy(c2f)
0155    img.fwd_model.coarse2fine = c2f;
0156    img.elem_data = ones(size(cmdl.elems,1),1);
0157    S3 = system_mat_1st_order(img.fwd_model,img);
0158    unit_test_cmp('c2f #1', S1.E, S3.E,1e-14);
0159 
0160    cmdl = mk_circ_tank(2,[],0);
0161    c2f = mk_coarse_fine_mapping( img.fwd_model, cmdl); % spy(c2f)
0162    c2f = c2f./(sum(c2f,2) * ones(1,size(c2f,2))); % FIX
0163    img.fwd_model.coarse2fine = c2f;
0164    img.elem_data = ones(size(cmdl.elems,1),1);
0165    S4 = system_mat_1st_order(img.fwd_model,img);
0166    unit_test_cmp('c2f #2', S1.E, S4.E,1e-14);
0167 
0168    test_2d_resistor
0169    test_3d_resistor
0170 
0171 function test_2d_resistor
0172    nn= 12;     % number of nodes
0173    ww=2;       % width = 4
0174    conduc=  .4;% conductivity in Ohm-meters
0175    current= 4;  % Amps
0176    z_contact= 9e-1;
0177    scale = .35;
0178    mdl=mk_grid_model([],3+scale*(1:ww), scale*(1:nn/ww));
0179    mdl= rmfield(mdl,'coarse2fine'); % don't calc this.
0180 
0181    mdl.gnd_node = 1;
0182    elec_nodes= [1:ww];
0183    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0184    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0185    stim.stim_pattern= [-1;1]*current;
0186    stim.meas_pattern= [-1,1];
0187    mdl.stimulation= stim;
0188    mdl.electrode= elec;
0189    mdl.normalize_measurements = 0;
0190    %show_fem(mdl);
0191 
0192    % analytical solution
0193    wid_len= max(mdl.nodes) - min(mdl.nodes);
0194    R = wid_len(2) / wid_len(1) / conduc + 2*z_contact/scale;
0195 
0196    V= current*R;
0197    %fprintf('Solver %s: %f\n', 'analytic', V);
0198 
0199    % AA_SOLVER
0200    mdl.solve = @fwd_solve_1st_order;
0201    mdl.system_mat = @system_mat_1st_order;
0202    img = mk_image(mdl, conduc);
0203    fsol= fwd_solve(img);
0204    %fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0205    unit_test_cmp('test 2d vs analytic', fsol.meas, V, 1e-11);
0206 
0207    % NP_SOLVER
0208    mdl.solve = @np_fwd_solve;
0209    mdl.system_mat = @np_calc_system_mat;
0210    img = mk_image(mdl, conduc);
0211    fsol= fwd_solve(img);
0212    %fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0213    unit_test_cmp('np_fwd_solve vs analytic', fsol.meas, V, 1e-11);
0214 
0215 function test_3d_resistor
0216    ll=5; % length
0217    ww=1; % width
0218    hh=1; % height
0219    conduc= .13;  % conductivity in Ohm-meters
0220    current= 4;  % Amps
0221    z_contact= 9e-1;
0222    scale = .46;
0223    nn=0;
0224    for z=0:ll; for x=0:ww; for y=0:hh
0225       nn=nn+1;
0226       mdl.nodes(nn,:) = [x,y,z];
0227    end; end; end
0228 
0229    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0230    mdl.nodes= mdl.nodes*scale;
0231    mdl= rmfield(mdl,'coarse2fine');
0232    mdl.boundary= find_boundary(mdl.elems);
0233    mdl.gnd_node = 1;
0234    mdl.normalize_measurements = 0;
0235    elec_nodes= [1:(ww+1)*(hh+1)];
0236    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0237    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0238    stim.stim_pattern= [-1;1]*current;
0239    stim.meas_pattern= [-1,1];
0240    mdl.stimulation= stim;
0241    mdl.electrode= elec;
0242 %  show_fem(mdl);
0243 
0244    % analytical solution
0245    R = ll / ww / hh / scale/ conduc + 2*z_contact/scale^2;
0246 
0247    V= current*R;
0248 %  fprintf('Solver %s: %f\n', 'analytic', V);
0249 
0250    mdl.solve = @fwd_solve_1st_order;
0251    mdl.system_mat = @system_mat_1st_order;
0252 
0253    img= mk_image(mdl, conduc);
0254 
0255    fsol= fwd_solve(img);
0256 %  fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0257    unit_test_cmp('test 2d vs analytic', fsol.meas, V, 1e-11);
0258 
0259    mdl.solve = @np_fwd_solve;
0260    mdl.system_mat = @np_calc_system_mat;
0261    mdl.misc.perm_sym= '{n}';
0262    img= mk_image(mdl, conduc);
0263 
0264    fsol= fwd_solve(img);
0265 %  fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0266    unit_test_cmp('np_fwd_solve vs analytic', fsol.meas, V, 1e-11);
system_mat_1st_order

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
