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