demo_real_test

 Perform tests based on the demo_real function

0001 function ok= demo_real_test
0002 % Perform tests based on the demo_real function
0003 
0004 % (C) 2005 Andy Adler + Nick Polydorides. License: GPL version 2 or version 3
0005 % $Id: demo_real_test.m 1535 2008-07-26 15:36:27Z aadler $
0006 
0007 isOctave= exist('OCTAVE_VERSION');
0008 
0009 datareal= 'datareal.mat';
0010 datacom=  'datacom.mat';
0011 drt=      'demo_real_test.mat';
0012 if isOctave
0013     datareal= file_in_loadpath(datareal);
0014     datacom=  file_in_loadpath(datacom);
0015     drt    =  file_in_loadpath(drt);
0016     page_screen_output= 0;
0017 end
0018 
0019 load(datareal);
0020 perm_sym= '{y}';
0021 
0022 [I,Ib] = set_3d_currents(protocol,elec,vtx,gnd_ind,no_pl);
0023 
0024 mat_ref = 1*ones(828,1);
0025 
0026 %Set the tolerance for the forward solver
0027 tol = 1e-5;
0028 
0029 [Eref,D,Ela,ppr] = fem_master_full(vtx,simp,mat_ref,gnd_ind,elec,zc,perm_sym);
0030 [Vref] = forward_solver(Eref,I,tol,ppr);
0031 [refH,refV,indH,indV,dfr]=get_3d_meas(elec,vtx,Vref,Ib,no_pl);
0032 dfr = dfr(1:2:end); %Taking just the horrizontal measurements
0033 
0034 mat=mat_ref;
0035 load(datacom,'A','B'); %Indices of the elements to represent the inhomogeneity
0036 sA = mat_ref(A(1))+0.15;
0037 sB = mat_ref(B(1))-0.20;
0038 mat(A) = sA;
0039 mat(B) = sB;
0040 
0041 [En,D,Ela,ppn] = fem_master_full(vtx,simp,mat,gnd_ind,elec,zc,perm_sym);
0042 [Vn] = forward_solver(En,I,tol,ppn,Vref);
0043 [voltageH,voltageV,indH,indV,dfv]=get_3d_meas(elec,vtx,Vn,Ib,no_pl);
0044 dfv = dfv(1:2:end);
0045 
0046 if size(dfr)~= size(dfv)
0047    error('Mismatched measurements')
0048 end
0049 
0050 [v_f] = m_3d_fields(vtx,32,indH,Eref,tol,gnd_ind);
0051 
0052 [J] = jacobian_3d(I,elec,vtx,simp,gnd_ind,mat_ref,zc,v_f,dfr,tol,perm_sym);
0053 
0054 [Reg] = iso_f_smooth(simp,vtx,3,1);
0055 
0056 tfac = 1e-8;
0057 
0058 sol = (J'*J +  tfac*Reg'*Reg)\J' * (voltageH - refH);
0059 
0060 Diag_Reg_012= [diag(Reg,0),[diag(Reg,1);0],[diag(Reg,2);0;0]];
0061 Jcolsby100=J(:,1:100:size(J,2));
0062 
0063 %% verifications
0064 
0065 load(drt);
0066 
0067 compare_tol( drt.voltageH, voltageH, 'voltageH' )
0068 compare_tol( drt.voltageV, voltageV, 'voltageV' )
0069 compare_tol( drt.sol, sol, 'sol' )
0070 compare_tol( drt.Jcolsby100, Jcolsby100, 'Jcolsby100' )
0071 compare_tol( drt.Diag_Reg_012, Diag_Reg_012, 'Diag_Reg_012' )
0072 
0073 ok=1;
0074 
0075 
0076 function compare_tol( cmp1, cmp2, errtext )
0077 % compare matrices and give error if not equal
0078 fprintf(2,'testing parameter: %s ...\n',errtext);
0079 
0080 tol= 2e-4;
0081 
0082 vd= mean(mean( abs(cmp1 - cmp2) ));
0083 vs= mean(mean( abs(cmp1) + abs(cmp2) ));
0084 if vd/vs > tol
0085    eidors_msg( ...
0086      'parameter %s exceeds tolerance %g (=%g)', errtext, tol, vd/vs,1 );
0087 end
0088