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.m 6502 2022-12-30 14:29:12Z aadler $ 0020 0021 0022 warning('EIDORS:deprecated','M_3D_FIELDS is deprecated as of 06-Jun-2012. '); 0023 0024 [vr,vc] = size(vtx); 0025 0026 Is_supl = zeros(vr,size(m_ind,1)); 0027 %no of electrodes x no of measurements (now currents)! 0028 0029 MC = []; 0030 0031 for i=1:size(m_ind,1) 0032 0033 m_n = zeros(el_no,1); 0034 0035 m_n(m_ind(i,1)) = 1; 0036 m_n(m_ind(i,2)) = -1; 0037 0038 MC = [MC,m_n]; 0039 0040 end 0041 0042 I = [Is_supl;MC]; 0043 I(gnd_ind,:) = 0; 0044 0045 % stupidity to be matlab 6+7 compatible 0046 if nargin < 7; v_f = zeros(size(E,1),size(I,2)); end 0047 if isempty(v_f);v_f = zeros(size(E,1),size(I,2)); end 0048 0049 maxiter=10000; % This should be high enough, but it may maybe this should 0050 % depend on the number of measurements? 0051 0052 0053 if isreal(E)==1 0054 0055 ver = eidors_obj('interpreter_version'); 0056 0057 if ver.isoctave % OCtave doesn't have Cholinc yet (as of 2.9.13) 0058 M= []; 0059 else 0060 opts.droptol = tol/10; 0061 if ver.ver < 7.012 0062 M = cholinc(E,opts.droptol); 0063 else 0064 opts.droptol = 1e-6; 0065 opts.type = 'ict'; %otherwise droptol is ignored opts.type = 'nofill'; 0066 0067 % M = ichol(E,opts); 0068 % ichol makes pcg even slower. It's better to use no pre-conditioner 0069 M = []; 0070 end 0071 end 0072 0073 for y=1:size(MC,2) 0074 %Set this line to suit your approximation needs. *************** 0075 %for more details use help pcg on Matlab's command window. 0076 [v_f(:,y),flag,relres,iter,resvec] = pcg(E,I(:,y), ... 0077 tol*norm(I(:,y)),maxiter,M',M,v_f(:,y)); 0078 end 0079 0080 else %is real 0081 0082 %Preconditioner 0083 % LUinc no longer available (matlab 6 or so) 0084 % [L,U] = luinc(E,tol/10); 0085 % [L,U] = ilu(E, struct('droptol',E/10)); 0086 opts.droptol = tol/10 0087 opts.type = 'nofill'; %otherwise droptol is ignored opts.type = 'nofill'; 0088 [L,U] = ilu(E, opts); 0089 0090 for y=1:size(MC,2) 0091 0092 [v_f(:,y),flag,relres,iter,resvec] = bicgstab(E,I(:,y), ... 0093 tol*norm(I(:,y)),maxiter,L,U); 0094 0095 end 0096 end %is complex 0097 0098 0099 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0100 % This is part of the EIDORS suite. 0101 % Copyright (c) N. Polydorides 2003 0102 % Copying permitted under terms of GNU GPL 0103 % See enclosed file gpl.html for details. 0104 % EIDORS 3D version 2.0 0105 % MATLAB version 5.3 R11 0106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%