system_mat_1st_order

PURPOSE ^

SYSTEM_MAT_1ST_ORDER: SS= system_mat_1st_order( fwd_model, img)

SYNOPSIS ^

function s_mat= system_mat_1st_order( fwd_model, img)

DESCRIPTION ^

 SYSTEM_MAT_1ST_ORDER: SS= system_mat_1st_order( fwd_model, img)
 Calc system matrix for Andy Adler's EIT code
 img.fwd_model = forward model
 img       = image background for system matrix calc
 s_mat.E = CC' * SS * conductivites * CC;
 where:
   SS  = Unconnected system Matrix
   CC  = Connectivity Matrix

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5574 2017-06-21 02:24:41Z 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; error(['system_mat_1st_order: provided elem_data ' ...
0095             ' (sz=%d) does not match c2f (sz=%d %d)'], sz_elem_data(1), sz_c2f);
0096      end
0097    else
0098      if sz_elem_data(1) ~= num_elems(fwd_model)
0099        error(['system_mat_1st_order: provided elem_data (sz=%d) does ' ...
0100           ' not match fwd_model (sz=%d)'], sz_elem_data(1), num_elems(fwd_model));
0101      end
0102    end
0103 
0104 
0105 % Need tests for each error case
0106 % Need tests for background
0107 
0108 function do_unit_test
0109    imdl=  mk_common_model('a2c2',16);
0110    img=  mk_image(imdl);
0111    S1 = system_mat_1st_order(img.fwd_model,img);
0112    unit_test_cmp('sys_mat1', S1.E(1,:), [4,-1,-1,-1,-1,zeros(1,36)],1e-14);
0113 
0114    img.elem_data([1:16]) = 2;
0115    img.calc_colours.clim = 4; show_fem(img,[0,0,3]);
0116 
0117    S2 = system_mat_1st_order(img.fwd_model,img);
0118    unit_test_cmp('sys_mat2', S2.E(1,:), 2*[4,-1,-1,-1,-1,zeros(1,36)],1e-14);
0119 
0120    imdl=  mk_common_model('a2C2',16);
0121    img=  mk_image(imdl);
0122    S1 = system_mat_1st_order(img.fwd_model,img);
0123    unit_test_cmp('sys_mat1', S1.E(1,:), [4,-1,-1,-1,-1,zeros(1,52)],1e-14);
0124 
0125    img.elem_data([1:16]) = 2;
0126    img.calc_colours.clim = 4; show_fem(img,[0,0,3]);
0127 
0128    S2 = system_mat_1st_order(img.fwd_model,img);
0129    unit_test_cmp('sys_mat2', S2.E(1,:), 2*[4,-1,-1,-1,-1,zeros(1,52)],1e-14);
0130    
0131    idx = 41-(0:15);
0132    unit_test_cmp('sys_mat4', S2.E(idx,idx),   S1.E(idx,idx),1e-14);
0133    idx = 1:5;
0134    unit_test_cmp('sys_mat3', S2.E(idx,idx), 2*S1.E(idx,idx),1e-14);
0135 
0136    img.elem_data([1:36]) = 2;
0137    S2 = system_mat_1st_order(img.fwd_model,img);
0138    idx = 1:5;
0139    unit_test_cmp('sys_mat3', S2.E(idx,idx), 2*S1.E(idx,idx),1e-14);
0140 
0141    img.elem_data = ones(63,1);
0142    try
0143      % this should give an error, since elem_data is wrong size
0144       S2 = system_mat_1st_order(img.fwd_model,img);
0145       unit_test_cmp('sys_mat: test for size error', 1,0);
0146    catch
0147       unit_test_cmp('sys_mat: test for size error', 1,1);
0148    end
0149 
0150    imdl=  mk_common_model('a2c2',16);
0151    cmdl = mk_circ_tank(4,[],0);
0152    c2f = mk_coarse_fine_mapping( img.fwd_model, cmdl); % spy(c2f)
0153    img.fwd_model.coarse2fine = c2f;
0154    img.elem_data = ones(size(cmdl.elems,1),1);
0155    S3 = system_mat_1st_order(img.fwd_model,img);
0156    unit_test_cmp('c2f #1', S1.E, S3.E,1e-14);
0157 
0158    cmdl = mk_circ_tank(2,[],0);
0159    c2f = mk_coarse_fine_mapping( img.fwd_model, cmdl); % spy(c2f)
0160    c2f = c2f./(sum(c2f,2) * ones(1,size(c2f,2))); % FIX
0161    img.fwd_model.coarse2fine = c2f;
0162    img.elem_data = ones(size(cmdl.elems,1),1);
0163    S4 = system_mat_1st_order(img.fwd_model,img);
0164    unit_test_cmp('c2f #2', S1.E, S4.E,1e-14);
0165 
0166    test_2d_resistor
0167    test_3d_resistor
0168 
0169 function test_2d_resistor
0170    nn= 12;     % number of nodes
0171    ww=2;       % width = 4
0172    conduc=  .4;% conductivity in Ohm-meters
0173    current= 4;  % Amps
0174    z_contact= 9e-1;
0175    scale = .35;
0176    mdl=mk_grid_model([],3+scale*(1:ww), scale*(1:nn/ww));
0177    mdl= rmfield(mdl,'coarse2fine'); % don't calc this.
0178 
0179    mdl.gnd_node = 1;
0180    elec_nodes= [1:ww];
0181    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0182    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0183    stim.stim_pattern= [-1;1]*current;
0184    stim.meas_pattern= [-1,1];
0185    mdl.stimulation= stim;
0186    mdl.electrode= elec;
0187    mdl.normalize_measurements = 0;
0188    %show_fem(mdl);
0189 
0190    % analytical solution
0191    wid_len= max(mdl.nodes) - min(mdl.nodes);
0192    R = wid_len(2) / wid_len(1) / conduc + 2*z_contact/scale;
0193 
0194    V= current*R;
0195    %fprintf('Solver %s: %f\n', 'analytic', V);
0196 
0197    % AA_SOLVER
0198    mdl.solve = @fwd_solve_1st_order;
0199    mdl.system_mat = @system_mat_1st_order;
0200    img = mk_image(mdl, conduc);
0201    fsol= fwd_solve(img);
0202    %fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0203    unit_test_cmp('test 2d vs analytic', fsol.meas, V, 1e-11);
0204 
0205    % NP_SOLVER
0206    mdl.solve = @np_fwd_solve;
0207    mdl.system_mat = @np_calc_system_mat;
0208    img = mk_image(mdl, conduc);
0209    fsol= fwd_solve(img);
0210    %fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0211    unit_test_cmp('np_fwd_solve vs analytic', fsol.meas, V, 1e-11);
0212 
0213 function test_3d_resistor
0214    ll=5; % length
0215    ww=1; % width
0216    hh=1; % height
0217    conduc= .13;  % conductivity in Ohm-meters
0218    current= 4;  % Amps
0219    z_contact= 9e-1;
0220    scale = .46;
0221    nn=0;
0222    for z=0:ll; for x=0:ww; for y=0:hh
0223       nn=nn+1;
0224       mdl.nodes(nn,:) = [x,y,z];
0225    end; end; end
0226 
0227    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0228    mdl.nodes= mdl.nodes*scale;
0229    mdl= rmfield(mdl,'coarse2fine');
0230    mdl.boundary= find_boundary(mdl.elems);
0231    mdl.gnd_node = 1;
0232    mdl.normalize_measurements = 0;
0233    elec_nodes= [1:(ww+1)*(hh+1)];
0234    elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0235    elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0236    stim.stim_pattern= [-1;1]*current;
0237    stim.meas_pattern= [-1,1];
0238    mdl.stimulation= stim;
0239    mdl.electrode= elec;
0240 %  show_fem(mdl);
0241 
0242    % analytical solution
0243    R = ll / ww / hh / scale/ conduc + 2*z_contact/scale^2;
0244 
0245    V= current*R;
0246 %  fprintf('Solver %s: %f\n', 'analytic', V);
0247 
0248    mdl.solve = @fwd_solve_1st_order;
0249    mdl.system_mat = @system_mat_1st_order;
0250 
0251    img= mk_image(mdl, conduc);
0252 
0253    fsol= fwd_solve(img);
0254 %  fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0255    unit_test_cmp('test 2d vs analytic', fsol.meas, V, 1e-11);
0256 
0257    mdl.solve = @np_fwd_solve;
0258    mdl.system_mat = @np_calc_system_mat;
0259    mdl.misc.perm_sym= '{n}';
0260    img= mk_image(mdl, conduc);
0261 
0262    fsol= fwd_solve(img);
0263 %  fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0264    unit_test_cmp('np_fwd_solve vs analytic', fsol.meas, V, 1e-11);

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005