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