0001 function s_mat= system_mat_2p5d_1st_order( fwd_model, img)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
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
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);
0049 FT = system_mat_2p5d_fields(img.fwd_model);
0050 [ES, ET] = calc_sigma(img, FS, FT);
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
0071 elem_sigma = kron( elem_data, ones(ED,1) );
0072 elem_sigma(end+1:lFS) = 1;
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;
0077 ET= spdiags(elem_sigma,0,lFT,lFT);
0078 else
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094 error('%d D anisotropic elements not implemented', elem_dim(fwd_model));
0095
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
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);
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;
0143
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
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
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
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;
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);