0001 function data =fwd_solve_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
0024
0025
0026 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0027
0028 if nargin == 1
0029 img= fwd_model;
0030 end
0031 fwd_model= img.fwd_model;
0032
0033 img = data_mapper(img);
0034 if ~ismember(img.current_params, supported_params)
0035 error('EIDORS:PhysicsNotSupported', '%s does not support %s', ...
0036 'FWD_SOLVE_2P5D_1ST_ORDER',img.current_params);
0037 end
0038
0039 img = convert_img_units(img, 'conductivity');
0040
0041 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0042
0043 k = [0 Inf]; try k = img.fwd_model.fwd_solve_2p5d_1st_order.k; end
0044 method = 'quadv'; try method = img.fwd_model.fwd_solve_2p5d_1st_order.method; end
0045
0046 img.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0047 gnd = fwd_model.gnd_node;
0048 img.fwd_model.system_mat_2p5d_1st_order.k = 0;
0049 img.fwd_model.system_mat_2p5d_1st_order.factory = 1;
0050 Sf = system_mat_2p5d_1st_order( img );
0051 if length(k) == 1
0052 v=2*potential_k(Sf(k), pp, gnd);
0053 else
0054 switch method
0055 case 'trapz'
0056
0057 n = 0;
0058 tol = 1e-8;
0059 k(isinf(k)) = 1/sqrt(tol);
0060 vf = zeros(pp.n_elec*pp.n_stim,length(k));
0061 for ki = k
0062 n = n + 1;
0063 vf(:,n) = potential_k(Sf(ki), pp, gnd);
0064 end
0065 v=2/pi*trapz(k,vf,2);
0066 case 'quadv'
0067
0068 tol = 1e-8;
0069 kend = min(1/sqrt(tol), max(k));
0070
0071
0072 v=2/pi*quadv(@(kk) potential_k(Sf(kk), pp, gnd), k(1), kend);
0073 case 'integral'
0074 tol = 1e-8;
0075 kend = min(1/sqrt(tol), max(k));
0076
0077
0078 v=2/pi*integral(@(kk) potential_k(Sf(kk), pp, gnd), k(1), kend, 'ArrayVAlued', true);
0079 otherwise
0080 error(['unrecognized method: ' method]);
0081 end
0082 end
0083
0084
0085 v2meas = get_v2meas(pp.n_elec, pp.n_stim, img.fwd_model.stimulation);
0086 data.meas = v2meas' * v;
0087 data.time= NaN;
0088 data.name= 'solved by fwd_solve_2p5d_1st_order';
0089 try; if img.fwd_solve.get_all_meas == 1
0090 data.volt = v(1:pp.n_node,:);
0091 end; end
0092 try; if img.fwd_solve.get_all_nodes== 1
0093 data.volt = v;
0094 end; end
0095
0096
0097 function v = potential_k(S, pp, gnd)
0098 vf_nodes = zeros(pp.n_node, pp.n_stim);
0099 idx = 1:size(S,1);
0100 idx( gnd ) = [];
0101 idx2 = find(any(pp.N2E));
0102 vf_nodes(idx,:) = left_divide( S(idx,idx), 1/2*pp.QQ(idx,:));
0103 vf_elecs = pp.N2E(:,idx2) * vf_nodes(idx2,:);
0104 v = vf_elecs(:);
0105
0106 function v2meas = get_v2meas(n_elec,n_stim,stim)
0107 v2meas = sparse(n_elec*n_stim,0);
0108 for i=1:n_stim
0109 meas_pat= stim(i).meas_pattern;
0110 n_meas = size(meas_pat,1);
0111 v2meas((i-1)*n_elec + 1: i*n_elec,end+(1:n_meas)) = meas_pat';
0112 end
0113
0114 function do_unit_test
0115 ne = 16;
0116 isOctave= exist('OCTAVE_VERSION');
0117
0118 imdl2 = mk_geophysics_model('h2a', ne);
0119 imdl2.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0120 img2 = mk_image(imdl2.fwd_model, 1);
0121
0122 img2.fwd_model.fwd_solve_2p5d_1st_order.k = [0 logspace(-2,0,12) 3];
0123 tic
0124 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0125 v2 = fwd_solve(img2);
0126 fprintf('fwd_solve 2D t= %f s\n', toc);
0127
0128 tic
0129 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'trapz';
0130 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0131 v2p5t = fwd_solve_2p5d_1st_order(img2);
0132 fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s, trapz)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0133 tic
0134 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'quadv';
0135 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0136 v2p5q = fwd_solve_2p5d_1st_order(img2);
0137 fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s, quadv)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0138 if isOctave
0139 fprintf('fwd_solve_2p5d_1st_order 2D t= --- s (-- k''s, integral) SKIPPED [not supported in Octave]\n');
0140 else
0141 tic
0142 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'integral';
0143 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0144 v2p5i = fwd_solve_2p5d_1st_order(img2);
0145 fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s, integral)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0146 end
0147 tic
0148 img2.fwd_model.fwd_solve_2p5d_1st_order.k = 0;
0149 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0150 v2p5k0 = fwd_solve_2p5d_1st_order(img2);
0151 fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0152
0153 imdl3 = mk_geophysics_model('h3a', ne);
0154 imdl3.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0155 img3 = mk_image(imdl3.fwd_model, 1);
0156 tic
0157 img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8;
0158 v3 = fwd_solve(img3);
0159 fprintf('fwd_solve 3D t= %f s\n', toc);
0160
0161 tic
0162 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0163 vh = fwd_solve_halfspace(img2);
0164 fprintf('fwd_solve_halfspace 3D t= %f s\n', toc);
0165 scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0166 unit_test_cmp('h2a 2p5d values TEST', uvh, ...
0167 [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0168 tol = norm(vh.meas)*0.025;
0169 unit_test_cmp('2D (h2a) vs analytic TEST', norm(vh.meas - v2.meas), 0, -tol);
0170 unit_test_cmp('2.5D (h2a + 2.5D, k=0) vs 2D TEST', norm(v2.meas - v2p5k0.meas), 0, tol);
0171 unit_test_cmp('2.5D (h2a + 2.5D trapz,2*tol) vs analytic TEST', norm(vh.meas - v2p5t.meas), 0, 2*tol);
0172 unit_test_cmp('2.5D (h2a + 2.5D quadv**) vs analytic TEST', norm(vh.meas - v2p5q.meas), 0, tol);
0173 if ~isOctave
0174 unit_test_cmp('2.5D (h2a + 2.5D integral) vs analytic TEST', norm(vh.meas - v2p5i.meas), 0, tol);
0175 end
0176 unit_test_cmp('3D (h3b) vs analytic TEST', norm(vh.meas - v3.meas), 0, 2*tol);
0177 pause(2);
0178 clf; h=plot([vh.meas, v2p5q.meas, v3.meas, v2.meas],':');
0179 legend('analytic', 'FEM 2.5D', 'FEM 3D', 'FEM 2D', 'Location','Best');
0180 set(gca,'box','off');
0181 set(h,'LineWidth',2);
0182 set(h(1),'Marker','o','MarkerSize',7);
0183 set(h(2),'Marker','square','MarkerSize',7);
0184 set(h(3),'Marker','diamond','MarkerSize',3);
0185 set(h(4),'Marker','+');
0186 pause(2);
0187 clf; h=plot([vh.meas, v2p5q.meas, v3.meas],':');
0188 legend('analytic', 'FEM 2.5D', 'FEM 3D','Location','Best');
0189 set(gca,'box','off');
0190 set(h,'LineWidth',2);
0191 set(h(1),'Marker','o','MarkerSize',7);
0192 set(h(2),'Marker','square','MarkerSize',7);
0193 set(h(3),'Marker','diamond','MarkerSize',3);