system_mat_higher_order

Assemble the total stiffness matrix : s_mat.E=At;
M Crabb - 29.06.2012
TODO - Sparse assignment of the matrices

0001 function [s_mat]=system_mat_higher_order(fwd_model,img)
0002 %Assemble the total stiffness matrix : s_mat.E=At;
0003 %M Crabb - 29.06.2012
0004 %TODO - Sparse assignment of the matrices
0005 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0006 
0007 if nargin == 1
0008    img= fwd_model;
0009 elseif  strcmp(getfield(warning('query','EIDORS:DeprecatedInterface'),'state'),'on')
0010    warning('EIDORS:DeprecatedInterface', ...
0011       ['Calling SYSTEM_MAT_HIGHER_ORDER with two arguments is deprecated and will cause' ...
0012        ' an error in a future version. First argument ignored.']);
0013 end
0014 fwd_model= img.fwd_model;
0015 
0016 %Find no. of electrodes and no. of ndoes
0017 elecstruc=fwd_model.electrode; nelecs=size(elecstruc,2);
0018 nodestruc=fwd_model.nodes; nnodes=size(nodestruc,1); 
0019 
0020 %Test - Point/Complete electrodes. Assume no mixed model so test first elect
0021 if isfield(elecstruc,'faces');
0022    error('Can''t process faces-type electrodes. See mat_idx_to_electrode for help');
0023 end
0024 if(size(elecstruc(1).nodes,2)==1 && size(elecstruc(1).nodes,1)==1) %POINT ELECTRODE
0025     %IF POINT ELECTRODE
0026     [At,elemstiff]=mc_calc_stiffness(fwd_model,img);
0027 else %COMPLETE ELECTRODE
0028     [Am,elemstiff]=mc_calc_stiffness(fwd_model,img);
0029     
0030      [Aw,Az,Ad]=mc_calc_complete(fwd_model);
0031      At=sparse(nnodes+nelecs,nnodes+nelecs);     
0032 %     [i,j,s] = find(Am);
0033 %     At=At+sparse(i,j,s,nnodes+nelecs,nnodes+nelecs);
0034 %     [i,j,s] = find(Az);
0035 %     At=At+sparse(i,j,s,nnodes+nelecs,nnodes+nelecs);
0036 %     [i,j,s] = find(Aw);
0037 %     At=At+sparse(i,j+nnodes,s,nnodes+nelecs,nnodes+nelecs);
0038 %     At=At+sparse(j+nnodes,i,s,nnodes+nelecs,nnodes+nelecs);
0039 %     [i,j,s] = find(Ad);
0040 %     At=At+sparse(i+nnodes,j+nnodes,s,nnodes+nelecs,nnodes+nelecs);
0041     At(1:nnodes,1:nnodes) = Am+Az;
0042     At(1:nnodes,nnodes+1:nnodes+nelecs) = Aw;
0043     At(nnodes+1:nnodes+nelecs,1:nnodes)=Aw';
0044     At(nnodes+1:nnodes+nelecs,nnodes+1:nnodes+nelecs)=Ad;
0045 end
0046 
0047 %Put in structure to be compatibile with eidors
0048 s_mat.E=sparse(At);
0049 
0050 %Store individual stiffness matrices for Jacobian
0051 s_mat.elemstiff=elemstiff;
0052 end
0053 
0054 %STIFFNESS MATRIX PART
0055 function [Agal,elemstiff]=mc_calc_stiffness(fwd_model,img)
0056 %Stiffness matrix, including piecewise conductivity, for EIT. The second
0057 %argument is for Jacobian, it gives discretied gradients in element.
0058 
0059 %If function called only with image, extract forward model
0060 if(nargin==1)
0061     img=fwd_model; fwd_model=img.fwd_model;
0062 end
0063 
0064 %Cache node structure and find no. of spatial dimensions and nodes
0065 %Cache element structure and find no. of elements
0066 nodestruc=fwd_model.nodes; nodedim=size(nodestruc,2); nnodes=size(nodestruc,1); 
0067 elemstruc=fwd_model.elems; nelems=size(elemstruc,1);
0068 
0069 %Find quadrature points/weights for integration by switching between cases
0070 eletype=fwd_model.approx_type; 
0071 if(strcmp(eletype,'tri3'))
0072     dim=2; order=0;
0073 elseif(strcmp(eletype,'tri6'))
0074     dim=2; order=2;
0075 elseif(strcmp(eletype,'tri10'))
0076     dim=2; order=4;
0077 elseif(strcmp(eletype,'tet4'))
0078     dim=3; order=0;
0079 elseif(strcmp(eletype,'tet10'))
0080     dim=3; order=2;
0081 else  
0082     error('Element type not recognised for integration rules');
0083 end
0084 [weight,xcoord,ycoord,zcoord]=gauss_points(dim,order);
0085 
0086 %Find derivative of shape function on domain element
0087 for kk=size(weight,2):-1:1
0088     dphi(:,:,kk) = element_d_shape_function(eletype,xcoord(kk),ycoord(kk),zcoord(kk));
0089 end
0090 
0091 %Initialise global stiffness matrix
0092 %Agal=sparse(nnodes,nnodes);
0093 
0094 %Loop over the elements and calculate local Am matrix
0095 for i=nelems:-1:1
0096     %Find the list of node numbers for each element
0097     eleminodelist=elemstruc(i,:);
0098     
0099     %List by row of coordinate on the element
0100     thise = nodestruc(eleminodelist,:);
0101     
0102     %Find the Jacobian of the mapping in 2D and 3D
0103     if(nodedim==2); jacobianelem = ... %2D Jacobian of mapping
0104             [thise(2,1)-thise(1,1),thise(2,2)-thise(1,2); ...
0105             thise(3,1)-thise(1,1),thise(3,2)-thise(1,2)];  
0106     elseif(nodedim==3); jacobianelem = ... %3D Jacobian of mapping
0107             [thise(2,1)-thise(1,1),thise(2,2)-thise(1,2),thise(2,3)-thise(1,3); ...
0108             thise(3,1)-thise(1,1),thise(3,2)-thise(1,2),thise(3,3)-thise(1,3); ...
0109             thise(4,1)-thise(1,1),thise(4,2)-thise(1,2),thise(4,3)-thise(1,3)];
0110     end
0111     
0112     %Find the magnitude of the Jacobian of the mapping
0113     % magjacelem=det(jacobianelem);
0114     magjacelem=abs(det(jacobianelem));
0115 
0116     %Initialise and find elemental stiffness matrices
0117     Ammat=0;
0118     for kk=1:size(weight,2)
0119 %         Ammat = Ammat + weight(kk)* ...
0120 %             (jacobianelem\dphi(:,:,kk))'* ...
0121 %             (jacobianelem\dphi(:,:,kk))*magjacelem;
0122         jdphitemp=jacobianelem\dphi(:,:,kk);
0123         Ammat = Ammat + weight(kk)* ...
0124             (jdphitemp'*jdphitemp)*magjacelem;
0125     end
0126 
0127     %SPEED UP
0128     %Can we get system_mat_fields here to speed Jacobian?
0129     
0130     %Store the Ammat without multiplication of conductivity for Jacobian
0131     elemstiff(i).elemstiff=Ammat;
0132    
0133     %This is element stiffness matrix (and multiply by its conductivity)
0134 
0135     stiff=Ammat*img.elem_data(i); 
0136     
0137     %Assemble global stiffness matrix (Silvester's book!!)
0138 %   Agal(elemstruc(i,:), elemstruc(i,:)) = Agal(elemstruc(i,:), elemstruc(i,:)) + stiff;
0139 %   [ii,jj]= meshgrid(eleminodelist); ll = length(eleminodelist)^2; %% meshgrid CRAZY SLOW
0140     ll = length(eleminodelist);
0141     ii = eleminodelist(ones(ll,1),:); jj=ii';
0142     idx = (i-1)*ll^2 + (1:ll^2);
0143     AI(idx) = ii(:);
0144     AJ(idx) = jj(:);
0145     Aval(idx) = stiff(:); 
0146 
0147 end
0148 %Agal2= sparse(AI, AJ, Aval, nnodes, nnodes); norm(Agal-Agal2,'fro')
0149 Agal= sparse(AI, AJ, Aval, nnodes, nnodes);
0150 end
0151 
0152 %COMPLETE ELECTRODE MATRICES
0153 function [Aw,Az,Ad]=mc_calc_complete(fwd_model)
0154 %Takes a forward model and calculates Az, Aw, Ad for complete electrode
0155 
0156 %Get the electrode structure, find number of electrodes
0157 %Get the boundary strucutre, find number of boundaries
0158 %Get the node structrue, find number of nodes and problem dim
0159 elecstruc=fwd_model.electrode; nelecs=size(elecstruc,2);boundstruc=fwd_model.boundary; nodestruc=fwd_model.nodes; 
0160 nnodes=size(nodestruc,1); nodedim=size(nodestruc,2);
0161 
0162 %Connect boundary/electrode -Put boundary into old matrix strucutre
0163 %for i=nbounds:-1:1
0164 %    boundstrucold(i,:)=boundstruc(i).nodes;
0165 %end
0166 
0167 %Find quadrature points/weights for integration by switching between cases
0168 eletype=fwd_model.approx_type; 
0169 if(strcmp(eletype,'tri3'))
0170     dim=1; order=2;
0171 elseif(strcmp(eletype,'tri6'))
0172     dim=1; order=4;
0173 elseif(strcmp(eletype,'tri10'))
0174     dim=1; order=6;
0175 elseif(strcmp(eletype,'tet4'))
0176     dim=2; order=2;
0177 elseif(strcmp(eletype,'tet10'))
0178     dim=2; order=4;
0179 else  
0180     error('Element type not recognised for integration rules');
0181 end
0182 [weight,xcoord,ycoord]=gauss_points(dim,order);
0183 
0184 %Find shape function on boundary element
0185 for kk=size(weight,2):-1:1
0186     phi(:,kk) = boundary_shape_function(eletype,xcoord(kk),ycoord(kk))';
0187 end
0188 
0189 %1. Initialise global Az/Aw/Ad matrices and assemble a la Silvester
0190 Az=sparse(nnodes,nnodes); %sparse updating non zero slow
0191 Aw=zeros(nnodes,nelecs); Ad=zeros(nelecs,nelecs);
0192 
0193 %Loop over the electrodes
0194 for ke=1:nelecs    
0195     %The boundary numbers and areas, outputs rows of mdl.boundary of electrode
0196     [bdy_idx,bdy_area]=find_electrode_bdy(boundstruc(:,1:nodedim),nodestruc,elecstruc(ke).nodes);
0197     
0198     %Store boundary numbers, and corresponding areas
0199     boundidx_ke=bdy_idx; area_ke=bdy_area;
0200     
0201     %Find contact impedance of electrode
0202     elecimped=elecstruc(ke).z_contact;   
0203            
0204     %Find total electrode area (absolute values)
0205     elecarea=0;
0206     for i=1:size(area_ke,2)
0207         elecarea = elecarea + abs(area_ke(i));
0208     end
0209     
0210     %Form the matrix Ad
0211     Ad(ke,ke)=elecarea/elecimped; 
0212     
0213     
0214     %Loop over boundarys and calculate Aw/Az matrices
0215     for ii=1:length(boundidx_ke)
0216         %List by row of coordinates of on the boundaryNodal coordinates on the boundary
0217         thisb=nodestruc(boundstruc(boundidx_ke(ii),:),:);
0218     
0219         %Find the magnitude Jacobian of the mapping in 2D/3D
0220         %NB:Scalings are consistent with reference element shape
0221         if(nodedim==2)
0222             %Jacobian = 0.5*|(x2-x1)| (x1,x2 vector of coords)
0223             diff21=thisb(2,:)-thisb(1,:);
0224             magjacbound=0.5*sqrt(diff21(1)^2+diff21(2)^2);
0225         elseif(nodedim==3)
0226             %Jacobian = |(x3-x1)x(x3-x2)| (x1,x2,x3 vector of coords)
0227             diffprod=cross(thisb(3,:)-thisb(1,:),thisb(3,:)-thisb(2,:));
0228             magjacbound=sqrt(diffprod(1)^2+diffprod(2)^2+diffprod(3)^2);
0229         end
0230 
0231         %Initialise Azlocmat/Awlocmat and find local matrices
0232         Azmat=0; Awmat=0;
0233         for kk=1:size(weight,2)            
0234             temphikk = phi(:,kk);
0235             Azmat = Azmat + weight(kk)* ...
0236                 (temphikk*temphikk')*magjacbound;
0237             Awmat = Awmat + weight(kk)* ...
0238                 (phi(:,kk))'*magjacbound;            
0239             %Azmat = Azmat + weight(kk)* ...
0240             %    (phi(:,kk))* ...
0241             %    (phi(:,kk))'*magjacbound;
0242             %Awmat = Awmat + weight(kk)* ...
0243             %    (phi(:,kk))'*magjacbound;
0244         end         
0245         
0246         %Node numbers for this boundary
0247         boundnodes=boundstruc(boundidx_ke(ii),:);
0248         
0249         %Assemble the matrices
0250         Az(boundnodes,boundnodes) = Az(boundnodes,boundnodes)+Azmat/elecimped;
0251         Aw(boundnodes,ke) = Aw(boundnodes,ke) - Awmat'/elecimped;
0252     end
0253        
0254 end
0255 
0256 end
0257 
0258 function do_unit_test
0259    for i=0:100 ; switch i
0260       case 0;img = mk_image( mk_common_model('c2C2',16),1);
0261              str = 'c2C2';
0262              tst= [4,[1,1,1]/3];
0263              ap0 = 'tri3'; ap1 = 'tri6';
0264       case 1;img = mk_image( mk_common_model('b3cr',16),1);
0265              str = 'b3cr';
0266              tst= [0.268642857142857,0.028,0.009333333333333,0.028];
0267              ap0 = 'tet4'; ap1 = 'tet10';
0268       case 2;img = mk_image( ng_mk_cyl_models(1,[6,0.5],.1),1);
0269              str = 'ng_mk_cyl1';
0270              tst= [0.131643427733397,0.000323024820911,0,0];
0271              ap0 = 'tet4'; ap1 = 'tet10';
0272       case 3;img = mk_image( mk_common_model('d3cr',16),1);
0273              str = 'd3cr';
0274              tst= [0.11125, 0.01, 0.003333333333333,0.01];
0275              ap0 = 'tet4'; ap1 = 'tet10';
0276       case 4;img = mk_image( ng_mk_cyl_models([1,1,.1],[6,0.5],.1),1);
0277              str = 'ng_mk_cyl2';
0278              tst= [0.061014587493642, 0.000551387020511, 0, 0];
0279              ap0 = 'tet4'; ap1 = 'tet10';
0280       case 5;break
0281       end
0282 
0283       img.fwd_model.approx_type = ap0;
0284       S1= system_mat_1st_order( img );
0285       Sl= system_mat_higher_order( img );
0286       unit_test_cmp(['Linear ',str], S1.E, Sl.E, 1e-13) 
0287 
0288       img.fwd_model.approx_type = ap1;
0289       img.fwd_model.system_mat = @system_mat_higher_order;
0290       [img.fwd_model.boundary,img.fwd_model.elems,img.fwd_model.nodes] = ...
0291             fem_1st_to_higher_order(img.fwd_model);
0292       Sh= calc_system_mat( img );
0293       unit_test_cmp(['Linear ',str], Sh.E(1,1:4), tst, 1e-13) 
0294    end
0295 end