Description of system_mat_2p5d_1st

0001 function s_mat= system_mat_2p5d_1st_order( fwd_model, img)
0002 % SYSTEM_MAT_2P5D_1ST_ORDER: 2.5D system matrix
0003 % SS= SYSTEM_MAT_2P5D_1ST_ORDER( FWD_MODEL, IMG)
0004 % fwd_model = forward model
0005 % img       = image background for system matrix calc
0006 % s_mat = CC' * (SS + k^2 * TT) * conductivites * CC;
0007 % where:
0008 %   SS  = Unconnected system Matrix
0009 %   TT  = Unconnected Fourier correction
0010 %   CC  = Connectivity Matrix
0011 %
0012 % Parameters
0013 %   fwd_model.system_mat_2_5D_1st_order.k
0014 %        the k-th harmonic 0..inf [default 0]
0015 %   fwd_model.system_mat_2_5D_1st_order.factory
0016 %        return a factory in s_mat.assemble(s_mat, k) to build
0017 %        the system matrix based on argument k [default 0]
0018 %
0019 % NOTE: This code is adapted from dev/a_adler/system_mat_2_5D_1st_order.m
0020 
0021 % (C) 2005, 2013, 2016 A. Adler, A. Boyle. Licensed under the GPL Version 2
0022 
0023 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 if nargin == 1
0026    img= fwd_model;
0027 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0028    warning('EIDORS:DeprecatedInterface', ...
0029       ['Calling SYSTEM_MAT_1ST_ORDER with two arguments is deprecated and will cause' ...
0030        ' an error in a future version. First argument ignored.']);
0031 end
0032 fwd_model= img.fwd_model;
0033 
0034 % check parametrizations
0035 if isfield(img,'current_params') && ...
0036      (~ischar(img.current_params) || ~strcmp(img.current_params,'conductivity'))
0037     if ischar(img.current_params)
0038        error('system_mat_2p5d_1st_order does not work for %s',img.current_params);
0039     else
0040        img.current_params
0041        error('system_mat_2p5d_1st_order does not work for the above');
0042     end
0043 end
0044 
0045 factory = 0; try factory = fwd_model.system_mat_2p5d_1st_order.factory; end
0046 k = 0; try k = fwd_model.system_mat_2p5d_1st_order.k; end
0047 
0048 FS = system_mat_fields(img.fwd_model);      % 2D solution
0049 FT = system_mat_2p5d_fields(img.fwd_model); % Fourier domain correction
0050 [ES, ET] = calc_sigma(img, FS, FT); % conductivity
0051 
0052 if factory
0053    s_mat = @(k) assemble(k,FS,FT,ES,ET);
0054 else
0055    s_mat.E = assemble(k,FS,FT,ES,ET);
0056 end
0057 
0058 
0059 function SS = assemble(k,FS,FT,ES,ET)
0060    assert(all(size(k) == [1 1]), 'expected scalar k');
0061    SS = FS'*ES*FS + k^2*FT'*ET*FT;
0062 
0063 function [ES,ET] = calc_sigma(img, FS, FT)
0064    lFS= size(FS,1);
0065    lFT= size(FT,1);
0066    ED = elem_dim(img.fwd_model);
0067    lNE= ED*num_elems(img.fwd_model);
0068    elem_data = check_elem_data(img.fwd_model, img);
0069    if size(elem_data,3) == 1
0070    % Scalar conductivity == isotropic
0071       elem_sigma = kron( elem_data, ones(ED,1) );
0072       elem_sigma(end+1:lFS) = 1; % add ones for CEM
0073       ES= spdiags(elem_sigma,0,lFS,lFS);
0074 
0075       elem_sigma = kron( elem_data, ones(ED+1,1) );
0076       elem_sigma(end+1:lFT) = 1; % add ones for CEM
0077       ET= spdiags(elem_sigma,0,lFT,lFT);
0078    else
0079 %      switch elem_dim(img.fwd_model)
0080 %        case 2;
0081 %          idx = 1:2:lNE;
0082 %          ES= sparse([idx,idx+1,idx,idx+1]', ...
0083 %                     [idx,idx,idx+1,idx+1]', elem_data(:), lFS,lFS);
0084 %
0085 %          ES(lNE+1:lFS,lNE+1:lFS) = speye(lFS-lNE);
0086 %        case 3;
0087 %          idx = 1:3:lNE;
0088 %          ES= sparse([idx,idx+1,idx+2,idx,idx+1,idx+2,idx,idx+1,idx+2]', ...
0089 %                     [idx,idx,idx,idx+1,idx+1,idx+1,idx+2,idx+2,idx+2]', ...
0090 %                      elem_data(:), lFS,lFS);
0091 %
0092 %          ES(lNE+1:lFS,lNE+1:lFS) = speye(lFS-lNE);
0093 %        otherwise;
0094           error('%d D anisotropic elements not implemented', elem_dim(fwd_model));
0095 %      end
0096    end
0097 
0098 function check_2d(img)
0099    assert(size(img.fwd_model.nodes,2) == 2, 'expecting 2D model for 2.5D reconstruction');
0100    % check that all electrodes are actually part of the mesh
0101    nn = unique(img.fwd_model.elems);
0102    for i=1:length(img.fwd_model.electrode)
0103       ne = img.fwd_model.electrode(i).nodes;
0104       for j=1:length(ne)
0105          assert(length(find(ne(j) == nn))>0, sprintf('elec#%d, node#%d: not connected to mesh',i,j));
0106       end
0107    end
0108 
0109 function elem_data = check_elem_data(fwd_model, img);
0110    elem_data = img.elem_data;
0111    if any(size(elem_data) == [1 1])
0112       elem_data = elem_data(:);
0113    end
0114    sz_elem_data = size(elem_data);
0115    if sz_elem_data(2) ~= 1;
0116       error('system_mat_2p5d_1st_order: can only solve one image (sz_elem_data=%)', ...
0117             sz_elem_data);
0118    end
0119 
0120    if isfield(fwd_model, 'coarse2fine');
0121      c2f = fwd_model.coarse2fine;
0122      sz_c2f = size(c2f);
0123      switch sz_elem_data(1)
0124        case sz_c2f(1); % Ok
0125        case sz_c2f(2); elem_data = c2f * elem_data;
0126           if isfield(fwd_model, 'background')
0127               elem_data = elem_data + fwd_model.background;
0128           end
0129 
0130        otherwise; error(['system_mat_2p5d_1st_order: provided elem_data ' ...
0131             ' (sz=%d) does not match c2f (sz=%d %d)'], sz_elem_data(1), sz_c2f);
0132      end
0133    else
0134      if sz_elem_data(1) ~= num_elems(fwd_model)
0135        error(['system_mat_2p5d_1st_order: provided elem_data (sz=%d) does ' ...
0136           ' not match fwd_model (sz=%d)'], sz_elem_data(1), num_elems(fwd_model));
0137      end
0138    end
0139 
0140 function do_unit_test()
0141    ne = 16;
0142    k=0; % for k=0, the result is a 2D reconstruction
0143    % 2D
0144    imdl2 = mk_geophysics_model('h2a', ne);
0145    imdl2.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0146    img2 = mk_image(imdl2.fwd_model, 1);
0147    check_2d(img2);
0148    v2 = fwd_solve(img2);
0149    % 2.5D
0150    img2.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0151    img2.fwd_model.system_mat_2p5d_1st_order.k = k;
0152    v2p5 = fwd_solve(img2);
0153    % 3D
0154    imdl3 = mk_geophysics_model('h3a', ne);
0155    imdl3.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0156    img3 = mk_image(imdl3.fwd_model, 1);
0157    v3 = fwd_solve(img3);
0158    % analytic half-space (dipole)
0159    vh = fwd_solve_halfspace(img2);
0160    scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0161    unit_test_cmp('h2a 2p5d values TEST', uvh, ...
0162    [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0163    tol = norm(vh.meas)*0.050; % 5% error tolerance
0164    unit_test_cmp('2D   (h2a)                vs analytic TEST', norm(vh.meas - v2.meas),   0, -tol);
0165    unit_test_cmp('2.5D (h2a + 2.5D sys_mat) vs 2D       TEST', norm(v2.meas - v2p5.meas), 0, tol);
0166    unit_test_cmp('3D   (h3a)                vs analytic TEST', norm(vh.meas - v3.meas),   0, tol);
0167 clf; h=plot([vh.meas, v2.meas, v2p5.meas, v3.meas],':');
0168      legend('analytic','FEM 2D', sprintf('FEM 2.5D k=%0.1g',k), 'FEM 3D','Location','Best');
0169      set(gca,'box','off');
0170      set(h,'LineWidth',2);
0171      set(h(1),'Marker','o','MarkerSize',7);
0172      set(h(2),'Marker','+');
0173      set(h(3),'Marker','square','MarkerSize',7);
0174      set(h(4),'Marker','diamond','MarkerSize',3);
system_mat_2p5d_1st_order

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
