calc_error_norms_for_square_domain

Get forward model of the img and the conductivity per element

0001 function [L2_tot_error,H1semi_tot_error,H1_tot_error,I_err,U_errS,U_errM,U_errSM,timing_solver,DOF]=error_2D_squ_CEM(img,eletype,plot_on)
0002 
0003 %Get forward model of the img and the conductivity per element
0004 img.fwd_model.approx_type=eletype;
0005 mdl=img.fwd_model;
0006 
0007 %Copy these images for the forward solution
0008 img2=img; mdl2=mdl;
0009 
0010 
0011 %Modify forward model to ensure elements are reordered and get the new
0012 %nodes and elements
0013 if ~isfield(mdl,'approx_type')    || ...
0014    strcmp(mdl.approx_type,'tri3') || ...
0015    strcmp(mdl.approx_type,'tet4')   
0016 else
0017 mdl.approx_type=eletype; [bound,elem,nodes]=fem_1st_to_higher_order(mdl); 
0018 mdl.boundary=bound; mdl.elems=elem; mdl.nodes=nodes;
0019 img.fwd_model=mdl;
0020 end
0021 
0022 %Calculate number of nodes and elecs
0023 nnodesF=length(mdl.nodes(:,1)); nelecsF=length(mdl.electrode);
0024 DOF = nnodesF + nelecsF;
0025 tic;
0026 %Calculate stiffness matrix
0027 A_mat = system_mat_higher_order(img);
0028 At=A_mat.E;
0029 ATL = At(1:nnodesF,1:nnodesF);
0030 ATR = At(1:nnodesF,nnodesF+1:nnodesF+nelecsF);
0031 ABL = At(nnodesF+1:nnodesF+nelecsF,1:nnodesF);
0032 ABR = At(nnodesF+1:nnodesF+nelecsF,nnodesF+1:nnodesF+nelecsF);
0033 
0034 %Form a newmatirx
0035 AtN=zeros(nnodesF+nelecsF,nnodesF+nelecsF);     
0036 AtN(1:nnodesF,1:nnodesF) = ATL;
0037 AtN(1:nnodesF,nnodesF+1:nnodesF+nelecsF) = 0;
0038 AtN(nnodesF+1:nnodesF+nelecsF,1:nnodesF)=ABL;
0039 AtN(nnodesF+1:nnodesF+nelecsF,nnodesF+1:nnodesF+nelecsF)=-eye(nelecsF);         
0040 
0041 %Form a RHS vector
0042 %We do +1 on E1 and -1 on E2
0043 Uvec=[1;-1]; 
0044 RHSvec(1:nnodesF,1) = -ATR*Uvec;
0045 RHSvec(nnodesF+1:nnodesF+nelecsF)=-ABR*Uvec;
0046 
0047 %Inverse
0048 uI = AtN \ RHSvec;
0049 volt_nodes=uI(1:nnodesF);
0050 timing_solver=toc;
0051 
0052 %Find elemental stiffness matrix on the REFINED MODEL
0053 [~,elemstiff]=mc_calc_stiffness2(mdl,img);
0054 
0055 
0056 %ANALYTIC Solution
0057 
0058 %Okay find the element matrix
0059 elemstruc=mdl.elems; nodestruc=mdl.nodes;
0060 
0061 %Find elemental stiffness matrix
0062 [~,elemstiff]=mc_calc_stiffness2(mdl,img);
0063 
0064 %Find no. of elements, and initialize vector of H1 errors
0065 nelems=size(elemstruc,1); H1_error=zeros(nelems,1);
0066 
0067 %Analytic solution by Fourier decomposition - Two electrodes only
0068 %We have infinite system of matrices U=SA
0069 %U is integral of potential, S is integral of cos products and A are coeffs
0070 
0071 %Impedance, COM and half width of each electrode
0072 z_c1=mdl.electrode(1).z_contact; x1=3*pi/10; w1=pi/10;
0073 z_c2=mdl.electrode(2).z_contact; x2=7*pi/10; w2=pi/10; 
0074 %We need ratio of Ul/zl
0075 uz1=Uvec(1)/z_c1; uz2=Uvec(2)/z_c2;
0076 
0077 %Number of terms to truncate the analytic and interior terms
0078 n_trunc=1000; n_int_trunc=225;
0079 
0080 %
0081 Sc=zeros(n_trunc+1,n_trunc+1);
0082 Uc=zeros(n_trunc+1,1); Ac=zeros(n_trunc+1,1); 
0083 
0084 %We first compute the sl(integer) term sl(
0085 sln1=zeros(6*n_trunc+1); sln2=zeros(6*n_trunc+1); 
0086 sln1(3*n_trunc+1)=2*w1; sln2(3*n_trunc+1)=2*w2;
0087 for i=1:3*n_trunc
0088    sln1(3*n_trunc+1+i) = ( sin(i*(x1+w1))  - sin(i*(x1-w1))  )/i; 
0089    sln1(3*n_trunc+1-i) = sln1(3*n_trunc+1+i);
0090    sln2(3*n_trunc+1+i) = ( sin(i*(x2+w2))  - sin(i*(x2-w2))  )/i; 
0091    sln2(3*n_trunc+1-i) = sln2(3*n_trunc+1+i);   
0092 end
0093 
0094 %Assemble matrices and then invert
0095 for i=1:n_trunc+1
0096    %First workout the t coefficnet
0097    Uc(i) = uz1*sln1(3*n_trunc+i) + uz2*sln2(3*n_trunc+i) ;   
0098    for j=1:n_trunc+1       
0099         Sc(i,j) = 0.5*(1/z_c1)*(sln1(j-i+3*n_trunc+1)+sln1(j+i-2+3*n_trunc+1))...
0100                 + 0.5*(1/z_c2)*(sln2(j-i+3*n_trunc+1)+sln2(j+i-2+3*n_trunc+1));
0101         if(i==j) %Diagonal term?
0102             Sc(i,j) = Sc(i,j) + tanh((i-1)*pi)*(i-1)*pi/2;
0103         end        
0104    end
0105 end
0106 
0107 %Solve and rescle
0108 Ac=Sc\Uc;
0109 
0110 %Loop over the nodes and find the boundary
0111 volt_analytic=zeros(length(img.fwd_model.nodes(:,1)),1);
0112 for in=1:length(img.fwd_model.nodes(:,1));
0113    xin=img.fwd_model.nodes(in,1);
0114    yin=img.fwd_model.nodes(in,2);
0115    volt_analytic(in)=Ac(1); %constant
0116    if(yin==0)
0117        for k=1:n_trunc
0118            volt_analytic(in) = volt_analytic(in) + ...
0119                (Ac(k+1))*cos(k*xin);
0120        end
0121    else
0122        for k=1:n_int_trunc
0123            volt_analytic(in) = volt_analytic(in) + ...
0124                (Ac(k+1))*cos(k*xin)*cosh(k*(pi-yin))/cosh(k*pi);
0125        end
0126    end
0127 end
0128 
0129 %Determine the currents
0130 %Current 1 and 2 constant part
0131 I_analytic(1) = 2*w1/z_c1*(Uvec(1) - Ac(1)); 
0132 I_analytic(2) = 2*w2/z_c2*(Uvec(2) - Ac(1));
0133 for kkk=1:n_trunc
0134    I_analytic(1) = I_analytic(1) - Ac(kkk+1)/(kkk*z_c1)*( sin(kkk*(x1+w1)) - sin(kkk*(x1-w1)) );
0135    I_analytic(2) = I_analytic(2) - Ac(kkk+1)/(kkk*z_c2)*( sin(kkk*(x2+w2)) - sin(kkk*(x2-w2)) );   
0136 end
0137 I_FEM(1) = uI(nnodesF+1);
0138 I_FEM(2) = uI(nnodesF+2);
0139 
0140 %Calculate the 2-norm error
0141 I_err = norm(I_FEM-I_analytic,2);
0142 
0143 %Now solve with curremts
0144 %Ok now lets's solve with the analytic currents
0145 IRHS = zeros(nnodesF+nelecsF,1);
0146 IRHS(nnodesF+1)=I_analytic(1);
0147 IRHS(nnodesF+2)=I_analytic(2);
0148 
0149 %Create index vector and eliminate ground node equation from index
0150 groundnode=mdl.gnd_node; idx=1:size(At,1); idx(groundnode)=[];
0151 
0152 %Ok now create empty potentials
0153 uU = zeros(nnodesF+nelecsF,1);
0154 uU(idx) = left_divide(At(idx,idx),IRHS(idx,:));
0155 %Subtract zero of potential - x=pi/2 i.e. from symmetry applied potential
0156 uU = uU-0.5*(uU(13)+uU(19));
0157 
0158 %Now the potential on stim electrodes
0159 Uvec_FEM(1) = uU(nnodesF+1);
0160 Uvec_FEM(2) = uU(nnodesF+2);
0161 
0162 %Error just on stim electrodes
0163 U_errS=norm(Uvec-Uvec_FEM',2);
0164 
0165 %Error just on meas (point) electrodes
0166 bound_nodes_not_elecs=img.fwd_model.bound_nodes_not_elecs;
0167 UvecM=volt_analytic(bound_nodes_not_elecs);
0168 UvecM_FEM=uU(bound_nodes_not_elecs);
0169 U_errM = norm(UvecM-UvecM_FEM);
0170 
0171 %Error on stim and meas(point) electrodes
0172 UvecSM=[Uvec',UvecM'];
0173 UvecSM_FEM=[Uvec_FEM,UvecM_FEM'];
0174 U_errSM=norm(UvecSM-UvecSM_FEM);
0175 
0176 %Plot the solution
0177 if(plot_on==1 && (z_c1==0.00001 || z_c2==1000) && strcmp(mdl.approx_type,'tri10'))
0178 figure; plot3(img.fwd_model.nodes(:,1),img.fwd_model.nodes(:,2),volt_nodes,'r*')
0179 hold on; plot3(img.fwd_model.nodes(:,1),img.fwd_model.nodes(:,2),volt_analytic,'b*');
0180 title('Nodal analytic and FEM solutions to CEM');
0181 legend('FEM solution','Analytic solution');
0182 xlabel('x'); ylabel('y'); zlabel('u(x,y)')
0183 
0184 if(z_c1==0.00001)
0185     print_convert error_rates_contact_impedance02a.png;
0186 elseif(z_c1==1000)
0187     print_convert error_rates_contact_impedance02b.png;
0188 end
0189 end
0190 
0191 %Get the elements we want
0192 v=1:nelems;
0193 
0194 %Get the basis and gradients at the knot points
0195 eletype=mdl.approx_type; 
0196 if(strcmp(eletype,'tri3'))
0197     dim=2; order1=0; order2=2;
0198 elseif(strcmp(eletype,'tri6'))    
0199     dim=2; order1=2; order2=4;
0200 elseif(strcmp(eletype,'tri10'))
0201     dim=2; order1=4; order2=7;
0202 elseif(strcmp(eletype,'tet4'))
0203     dim=3; order1=0; order2=2;
0204 elseif(strcmp(eletype,'tet10'))
0205     dim=3; order1=2; order2=4;
0206 else  
0207     error('Element type not recognised for integration rules');
0208 end
0209 %Find shape function gradient at knot points
0210 [weight1,xcoord1,ycoord1,zcoord1]=gauss_points(dim,order1);
0211 for kk=1:size(weight1,2)
0212     dphi_sol(:,:,kk) = element_d_shape_function(eletype,xcoord1(kk),ycoord1(kk),zcoord1(kk));
0213 end
0214 %Find shape functions at knot points (these are higher order ones now)
0215 [weight2,xcoord2,ycoord2,zcoord2]=gauss_points(dim,order2);
0216 for kk=1:size(weight2,2)
0217     phi_sol(:,kk) = element_shape_function(eletype,xcoord2(kk),ycoord2(kk),zcoord2(kk));
0218 end
0219 %Initialsie the erros
0220 H1semi_error=zeros(nelems,1);
0221 L2_error=zeros(nelems,1);
0222 
0223 %Loop through elems and calculate the error
0224 for j=1:length(v)  
0225     %Get the element
0226     i=v(j);    
0227     %Get the node numbers
0228     eleminodelist=elemstruc(i,:);
0229     
0230     %List by row of coordinate on the element
0231     thise = nodestruc(eleminodelist,:);
0232     
0233     %Find the Jacobian of the mapping in 2D and 3D
0234     jacobianelem = [thise(2,1)-thise(1,1),thise(2,2)-thise(1,2); ...
0235             thise(3,1)-thise(1,1),thise(3,2)-thise(1,2)];  
0236     
0237     %Find the magnitude of the Jacobian of the mapping
0238     % magjacelem=det(jacobianelem);
0239     magjacelem=abs(det(jacobianelem));
0240     
0241     %L2 error of solution
0242     %Loop over the weights and evaluate the sensitivty contribution
0243     for ii=1:size(weight2,2)
0244         %Find the isoparametric map basis functions at the local point
0245         map=element_shape_function('tri3',xcoord2(ii),ycoord2(ii),zcoord2(ii));
0246         
0247         %Map our local coordinate to global coordinates using map
0248         cart(1)=thise(1,1)*map(1)+thise(2,1)*map(2)+thise(3,1)*map(3);
0249         cart(2)=thise(1,2)*map(1)+thise(2,2)*map(2)+thise(3,2)*map(3);      
0250         
0251         %Find the analytic solution at this point
0252         analytic_sol = Ac(1);        
0253         if(cart(2)==0)
0254            for k=1:n_trunc
0255            analytic_sol = analytic_sol + (Ac(k+1))*cos(k*cart(1));
0256            end
0257         else
0258            for k=1:n_int_trunc
0259            analytic_sol = analytic_sol + ...
0260                (Ac(k+1))*cos(k*cart(1))*cosh(k*(pi-cart(2)))/cosh(k*pi);
0261            end
0262         end
0263         
0264         %Find the fem solution
0265         elemi_volt_nodes=volt_nodes(eleminodelist)'; %Vector
0266         fem_sol=elemi_volt_nodes*phi_sol(:,ii);
0267         
0268         %Difference in solution square
0269         diff_sol=(fem_sol-analytic_sol)^2;
0270                                 
0271         %Compute the difference and multiply by weight
0272         diff_sol=diff_sol*weight2(ii);
0273         
0274         %Add the contribution to the elemental sensitivity
0275         L2_error(i)=L2_error(i)+diff_sol;        
0276     end
0277     
0278     
0279     %H1 error of solution
0280     %Loop over the weights and evaluate the sensitivty contribution
0281     for ii=1:size(weight1,2)
0282         %Find the isoparametric map basis functions at the local point
0283         map=element_shape_function('tri3',xcoord1(ii),ycoord1(ii),zcoord1(ii));
0284         
0285         %Map our local coordinate to global coordinates using map
0286         cart(1)=thise(1,1)*map(1)+thise(2,1)*map(2)+thise(3,1)*map(3);
0287         cart(2)=thise(1,2)*map(1)+thise(2,2)*map(2)+thise(3,2)*map(3);      
0288                       
0289         %Find the analytic solution at this point
0290         analytic_sol(1) = 0; analytic_sol(2)=0;        
0291         if(cart(2)==0)
0292            for k=1:n_trunc
0293            analytic_sol(1) = analytic_sol(1) - ...
0294                k*(Ac(k+1))*sin(k*cart(1));
0295            analytic_sol(2) = analytic_sol(2) - ...
0296                k*(Ac(k+1))*cos(k*cart(1))*tanh(k*pi);           
0297            end
0298         else
0299            for k=1:n_int_trunc
0300            analytic_sol(1) = analytic_sol(1) - ...
0301                k*(Ac(k+1))*sin(k*cart(1))*cosh(k*(pi-cart(2)))/cosh(k*pi);
0302            analytic_sol(2) = analytic_sol(2) - ...
0303                k*(Ac(k+1))*cos(k*cart(1))*sinh(k*(pi-cart(2)))/cosh(k*pi);  
0304            end
0305         end        
0306         
0307         
0308         %Find the fem solution gradient
0309         elemi_volt_nodes=volt_nodes(eleminodelist)'; %Vector of local sols
0310         fem_sol=elemi_volt_nodes*(jacobianelem\dphi_sol(:,:,ii))'; %Gradient
0311         
0312         %Difference in solution square
0313         diff_sol=(fem_sol-analytic_sol)*(fem_sol-analytic_sol)';
0314                                 
0315         %Compute the difference and multiply by weight
0316         diff_sol=diff_sol*weight1(ii);
0317         
0318         %Add the contribution to the elemental sensitivity
0319         H1semi_error(i)=H1semi_error(i)+diff_sol;        
0320     end    
0321     
0322     %We have the sensitivity on reference and multiply by Jacobian
0323     H1semi_error(i)=H1semi_error(i)*magjacelem;   
0324     L2_error(i)=L2_error(i)*magjacelem; 
0325             
0326     %Total error is the sum
0327     H1_error(i)=H1semi_error(i)+L2_error(i);
0328 end
0329 
0330 
0331 %Now find the total H1_error and square root to get norm
0332 H1_tot_error=sum(H1_error); H1_tot_error=sqrt(H1_tot_error);
0333 H1semi_tot_error=sum(H1semi_error); H1semi_tot_error=sqrt(H1semi_tot_error);
0334 L2_tot_error=sum(L2_error); L2_tot_error=sqrt(L2_tot_error);
0335 
0336 
0337 end