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