inv_solve_dual_mesh

 INV_SOLVE_DUAL_MESH using a coarse and fine mesh
 img= inv_solve_dual_mesh( inv_model, data1)
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data1      => EIT data object

0001 function img= inv_solve_dual_mesh( inv_model, voltage)
0002 % INV_SOLVE_DUAL_MESH using a coarse and fine mesh
0003 % img= inv_solve_dual_mesh( inv_model, data1)
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data1      => EIT data object
0007 %
0008 
0009 % (C) 2005 David Stephenson. License: GPL version 2 or version 3
0010 % $Id: inv_solve_dual_mesh.m 3060 2012-06-06 16:37:38Z aadler $
0011 
0012 warning('EIDORS:deprecated','INV_SOLVE_DUAL_MESH is deprecated as of 06-Jun-2012. INV_SOLVE now does this.');
0013 
0014 M_dense= inv_model.fwd_model;
0015 % Load parameters
0016 M_coarse= inv_model.inv_solve_dual_mesh.coarse_mdl;
0017 %mapper_func= inv_model.inv_solve_dual_mesh.mapper_func;
0018 
0019 [index_simp]=edge_refined_elem_mapper( M_coarse, M_dense);
0020 c2d_elem= index_simp(:,1);
0021 [index_vtx]=edge_refined_node_mapper( M_coarse, M_dense);
0022 d2c_node= index_vtx(:,1);
0023 
0024 %FIXME: create parameters
0025 tol_for= 1e-5;
0026 tol_inv= 1e-5;
0027 
0028 tfac= calc_hyperparameter(inv_model);
0029 IM_coarse= inv_model; IM_coarse.fwd_model= M_coarse;
0030 RtR= calc_RtR_prior( IM_coarse );
0031 
0032 % Initial conductivity estimate
0033 mat_ref_coarse = ones(size(M_coarse.elems,1),1);
0034 mat_ref_coarse(:)= inv_model.jacobian_bkgnd.value;
0035 mat_ref_dense= mat_ref_coarse( c2d_elem );
0036 % Initial estimate - homogeneous background coarse mesh
0037 sol_upd_dense= mat_ref_dense;
0038 sol_upd_coarse = mat_ref_coarse;
0039 
0040 index_simp=edge_refined_elem_mapper(M_coarse,M_dense);
0041 
0042 % create dense model
0043 img_dense= eidors_obj('image', 'dense image');
0044 img_dense.elem_data= mat_ref_dense;
0045 img_dense.fwd_model= M_dense;
0046 pdense= np_fwd_parameters( M_dense );
0047 vtx_dense      = pdense.vtx;
0048 simp_dense     = pdense.simp;
0049 gnd_ind_dense  = pdense.gnd_ind;
0050 elec_dense     = pdense.elec;
0051 I_dense        = pdense.I;
0052 Ib_dense       = pdense.Ib;
0053 df_dense       = pdense.df;
0054 elec_dense     = pdense.elec;
0055 zc             = pdense.zc;
0056 perm_sym       = pdense.perm_sym;
0057 no_pl          = 1; % FIXME
0058 el_no          = pdense.n_elec;
0059 
0060 % create coarse model
0061 img_coarse= eidors_obj('image', 'coarse image');
0062 img_coarse.elem_data= mat_ref_coarse;
0063 img_coarse.fwd_model= M_coarse;
0064 pcoarse= np_fwd_parameters( M_coarse );
0065 vtx_coarse      = pcoarse.vtx;
0066 simp_coarse     = pcoarse.simp;
0067 gnd_ind_coarse  = pcoarse.gnd_ind;
0068 elec_coarse     = pcoarse.elec;
0069 I_coarse        = pcoarse.I;
0070 Ib_coarse       = pcoarse.Ib;
0071 df_coarse       = pcoarse.df;
0072 elec_coarse     = pcoarse.elec;
0073 
0074 %Set the tolerance for the forward solver
0075 for iter= 1:inv_model.parameters.max_iterations;
0076    [E_dense,D_dense,Ela_dense,pp_dense] = ...
0077    fem_master_full(vtx_dense,simp_dense,sol_upd_dense,gnd_ind_dense,elec_dense,zc,perm_sym);
0078 
0079    if iter==1
0080       %sprintf('Current fields for iteration %d',i)
0081       [V_dense] = forward_solver(E_dense,I_dense,tol_for,pp_dense);
0082       [viH,viV,indH_dense,indV,df] = get_3d_meas(elec_dense,vtx_dense,V_dense,Ib_dense,no_pl);
0083       dfv = df(1:2:end);
0084       vi = (viH);
0085       %sprintf('Measurement fields for iteration %d',i)
0086       [v_f_dense] = m_3d_fields(vtx_dense,el_no,indH_dense,E_dense,tol_for,gnd_ind_dense);
0087    else
0088       %sprintf('Current fields for iteration %d',i)
0089       [V_dense] = forward_solver(E_dense,I_dense,tol_for,pp_dense,V_dense);
0090       [viH,viV,indH_dense,indV,df] = get_3d_meas(elec_dense,vtx_dense,V_dense,Ib_dense,no_pl);
0091       dfv = df(1:2:end);
0092       vi = (viH);
0093       %sprintf('Measurement fields for iteration %d',i)
0094       [v_f_dense] = m_3d_fields(vtx_dense,el_no,indH_dense,E_dense,tol_for,gnd_ind_dense,v_f_dense);
0095    end
0096 
0097 
0098    %% DENSE MESH TO COARSE MESH HERE for scaler potentials
0099 
0100    [V_coarse] = mapvd(V_dense,index_vtx,elec_dense,vtx_coarse,vtx_dense);
0101 
0102    [v_f_coarse] = mapvm(v_f_dense,index_vtx,indH_dense,vtx_coarse,vtx_dense);
0103 
0104    [viH,viV,indH_coarse,indV,df_coarse] = get_3d_meas(elec_coarse,vtx_coarse,V_coarse,Ib_coarse,no_pl); % for dfv_coarse in Jacobian calculation
0105 
0106    dfv_coarse = df_coarse(1:2:end);
0107 
0108    [E_coarse,D_coarse,Ela_coarse,pp_coarse] = fem_master_full(vtx_coarse,simp_coarse,mat_ref_coarse,gnd_ind_coarse,elec_coarse,zc,perm_sym);  % for D in Jacobian calculation
0109 
0110    [J_coarse] = jacobian_3d_with_fields(V_coarse,Ela_coarse,D_coarse,I_coarse,elec_coarse,vtx_coarse,simp_coarse,gnd_ind_coarse,sol_upd_coarse,zc,v_f_coarse,dfv_coarse,tol_for,perm_sym);
0111 
0112    sol_coarse = (J_coarse.'*J_coarse + tfac^2*RtR)\ (J_coarse.' * (voltage - vi));
0113 
0114    % COARSE MESH TO DENSE MESH HERE for conductivity
0115 
0116    sol_dense=sol_coarse(c2d_elem);
0117 
0118    sol_upd_dense = sol_upd_dense + sol_dense;
0119    sol_upd_coarse = sol_upd_coarse + sol_coarse;
0120 
0121    sol_coarse = sol_upd_coarse;
0122    sol_dense = sol_upd_dense;
0123 
0124    sol_array(:,iter)=sol_coarse;
0125 
0126    sprintf('Error norm at iteration %d is %f',iter,norm(voltage - vi))
0127 
0128 end
0129 
0130 
0131 % create a data structure to return
0132 img.name= 'solved by inv_solve_trunc_iterative';
0133 img.elem_data = sol_array;
0134 img.fwd_model= M_dense;
0135 
0136 function voltH = elec_volts( Vfwd, fwd_model)
0137    p = np_fwd_parameters( fwd_model);
0138    Velec=Vfwd( p.n_node+(1:p.n_elec),:);
0139    voltH = zeros( p.n_meas, 1 );
0140    idx=0;
0141    for i=1:p.n_stim
0142       meas_pat= fwd_model.stimulation(i).meas_pattern;
0143       n_meas  = size(meas_pat,1);
0144       voltH( idx+(1:n_meas) ) = meas_pat*Velec(:,i);
0145       idx= idx+ n_meas;
0146    end
0147 
0148 function [V_coarse] = mapvd(V_dense,index_vtx,elec,vtx_coarse,vtx_dense);
0149 
0150 % This function maps the scaler potential field
0151 % from the dense mesh to the coarse
0152 % mesh using verticies extraction.
0153 
0154    for j=1:size(elec,1);
0155 
0156       for i=1:size(vtx_coarse,1);
0157 
0158          V_coarse(i,j)=V_dense((index_vtx(i,1)),j);
0159 
0160       end
0161 
0162    end
0163 
0164    V_coarse_elec=V_dense((size(vtx_dense,1)+1):(size(V_dense,1)),:);
0165 
0166    V_coarse=[V_coarse;V_coarse_elec];
0167 
0168 % This function maps the scaler potential field
0169 % from the dense mesh to the coarse
0170 % mesh using verticies extraction.
0171 function [v_f_coarse] = mapvm(v_f_dense,index_vtx,indH_dense,vtx_coarse,vtx_dense);
0172 
0173 
0174    for j=1:size(indH_dense,1);
0175 
0176       for i=1:size(vtx_coarse,1);
0177 
0178          v_f_coarse(i,j)=v_f_dense((index_vtx(i,1)),j);
0179 
0180       end
0181 
0182    end
0183 
0184    v_f_coarse_elec=v_f_dense((size(vtx_dense,1)+1):(size(v_f_dense,1)),:);
0185 
0186    v_f_coarse=[v_f_coarse;v_f_coarse_elec];