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 data.meas = pp.v2meas * v;
0086 data.time= NaN;
0087 data.name= 'solved by fwd_solve_2p5d_1st_order';
0088 try; if img.fwd_solve.get_all_meas == 1
0089 data.volt = v(1:pp.n_node,:);
0090 end; end
0091 try; if img.fwd_solve.get_all_nodes== 1
0092 data.volt = v;
0093 end; end
0094
0095
0096 function v = potential_k(S, pp, gnd)
0097 vf_nodes = zeros(pp.n_node, pp.n_stim);
0098 idx = 1:size(S,1);
0099 idx( gnd ) = [];
0100 idx2 = find(any(pp.N2E));
0101 vf_nodes(idx,:) = left_divide( S(idx,idx), 1/2*pp.QQ(idx,:));
0102 vf_elecs = pp.N2E(:,idx2) * vf_nodes(idx2,:);
0103 v = vf_elecs(:);
0104
0105 function do_unit_test
0106 ne = 16;
0107
0108 imdl2 = mk_geophysics_model('h2a', ne);
0109 imdl2.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0110 img2 = mk_image(imdl2.fwd_model, 1);
0111
0112 img2.fwd_model.fwd_solve_2p5d_1st_order.k = [0 logspace(-2,0,12) 3];
0113 tic
0114 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0115 v2 = fwd_solve(img2);
0116 fprintf('fwd_solve 2D t= %f s\n', toc);
0117
0118 tic
0119 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'trapz';
0120 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0121 v2p5t = fwd_solve_2p5d_1st_order(img2);
0122 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));
0123 tic
0124 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'quadv';
0125 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0126 v2p5q = fwd_solve_2p5d_1st_order(img2);
0127 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));
0128
0129 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'integral';
0130 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0131 v2p5i = fwd_solve_2p5d_1st_order(img2);
0132 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));
0133 img2.fwd_model.fwd_solve_2p5d_1st_order.k = 0;
0134 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0135 v2p5k0 = fwd_solve_2p5d_1st_order(img2);
0136 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));
0137
0138 imdl3 = mk_geophysics_model('h3a', ne);
0139 imdl3.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0140 img3 = mk_image(imdl3.fwd_model, 1);
0141 tic
0142 img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8;
0143 v3 = fwd_solve(img3);
0144 fprintf('fwd_solve 3D t= %f s\n', toc);
0145
0146 tic
0147 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0148 vh = fwd_solve_halfspace(img2);
0149 fprintf('fwd_solve_halfspace 3D t= %f s\n', toc);
0150 scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0151 unit_test_cmp('h2a 2p5d values TEST', uvh, ...
0152 [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0153 tol = norm(vh.meas)*0.025;
0154 unit_test_cmp('2D (h2a) vs analytic TEST', norm(vh.meas - v2.meas), 0, -tol);
0155 unit_test_cmp('2.5D (h2a + 2.5D, k=0) vs 2D TEST', norm(v2.meas - v2p5k0.meas), 0, tol);
0156 unit_test_cmp('2.5D (h2a + 2.5D trapz,2*tol) vs analytic TEST', norm(vh.meas - v2p5t.meas), 0, 2*tol);
0157 unit_test_cmp('2.5D (h2a + 2.5D quadv**) vs analytic TEST', norm(vh.meas - v2p5q.meas), 0, tol);
0158 unit_test_cmp('2.5D (h2a + 2.5D integral) vs analytic TEST', norm(vh.meas - v2p5i.meas), 0, tol);
0159 unit_test_cmp('3D (h3b) vs analytic TEST', norm(vh.meas - v3.meas), 0, 2*tol);
0160 pause(2);
0161 clf; h=plot([vh.meas, v2p5q.meas, v3.meas, v2.meas],':');
0162 legend('analytic', 'FEM 2.5D', 'FEM 3D', 'FEM 2D', 'Location','NorthEast');
0163 set(gca,'box','off');
0164 set(h,'LineWidth',2);
0165 set(h(1),'Marker','o','MarkerSize',7);
0166 set(h(2),'Marker','square','MarkerSize',7);
0167 set(h(3),'Marker','diamond','MarkerSize',3);
0168 set(h(4),'Marker','+');
0169 pause(2);
0170 clf; h=plot([vh.meas, v2p5q.meas, v3.meas],':');
0171 legend('analytic', 'FEM 2.5D', 'FEM 3D','Location','NorthEast');
0172 set(gca,'box','off');
0173 set(h,'LineWidth',2);
0174 set(h(1),'Marker','o','MarkerSize',7);
0175 set(h(2),'Marker','square','MarkerSize',7);
0176 set(h(3),'Marker','diamond','MarkerSize',3);