calc_mass_stiffness

STIFFNESS MATRIX PART
Stiffness matrix, including piecewise conductivity, for EIT. The second
argument is for Jacobian, it gives discretied gradients in element.

0001 %STIFFNESS MATRIX PART
0002 function [Agal,elemstiff]=mc_calc_stiffness2(fwd_model,img)
0003 %Stiffness matrix, including piecewise conductivity, for EIT. The second
0004 %argument is for Jacobian, it gives discretied gradients in element.
0005 
0006 %If function called only with image, extract forward model
0007 if(nargin==1)
0008     img=fwd_model; fwd_model=img.fwd_model;
0009 end
0010 
0011 %Cache node structure and find no. of spatial dimensions and nodes
0012 %Cache element structure and find no. of elements
0013 nodestruc=fwd_model.nodes; nodedim=size(nodestruc,2); nnodes=size(nodestruc,1); 
0014 elemstruc=fwd_model.elems; nelems=size(elemstruc,1);
0015 
0016 %Find quadrature points/weights for integration by switching between cases
0017 eletype=fwd_model.approx_type; 
0018 if(strcmp(eletype,'tri3'))
0019     dim=2; order1=0; order2=2;
0020 elseif(strcmp(eletype,'tri6'))    
0021     dim=2; order1=2; order2=4;
0022 elseif(strcmp(eletype,'tri10'))
0023     dim=2; order1=4; order2=7;
0024 elseif(strcmp(eletype,'tet4'))
0025     dim=3; order1=0; order2=2;
0026 elseif(strcmp(eletype,'tet10'))
0027     dim=3; order1=2; order2=4;
0028 else  
0029     error('Element type not recognised for integration rules');
0030 end
0031 [weight1,xcoord1,ycoord1,zcoord1]=gauss_points(dim,order1);
0032 for kk=1:size(weight1,2)
0033     dphi(:,:,kk) = element_d_shape_function(eletype,xcoord1(kk),ycoord1(kk),zcoord1(kk));
0034 end
0035 
0036 [weight2,xcoord2,ycoord2,zcoord2]=gauss_points(dim,order2);
0037 for kk=1:size(weight2,2)
0038     phi(:,kk) = element_shape_function(eletype,xcoord2(kk),ycoord2(kk),zcoord2(kk));
0039 end
0040 
0041 %Initialise global stiffness matrix
0042 Agal=zeros(nnodes,nnodes); %sparse updating non zero slow
0043 
0044 %Loop over the elements and calculate local Am matrix
0045 for i=nelems:-1:1
0046     %Find the list of node numbers for each element
0047     eleminodelist=elemstruc(i,:);
0048     
0049     %List by row of coordinate on the element
0050     thise = nodestruc(eleminodelist,:);
0051     
0052     %Find the Jacobian of the mapping in 2D and 3D
0053     if(nodedim==2); jacobianelem = ... %2D Jacobian of mapping
0054             [thise(2,1)-thise(1,1),thise(2,2)-thise(1,2); ...
0055             thise(3,1)-thise(1,1),thise(3,2)-thise(1,2)];  
0056     elseif(nodedim==3); jacobianelem = ... %3D Jacobian of mapping
0057             [thise(2,1)-thise(1,1),thise(2,2)-thise(1,2),thise(2,3)-thise(1,3); ...
0058             thise(3,1)-thise(1,1),thise(3,2)-thise(1,2),thise(3,3)-thise(1,3); ...
0059             thise(4,1)-thise(1,1),thise(4,2)-thise(1,2),thise(4,3)-thise(1,3)];
0060     end
0061     
0062     %Find the magnitude of the Jacobian of the mapping
0063     % magjacelem=det(jacobianelem);
0064     magjacelem=abs(det(jacobianelem));
0065            
0066     %Elemental stiffness matrices
0067     Ammat=0;
0068     for kk=1:size(weight1,2)
0069         Ammat = Ammat + weight1(kk)* ...
0070             (jacobianelem\dphi(:,:,kk))'* ...
0071             (jacobianelem\dphi(:,:,kk))*magjacelem;
0072     end
0073 
0074     %Element mass matrices
0075     Mmmat=0;
0076     for kk=1:size(weight2,2)
0077         Mmmat = Mmmat + weight2(kk)* ...
0078             (phi(:,kk))* ...
0079             (phi(:,kk))' * magjacelem;
0080     end
0081     
0082     %Store the Ammat without multiplication of conductivity for Jacobian
0083     elemstiff(i).elemstiff = Ammat;
0084     elemstiff(i).elemmass  = Mmmat;
0085    
0086     %This is element stiffness matrix (and multiply by its conductivity)
0087     stiff=Ammat*img.elem_data(i); 
0088     
0089     %Assemble global stiffness matrix (Silvester's book!!)
0090     Agal(elemstruc(i,:), elemstruc(i,:)) = Agal(elemstruc(i,:), elemstruc(i,:)) + stiff;
0091 
0092 end
0093 
0094 end