test_3d_resistor

 Create 3D model of a Rectangular resistor
 $Id: test_3d_resistor.m 5865 2018-12-02 19:59:46Z aadler $

0001 % Create 3D model of a Rectangular resistor
0002 % $Id: test_3d_resistor.m 5865 2018-12-02 19:59:46Z aadler $
0003 
0004 ll=5*1; % length
0005 ww=1*2; % width
0006 hh=1*3; % height
0007 conduc= .13;  % conductivity in Ohm-meters
0008 current= 4;  % Amps
0009 z_contact= 1e-1;
0010 scale = .46;
0011 mdl= eidors_obj('fwd_model','3D rectangle');
0012 nn=0;
0013 for z=0:ll; for x=0:ww; for y=0:hh
0014    nn=nn+1;
0015    mdl.nodes(nn,:) = [x,y,z];
0016 end; end; end
0017 
0018 if 0
0019 % Matlab's delaunayn triangularization is screwed up
0020    mdl.elems = delaunayn(mdl.nodes);
0021 elseif 0
0022 % Matlab's delaunay3 triangularization is screwed up
0023    mdl.elems = delaunay3(mdl.nodes(:,1), mdl.nodes(:,2), mdl.nodes(:,3));
0024 elseif 0
0025    elem1= [ 1 2 3 5; 5 6 2 3; 5 6 7 3; ...
0026             4 2 3 8; 8 6 2 3; 8 6 7 3;];
0027    mdl.elems=[];
0028    for i=0:ll-1;
0029        mdl.elems= [mdl.elems; elem1+4*i];
0030    end
0031 else
0032    mdl= mk_grid_model([],0:ww,0:hh,0:ll);
0033    mdl.nodes= mdl.nodes*scale;
0034    mdl= rmfield(mdl,'coarse2fine');
0035 end
0036 mdl.boundary= find_boundary(mdl.elems);
0037 mdl.gnd_node = 1;
0038 elec_nodes= [1:(ww+1)*(hh+1)];
0039 elec(1).nodes= elec_nodes;      elec(1).z_contact= z_contact;
0040 elec(2).nodes= nn-elec_nodes+1; elec(2).z_contact= z_contact;
0041 stim.stim_pattern= [-1;1]*current;
0042 stim.meas_pattern= [-1,1];
0043 mdl.stimulation= stim;
0044 mdl.electrode= elec;
0045 show_fem(mdl);
0046 mdl = mdl_normalize(mdl,0);
0047 
0048 
0049 
0050 mdl.solve = @fwd_solve_1st_order;
0051 mdl.system_mat = @system_mat_1st_order;
0052 
0053 n_el = size(mdl.elems,1);
0054 img= eidors_obj('image','3D rectangle', ...
0055       'elem_data', ones(n_el,1) * conduc, ...
0056       'fwd_model', mdl); 
0057 
0058 fsol= fwd_solve(img);
0059 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0060 
0061 % analytical solution
0062 Block_R =  ll / ww / hh / scale/ conduc;
0063 Contact_R = z_contact/(ww*hh)/scale^2;
0064 % Contact R reflects z_contact / (width/scale)^2. Here we need to use
0065 %  the scale, since this is not reflected in the size of the
0066 %  FEM as created by the grid. This is different to the test_2d_resistor,
0067 %  where the FEM is created scaled, so that the ww
0068 %  don't need to be scaled by the scale parameter.
0069 R = Block_R + 2*Contact_R;
0070 
0071 V= current*R;
0072 fprintf('Solver %s: %f\n', 'analytic', V);
0073 fprintf('Solver %s: %f\n', 'analytic (no z_contact)', V - 2*Contact_R*current);
0074 
0075 if 0 % OLD CODE
0076    mdl.solve = @np_fwd_solve;
0077    mdl.jacobian = @np_calc_jacobian;
0078 else
0079    mdl.solve = @fwd_solve_1st_order;
0080    mdl.jacobian= @jacobian_adjoint;
0081 end
0082 mdl.misc.perm_sym= '{n}';
0083 mdl.normalize_measurements = 0;
0084 img= eidors_obj('image','3D rectangle', ...
0085       'elem_data', ones(size(mdl.elems,1),1) * conduc, ...
0086       'fwd_model', mdl); 
0087 
0088 fsol= fwd_solve(img);
0089 fprintf('Solver %s: %f\n', fsol.name, fsol.meas);
0090 
0091 
0092 % NOW CALCULATE THE ANALYTICAL JACOBIAN
0093 if 0 % OLD CODE
0094    mdl.solve = @np_fwd_solve;
0095    mdl.jacobian = @np_calc_jacobian;
0096 else
0097    mdl.solve = @fwd_solve_1st_order;
0098    mdl.jacobian= @jacobian_adjoint;
0099 end
0100 img.fwd_model = mdl;
0101 img.elem_data= ones(n_el,1) * conduc ;
0102 Jnp= calc_jacobian(img);
0103 
0104 mdl.jacobian = @jacobian_perturb;
0105 Jp1= calc_jacobian(img);
0106 
0107 img.elem_data= ones(n_el,1) * conduc ;
0108 Jp2= zeros(size(Jnp));
0109 delta= 1e-8;
0110 for i=1:n_el
0111    fsol_h= fwd_solve(img);
0112    img.elem_data(i) = img.elem_data(i) + delta;
0113    fsol_i= fwd_solve(img);
0114    Jp2(:,i) = (fsol_i.meas - fsol_h.meas)/delta; 
0115 end
0116    
0117 mdl.solve = @fwd_solve_1st_order;
0118 mdl.jacobian = @jacobian_adjoint;
0119 Jaa= calc_jacobian(img);
0120 
0121 fprintf('Jacobians: Cols by Jaa, Jnp, Jp1, Jp2\n');
0122 JJ = [Jaa;Jnp;Jp1;Jp2];
0123 disp(JJ(:,1:6))
0124 
0125 % UNIT_TEST
0126 unit_test_cmp('Jaa vs Jnp', Jaa, Jnp, 1e-6);
0127 unit_test_cmp('Jaa vs Jp1', Jaa, Jp1, 1e-6);
0128 unit_test_cmp('Jaa vs Jp2', Jaa, Jp2, 1e-3);
0129 
0130 
0131 function cleaner_code
0132 % Cleaner code using newer EIDORS helper functions (AA - dec 2018)
0133 ll=5; ww=1; hh=1; z_contact= 1e-1; conduc= .13; current= 4; scale = .46;
0134 mdl= rmfield( mk_grid_model([],0:ww,0:hh,0:ll), 'coarse2fine');
0135 mdl.boundary= find_boundary(mdl.elems); mdl.gnd_node = 1;
0136 mdl.electrode(1)=struct('z_contact', z_contact, 'nodes', ...
0137      find(mdl.nodes(:,3)==0));
0138 mdl.electrode(2)=struct('z_contact', z_contact, 'nodes', ...
0139      find(mdl.nodes(:,3)==ll));
0140 mdl.stimulation= stim_meas_list([1,2,1,2],2,current);
0141 mdl.nodes = mdl.nodes * scale;
0142 show_fem(mdl);
0143 Block_R =  ll / ww / hh / scale/ conduc;
0144 Contact_R = z_contact/(ww*hh)/scale^2;
0145 Vs = current * ( Block_R + 2*Contact_R );
0146 vh = fwd_solve(mk_image(mdl, conduc));
0147 unit_test_cmp('3D Resistor',Vs, vh.meas, 1e-10);