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 isOctave= exist('OCTAVE_VERSION');
0108
0109 imdl2 = mk_geophysics_model('h2a', ne);
0110 imdl2.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0111 img2 = mk_image(imdl2.fwd_model, 1);
0112
0113 img2.fwd_model.fwd_solve_2p5d_1st_order.k = [0 logspace(-2,0,12) 3];
0114 tic
0115 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0116 v2 = fwd_solve(img2);
0117 fprintf('fwd_solve 2D t= %f s\n', toc);
0118
0119 tic
0120 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'trapz';
0121 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0122 v2p5t = fwd_solve_2p5d_1st_order(img2);
0123 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));
0124 tic
0125 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'quadv';
0126 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0127 v2p5q = fwd_solve_2p5d_1st_order(img2);
0128 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));
0129 if isOctave
0130 fprintf('fwd_solve_2p5d_1st_order 2D t= --- s (-- k''s, integral) SKIPPED [not supported in Octave]\n');
0131 else
0132 tic
0133 img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'integral';
0134 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0135 v2p5i = fwd_solve_2p5d_1st_order(img2);
0136 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));
0137 end
0138 tic
0139 img2.fwd_model.fwd_solve_2p5d_1st_order.k = 0;
0140 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0141 v2p5k0 = fwd_solve_2p5d_1st_order(img2);
0142 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));
0143
0144 imdl3 = mk_geophysics_model('h3a', ne);
0145 imdl3.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0146 img3 = mk_image(imdl3.fwd_model, 1);
0147 tic
0148 img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8;
0149 v3 = fwd_solve(img3);
0150 fprintf('fwd_solve 3D t= %f s\n', toc);
0151
0152 tic
0153 img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8;
0154 vh = fwd_solve_halfspace(img2);
0155 fprintf('fwd_solve_halfspace 3D t= %f s\n', toc);
0156 scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0157 unit_test_cmp('h2a 2p5d values TEST', uvh, ...
0158 [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0159 tol = norm(vh.meas)*0.025;
0160 unit_test_cmp('2D (h2a) vs analytic TEST', norm(vh.meas - v2.meas), 0, -tol);
0161 unit_test_cmp('2.5D (h2a + 2.5D, k=0) vs 2D TEST', norm(v2.meas - v2p5k0.meas), 0, tol);
0162 unit_test_cmp('2.5D (h2a + 2.5D trapz,2*tol) vs analytic TEST', norm(vh.meas - v2p5t.meas), 0, 2*tol);
0163 unit_test_cmp('2.5D (h2a + 2.5D quadv**) vs analytic TEST', norm(vh.meas - v2p5q.meas), 0, tol);
0164 if ~isOctave
0165 unit_test_cmp('2.5D (h2a + 2.5D integral) vs analytic TEST', norm(vh.meas - v2p5i.meas), 0, tol);
0166 end
0167 unit_test_cmp('3D (h3b) vs analytic TEST', norm(vh.meas - v3.meas), 0, 2*tol);
0168 pause(2);
0169 clf; h=plot([vh.meas, v2p5q.meas, v3.meas, v2.meas],':');
0170 legend('analytic', 'FEM 2.5D', 'FEM 3D', 'FEM 2D', 'Location','Best');
0171 set(gca,'box','off');
0172 set(h,'LineWidth',2);
0173 set(h(1),'Marker','o','MarkerSize',7);
0174 set(h(2),'Marker','square','MarkerSize',7);
0175 set(h(3),'Marker','diamond','MarkerSize',3);
0176 set(h(4),'Marker','+');
0177 pause(2);
0178 clf; h=plot([vh.meas, v2p5q.meas, v3.meas],':');
0179 legend('analytic', 'FEM 2.5D', 'FEM 3D','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);