eidors2d_demo1

EidorsDemo1 Demonstrates the use of 2D EIT Package with linear basis
 EidorsDemo1 Demonstrates the use of 2D EIT Package for simulations with linear approksimation basis

0001 %EidorsDemo1 Demonstrates the use of 2D EIT Package with linear basis
0002 % EidorsDemo1 Demonstrates the use of 2D EIT Package for simulations with linear approksimation basis
0003 
0004 % M. Vauhkonen 28.3.2000,
0005 % University of Kuopio, Department of Applied Physics, PO Box 1627,
0006 % FIN-70211 Kuopio, Finland, email: Marko.Vauhkonen@uku.fi
0007 
0008 % Refactored for EIDORS 3.4 by Andy Adler - June 2009
0009 % $Id: eidors2d_demo1.m 3087 2012-06-07 18:11:42Z aadler $
0010 
0011 
0012 tgt_elems= [374,375,376,601,603,604, ...
0013             250,254,268,437,449,456];
0014 
0015 [fmdl1,fmdl2] = mv_mdl_meshdata;
0016 
0017 % Stimulation
0018 stim = mk_stim_patterns(16,1,'{trig}','{ad}',{},1);
0019 
0020 % Override with MV's stim pattern
0021 % Trigonometric current pattern.
0022 [jnk,T]=Current(length(fmdl1.electrode),0,'tri');
0023 
0024 for i=1:size(T,2)
0025    stim(i).stim_pattern= T(:,i);
0026    stim(i).meas_pattern(16,:) = [];
0027 end
0028 %stim(2:end) = [];
0029 
0030 fmdl2.stimulation = stim;
0031 fmdl1.stimulation = stim;
0032 
0033 Ne2= size(fmdl2.elems,1);
0034 Ne2= size(fmdl2.elems,1);
0035 
0036 
0037 % Create a sample image
0038 elem_data = 1/400*ones(Ne2,1);
0039 elem_data(tgt_elems) = 1/200;
0040 
0041 tgt_img= mk_image(fmdl2,elem_data);
0042 
0043 show_fem(tgt_img,[0,1,0])
0044 
0045 tgt_img.fwd_model.system_mat = @system_mat_1st_order;
0046 tgt_img.fwd_model.solve = @fwd_solve_1st_order;
0047  meas_aa = fwd_solve( tgt_img );
0048 
0049 pp= fwd_model_parameters( tgt_img.fwd_model );
0050 s_mat= calc_system_mat(tgt_img.fwd_model, tgt_img );
0051 [tgt_img.fwd_model.electrode.z_contact]= deal(50);
0052 v= zeros(pp.n_node,pp.n_stim);
0053 
0054 idx= 1:size(s_mat.E,1); idx( tgt_img.fwd_model.gnd_node ) = [];
0055 
0056 tol= 1e-5;
0057 v(idx,:)= left_divide( s_mat.E(idx,idx), pp.QQ(idx,:), tol);
0058 
0059 
0060 fmdl2.system_mat = @mv_calc_system_mat;
0061 fmdl2.solve = @mv_fwd_solve;
0062 tgt_img.fwd_model= fmdl2;
0063 meas_mv = fwd_solve( tgt_img );
0064 
0065 for i=1:10
0066 plot([v(1:1049,i),meas_mv.U.MeasField(:,i)])
0067 pause
0068 end
0069 
0070 return
0071 
0072 load meshdata % Data for two different meshes.
0073 
0074 NNode1=max(size(Node1));                      %The number of nodes
0075 NElement1=max(size(Element1));                %The number of element
0076 NNode2=max(size(Node2));                      %The number of nodes
0077 NElement2=max(size(Element2));                %The number of elements
0078 
0079 g1=reshape([Node1.Coordinate],2,NNode1)';
0080 H1=reshape([Element1.Topology],3,NElement1)';
0081 g2=reshape([Node2.Coordinate],2,NNode2)';
0082 H2=reshape([Element2.Topology],3,NElement2)';
0083 
0084 
0085 %disp('Choose a circular inhomogeneity. Left mouse button, center, right button, radius.')
0086 %Ind=ChooseCircle(Node2,Element2);  % Make data for an inhomogeneity.
0087 Ind = tgt_elems;
0088 sigma=1/400*ones(NElement2,1);            % Make a conductivity vector.
0089 sigma(Ind)=2/400;              % Conductivity of the inhomogeneity.
0090 
0091 clf,Plotinvsol(1./sigma,g2,H2);colorbar,title(['Your resistivity distribution']);drawnow
0092 disp('Press any key to continue...'),pause
0093 
0094 disp('Computes the simulated data.')
0095 L=16;                      % The number of electrodes.
0096 z=0.005*ones(L,1);              % Contact impedances.
0097 [II1,T]=Current(L,NNode2,'tri');      % Trigonometric current pattern.
0098 
0099 [Agrad,Kb,M,S,C]=FemMatrix(Node2,Element2,z);
0100 A=UpdateFemMatrix(Agrad,Kb,M,S,sigma);  % The system matrix.
0101 
0102 [U,p,r]=ForwardSolution(NNode2,NElement2,A,C,T,[],'real'); % Simulated data.
0103 Uel=U.Electrode(:);
0104 
0105 Agrad1=Agrad*Ind2;   % Group some of the element for the inverse computations
0106 
0107 
0108 %%             PROCEDURE TO SOLVE THE INVERSE PROBLEM           %%
0109 
0110 % Approximate the best homogenous resistivity.
0111 
0112 
0113 disp('Solves the full nonlinear inverse problem by regularised Gauss-Newton iteration.')
0114 
0115 disp('Initialisations...')
0116 
0117 A=UpdateFemMatrix(Agrad,Kb,M,S,ones(NElement2,1));  % The system matrix.
0118 Uref=ForwardSolution(NNode2,NElement2,A,C,T,[],'real',p,r);
0119 
0120 rho0=Uref.Electrode(:)\U.Electrode(:);
0121 
0122 A=UpdateFemMatrix(Agrad,Kb,M,S,1./rho0*ones(size(sigma)));  % The system matrix.
0123 Uref=ForwardSolution(NNode2,NElement2,A,C,T,[],'real',p,r);
0124 Urefel=Uref.Electrode(:);
0125 
0126 rho=rho0*ones(size(Agrad1,2),1);
0127 J=Jacobian(Node2,Element2,Agrad1,Uref.Current,Uref.MeasField, ...
0128            rho,'real');
0129 
0130 %Regularisation parameter and matrix
0131 
0132 alpha = 0.000005; 
0133 R=MakeRegmatrix(Element1);
0134 
0135 iter=5;
0136 
0137 disp('Iterations...')
0138 
0139 for ii=1:iter
0140  rho=rho+(J'*J+alpha*R'*R)\(J'*(Uel-Urefel)-alpha*R'*R*rho);
0141  rhobig=Ind2*rho;
0142  A=UpdateFemMatrix(Agrad,Kb,M,S,1./rhobig);  % The system matrix.
0143  Uref=ForwardSolution(NNode2,NElement2,A,C,T,[],'real',p,r);
0144  Urefel=Uref.Electrode(:);
0145  J=Jacobian(Node2,Element2,Agrad1,Uref.Current,Uref.MeasField,rho,'real');
0146  clf,Plotinvsol(rho,g1,H1);colorbar,title([num2str(ii) '. step']);drawnow;
0147 end
0148 
0149 
0150 
0151 
0152 
0153 
0154 
0155 
0156 
0157