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
0004 img.fwd_model.approx_type=eletype;
0005 mdl=img.fwd_model;
0006
0007
0008 img2=img; mdl2=mdl;
0009
0010
0011
0012
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
0023 nnodesF=length(mdl.nodes(:,1)); nelecsF=length(mdl.electrode);
0024 DOF = nnodesF + nelecsF;
0025 tic;
0026
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
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
0042
0043 Uvec=[1;-1];
0044 RHSvec(1:nnodesF,1) = -ATR*Uvec;
0045 RHSvec(nnodesF+1:nnodesF+nelecsF)=-ABR*Uvec;
0046
0047
0048 uI = AtN \ RHSvec;
0049 volt_nodes=uI(1:nnodesF);
0050 timing_solver=toc;
0051
0052
0053 [~,elemstiff]=mc_calc_stiffness2(mdl,img);
0054
0055
0056
0057
0058
0059 elemstruc=mdl.elems; nodestruc=mdl.nodes;
0060
0061
0062 [~,elemstiff]=mc_calc_stiffness2(mdl,img);
0063
0064
0065 nelems=size(elemstruc,1); H1_error=zeros(nelems,1);
0066
0067
0068
0069
0070
0071
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
0075 uz1=Uvec(1)/z_c1; uz2=Uvec(2)/z_c2;
0076
0077
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
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
0095 for i=1:n_trunc+1
0096
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)
0102 Sc(i,j) = Sc(i,j) + tanh((i-1)*pi)*(i-1)*pi/2;
0103 end
0104 end
0105 end
0106
0107
0108 Ac=Sc\Uc;
0109
0110
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);
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
0130
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
0141 I_err = norm(I_FEM-I_analytic,2);
0142
0143
0144
0145 IRHS = zeros(nnodesF+nelecsF,1);
0146 IRHS(nnodesF+1)=I_analytic(1);
0147 IRHS(nnodesF+2)=I_analytic(2);
0148
0149
0150 groundnode=mdl.gnd_node; idx=1:size(At,1); idx(groundnode)=[];
0151
0152
0153 uU = zeros(nnodesF+nelecsF,1);
0154 uU(idx) = left_divide(At(idx,idx),IRHS(idx,:));
0155
0156 uU = uU-0.5*(uU(13)+uU(19));
0157
0158
0159 Uvec_FEM(1) = uU(nnodesF+1);
0160 Uvec_FEM(2) = uU(nnodesF+2);
0161
0162
0163 U_errS=norm(Uvec-Uvec_FEM',2);
0164
0165
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
0172 UvecSM=[Uvec',UvecM'];
0173 UvecSM_FEM=[Uvec_FEM,UvecM_FEM'];
0174 U_errSM=norm(UvecSM-UvecSM_FEM);
0175
0176
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
0192 v=1:nelems;
0193
0194
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
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
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
0220 H1semi_error=zeros(nelems,1);
0221 L2_error=zeros(nelems,1);
0222
0223
0224 for j=1:length(v)
0225
0226 i=v(j);
0227
0228 eleminodelist=elemstruc(i,:);
0229
0230
0231 thise = nodestruc(eleminodelist,:);
0232
0233
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
0238
0239 magjacelem=abs(det(jacobianelem));
0240
0241
0242
0243 for ii=1:size(weight2,2)
0244
0245 map=element_shape_function('tri3',xcoord2(ii),ycoord2(ii),zcoord2(ii));
0246
0247
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
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
0265 elemi_volt_nodes=volt_nodes(eleminodelist)';
0266 fem_sol=elemi_volt_nodes*phi_sol(:,ii);
0267
0268
0269 diff_sol=(fem_sol-analytic_sol)^2;
0270
0271
0272 diff_sol=diff_sol*weight2(ii);
0273
0274
0275 L2_error(i)=L2_error(i)+diff_sol;
0276 end
0277
0278
0279
0280
0281 for ii=1:size(weight1,2)
0282
0283 map=element_shape_function('tri3',xcoord1(ii),ycoord1(ii),zcoord1(ii));
0284
0285
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
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
0309 elemi_volt_nodes=volt_nodes(eleminodelist)';
0310 fem_sol=elemi_volt_nodes*(jacobianelem\dphi_sol(:,:,ii))';
0311
0312
0313 diff_sol=(fem_sol-analytic_sol)*(fem_sol-analytic_sol)';
0314
0315
0316 diff_sol=diff_sol*weight1(ii);
0317
0318
0319 H1semi_error(i)=H1semi_error(i)+diff_sol;
0320 end
0321
0322
0323 H1semi_error(i)=H1semi_error(i)*magjacelem;
0324 L2_error(i)=L2_error(i)*magjacelem;
0325
0326
0327 H1_error(i)=H1semi_error(i)+L2_error(i);
0328 end
0329
0330
0331
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