ex_fwd2d_high_order

 ensure dev/m_crabb/forward_problem is on the path

0001 % ensure dev/m_crabb/forward_problem is on the path
0002 
0003 %Make an inverse model and extract forward model
0004 imdl = mk_common_model('c2C0',16);
0005 fmdl = imdl.fwd_model;
0006 
0007 %Make an image and get voltages using eidors default solver
0008 img0 = mk_image(fmdl,1);
0009 
0010 %Forward solver using eidors
0011 v0=fwd_solve(img0); v0e=v0.meas;
0012 
0013 %Change the solvers
0014 fmdl.solve = @mc_fwd_solve;
0015 fmdl.system_mat = @mc_calc_system_mat;
0016 fmdl.jacobian   = @mc_calc_jacobian;
0017 
0018 %Add solver to perform p_refinement
0019 fmdl.fem_modify = @mc_fem_modify;
0020 
0021 %Linear, quadratic and cubic
0022 fmdl.mc_type    = 'tri3'; % linear case
0023 img1 = mk_image(fmdl,1);
0024 img1.fwd_solve.get_all_meas = 1; %Internal voltage
0025 v1 = fwd_solve(img1); v1e=v1.meas;
0026 
0027 %Quadratic FEM
0028 fmdl.mc_type    = 'tri6';
0029 img2 = mk_image(fmdl,1);
0030 img2.fwd_solve.get_all_meas = 1; %Internal voltage
0031 v2 = fwd_solve(img2); v2e=v2.meas;
0032 
0033 %Cubic FEM
0034 fmdl.mc_type    = 'tri10';
0035 img3 = mk_image(fmdl,1);
0036 v3 = fwd_solve(img3); v3e=v3.meas;
0037 
0038 %Electrode voltages
0039 figure; plot([v0e,v1e,v2e,v3e,[v0e-v1e,v2e-v0e,v3e-v0e]*100]);
0040 legend('0','1','2','3','1-0','2-0','3-0')
0041 xlim([1,100]);
0042 
0043 
0044 %Internal voltage for linear
0045 v1all = v1.volt; 
0046 img1n = rmfield(img1,'elem_data');
0047 img1n.node_data = v1all(1:size(fmdl.nodes,1),1);
0048 figure; show_fem(img1n,[1,0,0]);
0049 
0050 %Internal voltage or quadratic
0051 v2all = v2.volt; 
0052 img2n = rmfield(img2,'elem_data');
0053 img2n.node_data = v2all(1:size(fmdl.nodes,1),1);
0054 figure; show_fem(img2n,[1,0,0]);
0055 
0056 %Internal voltage difference of linear-quadratic
0057 img12n=img1n; img12n.node_data=v1all(1:size(fmdl.nodes,1),1)-v2all(1:size(fmdl.nodes,1),1);
0058 figure; show_fem(img12n,[1,0,0]);