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     meas_pat=stimstruc(ii).meas_pattern; %Measurement patterns
0122     n_meas=size(meas_pat,1); %Number of measures
0123     vmeaselec(idx + (1:n_meas) ) = meas_pat*velec(:,ii); %Diff data
0124     idx=idx+n_meas; %Increase counter
0125 end
0126 
0127 %Return the electrode voltages in data structure
0128 data.meas= vmeaselec;
0129 data.time= NaN; % unknown
0130 data.name= 'solved by fwd_solve_higher_order';
0131 data.quantity = 'voltage';
0132 try; if img.fwd_solve.get_all_meas == 1
0133    data.volt = nodeunknowns(1:nnodes,:); % but not on CEM nodes
0134 end; end
0135 try; if img.fwd_solve.get_all_nodes== 1
0136    data.volt = nodeunknowns;             % all, including CEM nodes
0137 end; end
0138 
0139 
0140 end
0141 
0142 function do_unit_test
0143    tol = 1e-13;
0144    vve = do_unit_test_2D;
0145    vve_ref = [
0146    0.898225115241117   0.921510761628436   0.928516596253320
0147    0.406398239528486   0.412406966923347   0.413938179536629
0148    0.248067950415984   0.250212111398160   0.250609888163540
0149    0.179592593985273   0.179643951982781   0.179734695991927];
0150 
0151    unit_test_cmp('2D: 1st order',vve(1:4,1),vve_ref(1:4,1),tol);
0152    unit_test_cmp('2D: 2nd order',vve(1:4,2),vve_ref(1:4,2),tol);
0153    unit_test_cmp('2D: 3rd order',vve(1:4,3),vve_ref(1:4,3),tol);
0154    vve = do_unit_test_3D;
0155    vve_ref = [
0156    1.404189968566952   1.410900674463290
0157    0.403207625837809   0.402992578774667
0158    0.198517193844915   0.201211971071238
0159    0.133852904079284   0.133841105904217];
0160    unit_test_cmp('3D: 1st order',vve(1:4,1),vve_ref(1:4,1),tol);
0161    unit_test_cmp('3D: 2nd order',vve(1:4,2),vve_ref(1:4,2),tol);
0162 
0163 end
0164 function [vve] = do_unit_test_2D
0165    imdl = mk_common_model('c2C',16); img = mk_image(imdl.fwd_model,1);
0166    vv=fwd_solve(img);      v0e=vv.meas;
0167 
0168    %High-order EIDORS solver %Change default eidors solvers
0169    img.fwd_model.solve = @fwd_solve_higher_order;
0170    img.fwd_model.system_mat = @system_mat_higher_order;
0171 
0172    vve=[]; JJ4=[];
0173    for i= 1:3; switch i;
0174       case 1; img.fwd_model.approx_type = 'tri3'; % linear
0175       case 2; img.fwd_model.approx_type = 'tri6'; % quadratic
0176       case 3; img.fwd_model.approx_type = 'tri10'; % cubic;
0177       end %switch
0178       vv=fwd_solve(img);      vve(:,i)=vv.meas;
0179    end
0180 
0181    subplot(321);
0182    plot([v0e,vve,(v0e*[1,1,1]-vve)*10]);
0183    legend('Default','linear','quadratic','cubic','(1-0)x10','(2-0)x10','(3-0)x10');
0184    xlim([1,100]);
0185 end
0186 function vve=do_unit_test_3D;
0187    imdl = mk_common_model('b3cr',16); img = mk_image(imdl.fwd_model,1);
0188    vv=fwd_solve(img);      v0e=vv.meas;
0189 
0190    %High-order EIDORS solver %Change default eidors solvers
0191    img.fwd_model.solve = @fwd_solve_higher_order;
0192    img.fwd_model.system_mat = @system_mat_higher_order;
0193 
0194    vve=[]; JJ4=[];
0195    for i= 1:2; switch i;
0196       case 1; img.fwd_model.approx_type = 'tet4'; % linear
0197       case 2; img.fwd_model.approx_type = 'tet10'; % quadratic
0198       end %switch
0199       vv=fwd_solve(img);      vve(:,i)=vv.meas;
0200    end
0201 
0202    subplot(322);
0203    plot([v0e,vve,(v0e*[1,1]-vve)*10]);
0204    legend('Default','linear','quadratic','(1-0)x10','(2-0)x10');
0205    xlim([1,100]);
0206 
0207 end