fwd_solve_higher_order

Solve for voltages (nodes/electrodes) for a forward model. 
M Crabb - 29.06.2012
 
 Example (2D):
  imdl = mk_common_model('c2C',16); img=mk_image(imdl.fwd_model,1);
  img.fwd_model.solve = @fwd_solve_higher_order;
  img.fwd_model.system_mat = @system_mat_higher_order;
  
  vve=[]; JJ4=[];
  for i= 1:3; switch i;
     case 1; img.fwd_model.approx_type = 'tri3'; % linear
     case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
     case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
     end %switch
     vv=fwd_solve(img);      vve(:,i)=vv.meas;
  end

 Example (3D):
  imdl = mk_common_model('b3cr',16);  img=mk_image(imdl.fwd_model,1);
  img.fwd_model.solve = @fwd_solve_higher_order;
  img.fwd_model.system_mat = @system_mat_higher_order;
  
  vve=[]; JJ4=[];
  for i= 1:2; switch i;
     case 1; img.fwd_model.approx_type = 'tet4'; % linear
     case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
     end %switch
     vv=fwd_solve(img);      vve(:,i)=vv.meas;
  end

0001 function[data] = fwd_solve_higher_order(fwd_model,img)
0002 %Solve for voltages (nodes/electrodes) for a forward model.
0003 %M Crabb - 29.06.2012
0004 %
0005 % Example (2D):
0006 %  imdl = mk_common_model('c2C',16); img=mk_image(imdl.fwd_model,1);
0007 %  img.fwd_model.solve = @fwd_solve_higher_order;
0008 %  img.fwd_model.system_mat = @system_mat_higher_order;
0009 %
0010 %  vve=[]; JJ4=[];
0011 %  for i= 1:3; switch i;
0012 %     case 1; img.fwd_model.approx_type = 'tri3'; % linear
0013 %     case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
0014 %     case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
0015 %     end %switch
0016 %     vv=fwd_solve(img);      vve(:,i)=vv.meas;
0017 %  end
0018 %
0019 % Example (3D):
0020 %  imdl = mk_common_model('b3cr',16);  img=mk_image(imdl.fwd_model,1);
0021 %  img.fwd_model.solve = @fwd_solve_higher_order;
0022 %  img.fwd_model.system_mat = @system_mat_higher_order;
0023 %
0024 %  vve=[]; JJ4=[];
0025 %  for i= 1:2; switch i;
0026 %     case 1; img.fwd_model.approx_type = 'tet4'; % linear
0027 %     case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
0028 %     end %switch
0029 %     vv=fwd_solve(img);      vve(:,i)=vv.meas;
0030 %  end
0031 
0032 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0033 
0034 if nargin == 1
0035    img= fwd_model;
0036 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0037    warning('EIDORS:DeprecatedInterface', ...
0038       ['Calling FWD_SOLVE_HIGHER_ORDER with two arguments is deprecated and will cause' ...
0039        ' an error in a future version. First argument ignored.']);
0040 end
0041 
0042 fwd_model= img.fwd_model;
0043 
0044 %Modify the forward model to be of my type
0045 if ~isfield(fwd_model,'approx_type')    || ...
0046    strcmp(fwd_model.approx_type,'tri3') || ...
0047    strcmp(fwd_model.approx_type,'tet4')   
0048     %Do nothing
0049 else
0050     [bound,elem,nodes] = fem_1st_to_higher_order(fwd_model);
0051     fwd_model.boundary=bound; fwd_model.elems=elem; fwd_model.nodes=nodes;
0052     %We need to update fwd_model of img too for system_mat
0053     img.fwd_model=fwd_model;
0054 end
0055 
0056 %Calculate the total stiffness matrix
0057 s_mat = calc_system_mat(img); At=s_mat.E;
0058 
0059 %Find electrode stucture and no.of electrodes and initialize vector
0060 %Find stim strucutre and no. stimulations
0061 %Find node structure and find no.nodes
0062 elecstruc=fwd_model.electrode; nelecs=size(elecstruc,2);
0063 stimstruc=fwd_model.stimulation; nstims=size(stimstruc,2); 
0064 nodestruc=fwd_model.nodes; nnodes=size(nodestruc,1); 
0065 
0066 %Find total number of measurements
0067 nmeass=0;
0068 for k=1:nstims
0069     stimkmeasmatrix = stimstruc(k).meas_pattern;
0070     nmeass=nmeass+size(stimkmeasmatrix,1);
0071 end
0072 
0073 %Complete or Point? Check first electrode (no mized models) and this changes
0074 %the index vector of what 'node' number corresponds to an electrode
0075 elecnode=zeros(1,nelecs);
0076 if(size(elecstruc(1).nodes,2)==1 && size(elecstruc(1).nodes,1)==1) %POINT
0077     for i=1:nelecs
0078         %POINT - Assign electrode index at correct node
0079         elecnode(i)=elecstruc(i).nodes;
0080     end
0081     %Assign correct size unknowns and right hand side matrix
0082     rhsdata=zeros(nnodes,nstims); 
0083     nodeunknowns=zeros(nnodes,nstims); 
0084 else
0085     for i=1:nelecs
0086         %COMPLETE - Assign electrode at bottom of list
0087         elecnode(i)=nnodes+i;
0088     end
0089     %Assign correct size right hand side matrix
0090     rhsdata=zeros(nnodes+nelecs,nstims); 
0091     nodeunknowns=zeros(nnodes+nelecs,nstims); 
0092 end
0093 
0094 %Assign currents at correct point in rhs matrix using index vector
0095 for ii=1:nstims
0096     %The vector of current values for stimulation
0097     curnode=stimstruc(ii).stim_pattern;
0098     for i=1:nelecs
0099         rhsdata(elecnode(i),ii)=curnode(i);
0100     end
0101 end
0102 
0103 %Create index vector and eliminate ground node equation from index
0104 groundnode=fwd_model.gnd_node; idx=1:size(At,1); idx(groundnode)=[];
0105 
0106 %Solve the simulated linear system with index
0107 nodeunknowns(idx,:)=left_divide(At(idx,idx),rhsdata(idx,:));
0108 
0109 
0110 %Find electrode voltages and store in matrix
0111 %Calculate electrode voltages using index vector elecnode(i)
0112 velec=zeros(nelecs,nstims);
0113 for i=1:nelecs
0114     %This is the indexed entries in nodeunknowns
0115     velec(i,:)=nodeunknowns(elecnode(i),:);
0116 end
0117 
0118 %Get the measured voltages
0119 vmeaselec=zeros(nmeass,1); idx=0;
0120 for ii=1:nstims
0121     % Use conj on meas_pat: see Chap 5 Adler&Holder 2021
0122     meas_pat=conj(stimstruc(ii).meas_pattern); %Measurement patterns
0123     n_meas=size(meas_pat,1); %Number of measures
0124     vmeaselec(idx + (1:n_meas) ) = meas_pat*velec(:,ii); %Diff data
0125     idx=idx+n_meas; %Increase counter
0126 end
0127 
0128 %Return the electrode voltages in data structure
0129 data.meas= vmeaselec;
0130 data.time= NaN; % unknown
0131 data.name= 'solved by fwd_solve_higher_order';
0132 data.quantity = 'voltage';
0133 try; if img.fwd_solve.get_all_meas == 1
0134    data.volt = nodeunknowns(1:nnodes,:); % but not on CEM nodes
0135 end; end
0136 try; if img.fwd_solve.get_all_nodes== 1
0137    data.volt = nodeunknowns;             % all, including CEM nodes
0138 end; end
0139 
0140 
0141 end
0142 
0143 function do_unit_test
0144    tol = 1e-13;
0145    vve = do_unit_test_2D;
0146    vve_ref = [
0147    0.898225115241117   0.921510761628436   0.928516596253320
0148    0.406398239528486   0.412406966923347   0.413938179536629
0149    0.248067950415984   0.250212111398160   0.250609888163540
0150    0.179592593985273   0.179643951982781   0.179734695991927];
0151 
0152    unit_test_cmp('2D: 1st order',vve(1:4,1),vve_ref(1:4,1),tol);
0153    unit_test_cmp('2D: 2nd order',vve(1:4,2),vve_ref(1:4,2),tol);
0154    unit_test_cmp('2D: 3rd order',vve(1:4,3),vve_ref(1:4,3),tol);
0155    vve = do_unit_test_3D;
0156    vve_ref = [
0157    1.404189968566952   1.410900674463290
0158    0.403207625837809   0.402992578774667
0159    0.198517193844915   0.201211971071238
0160    0.133852904079284   0.133841105904217];
0161    unit_test_cmp('3D: 1st order',vve(1:4,1),vve_ref(1:4,1),tol);
0162    unit_test_cmp('3D: 2nd order',vve(1:4,2),vve_ref(1:4,2),tol);
0163 
0164 end
0165 function [vve] = do_unit_test_2D
0166    imdl = mk_common_model('c2C',16); img = mk_image(imdl.fwd_model,1);
0167    vv=fwd_solve(img);      v0e=vv.meas;
0168 
0169    %High-order EIDORS solver %Change default eidors solvers
0170    img.fwd_model.solve = @fwd_solve_higher_order;
0171    img.fwd_model.system_mat = @system_mat_higher_order;
0172 
0173    vve=[]; JJ4=[];
0174    for i= 1:3; switch i;
0175       case 1; img.fwd_model.approx_type = 'tri3'; % linear
0176       case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
0177       case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
0178       end %switch
0179       vv=fwd_solve(img);      vve(:,i)=vv.meas;
0180    end
0181 
0182    subplot(321);
0183    plot([v0e,vve,(v0e*[1,1,1]-vve)*10]);
0184    legend('Default','linear','quadratic','cubic','(1-0)x10','(2-0)x10','(3-0)x10');
0185    xlim([1,100]);
0186 end
0187 function vve=do_unit_test_3D;
0188    imdl = mk_common_model('b3cr',16); img = mk_image(imdl.fwd_model,1);
0189    vv=fwd_solve(img);      v0e=vv.meas;
0190 
0191    %High-order EIDORS solver %Change default eidors solvers
0192    img.fwd_model.solve = @fwd_solve_higher_order;
0193    img.fwd_model.system_mat = @system_mat_higher_order;
0194 
0195    vve=[]; JJ4=[];
0196    for i= 1:2; switch i;
0197       case 1; img.fwd_model.approx_type = 'tet4'; % linear
0198       case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
0199       end %switch
0200       vv=fwd_solve(img);      vve(:,i)=vv.meas;
0201    end
0202 
0203    subplot(322);
0204    plot([v0e,vve,(v0e*[1,1]-vve)*10]);
0205    legend('Default','linear','quadratic','(1-0)x10','(2-0)x10');
0206    xlim([1,100]);
0207 
0208 end