function [v_f] = m_3d_fields(vtx,el_no,m_ind,E,tol,gnd_ind,v_f); This function calculates the measurement fields using preconditioned conjugate gradients. vtx = The vertices el_no = The total number of electrodes in the system m_ind = The measurements matrix (indices of electrode pairs) E = The full rank system matrix tol = The tolerance in the forward solution gnd_ind = The ground index v_f = The measurements fields
0001 function [v_f] = m_3d_fields(vtx,el_no,m_ind,E,tol,gnd_ind,v_f); 0002 %function [v_f] = m_3d_fields(vtx,el_no,m_ind,E,tol,gnd_ind,v_f); 0003 % 0004 %This function calculates the measurement fields using preconditioned conjugate gradients. 0005 % 0006 % 0007 % 0008 %vtx = The vertices 0009 %el_no = The total number of electrodes in the system 0010 %m_ind = The measurements matrix (indices of electrode pairs) 0011 %E = The full rank system matrix 0012 %tol = The tolerance in the forward solution 0013 %gnd_ind = The ground index 0014 %v_f = The measurements fields 0015 0016 %FIXME: This code should call forward_solver, it can then decide 0017 %what the best solver strategy is. Right now cgls is slower than \ 0018 0019 % (C) Nick Polydorides GPL v2 or v3. $Id: m_3d_fields.html 2819 2011-09-07 16:43:11Z aadler $ 0020 0021 0022 [vr,vc] = size(vtx); 0023 0024 Is_supl = zeros(vr,size(m_ind,1)); 0025 %no of electrodes x no of measurements (now currents)! 0026 0027 MC = []; 0028 0029 for i=1:size(m_ind,1) 0030 0031 m_n = zeros(el_no,1); 0032 0033 m_n(m_ind(i,1)) = 1; 0034 m_n(m_ind(i,2)) = -1; 0035 0036 MC = [MC,m_n]; 0037 0038 end 0039 0040 I = [Is_supl;MC]; 0041 I(gnd_ind,:) = 0; 0042 0043 % stupidity to be matlab 6+7 compatible 0044 if nargin < 7; v_f = zeros(size(E,1),size(I,2)); end 0045 if isempty(v_f);v_f = zeros(size(E,1),size(I,2)); end 0046 0047 maxiter=10000; % This should be high enough, but it may maybe this should 0048 % depend on the number of measurements? 0049 0050 0051 if isreal(E)==1 0052 0053 ver = eidors_obj('interpreter_version'); 0054 0055 if ver.isoctave % OCtave doesn't have Cholinc yet (as of 2.9.13) 0056 M= []; 0057 else 0058 opts.droptol = tol/10; 0059 if ver.ver < 7.012 0060 M = cholinc(E,opts.droptol); 0061 else 0062 opts.droptol = 1e-6; 0063 opts.type = 'ict'; %otherwise droptol is ignored opts.type = 'nofill'; 0064 0065 % M = ichol(E,opts); 0066 % ichol makes pcg even slower. It's better to use no pre-conditioner 0067 M = []; 0068 end 0069 end 0070 0071 for y=1:size(MC,2) 0072 %Set this line to suit your approximation needs. *************** 0073 %for more details use help pcg on Matlab's command window. 0074 [v_f(:,y),flag,relres,iter,resvec] = pcg(E,I(:,y), ... 0075 tol*norm(I(:,y)),maxiter,M',M,v_f(:,y)); 0076 end 0077 0078 else %is real 0079 0080 %Preconditioner 0081 % OCtave doesn't have Cholinc yet (as of 2.9.13) 0082 if exist('OCTAVE_VERSION') 0083 L= []; U=[]; 0084 else 0085 [L,U] = luinc(E,tol/10); 0086 end 0087 0088 for y=1:size(MC,2) 0089 0090 [v_f(:,y),flag,relres,iter,resvec] = bicgstab(E,I(:,y), ... 0091 tol*norm(I(:,y)),maxiter,L,U); 0092 0093 end 0094 end %is complex 0095 0096 0097 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0098 % This is part of the EIDORS suite. 0099 % Copyright (c) N. Polydorides 2003 0100 % Copying permitted under terms of GNU GPL 0101 % See enclosed file gpl.html for details. 0102 % EIDORS 3D version 2.0 0103 % MATLAB version 5.3 R11 0104 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%